# XERUS:
An open-source tool for quick XRD phase identification and refinement
automation

Pedro Baptista de Castro<EMAIL_ADDRESS>National Institute for
Materials Science, 1-2-1 Sengen, Tsukuba, Japan  University of Tsukuba, 1-1-1
Tennoidai, Tsukuba, Japan  Kensei Terashima<EMAIL_ADDRESS>National Institute for Materials Science, 1-2-1 Sengen, Tsukuba, Japan  Miren
Garbine Esparza Echevarria  National Institute for Materials Science, 1-2-1
Sengen, Tsukuba, Japan  University of Tsukuba, 1-1-1 Tennoidai, Tsukuba, Japan
Hiroyuki Takeya  National Institute for Materials Science, 1-2-1 Sengen,
Tsukuba, Japan  Yoshihiko Takano  National Institute for Materials Science,
1-2-1 Sengen, Tsukuba, Japan  University of Tsukuba, 1-1-1 Tennoidai, Tsukuba,
Japan

###### Abstract

Analysis of XRD diffraction patterns is one of the keystones of materials
science and materials research. With the advancement of data-driven methods
for materials design, candidate materials can be quickly screened for the
study of a desiredphysical property. Efficient methods to automatically
analyze and identify phases present in a given pattern, are paramount for the
success of this new paradigm. To aid this process, the open source python
package Xray Estimation and Refinement Using Similarity (XERUS) for semi-
automatic/automatic phase identification is presented. XERUS takes advantages
of open crystal structure databases, not relying on proprietary databases, to
obtain crystal structures on the fly, being then chemical space agnostic. By
wrapping around GSASII, it can automatically simulate patterns and calculate
similarity measures used for phase identification. Our approach is simple and
quick but also applicable to multiphase identification, by coupling the
similarity calculations with quick refinements followed by an iterative peak
removal process. XERUS is shown in action in four different experimental
datasets and also it is benchmarked against a recently proposed deep learning
method for a mixture dataset covering the Li-Mn-O-F chemical space. XERUS will
be freely available on https://www.github.com/pedrobcst/Xerus/

Python, XRD, Automatic, Crystal Structure

##  1 Introduction

The establishment of high through-put computational databases such as the
Materials Project (MP)[1], the Open Quantum Materials Database[2] (OQMD), and
AFLOW[3], ignited a new age in materials science, where the conventional
trial-and-error approach for materials design is replaced by a data-centric
approach. More recently, artificial intelligence-based methods are also
becoming a key driver in this new paradigm of materials design[4], allowing
not only the creation of models for physical property prediction with machine
learning tools, but also allowing the extraction and creation of experimental
datasets[5, 6].
These new tools, that allow for the quick screening of candidate materials for
synthesis and experimental evaluation, also require methods than can quickly
identify the synthesis success of the proposed candidate. This often involves
analyzing the Xray Diffraction (XRD) of the obtained materials, a process that
is currently one of the primary bottlenecks of data-driven materials search
and automation since it requires time consuming human analysis[7].
To solve this issue, several methods for automating XRD analysis has been
proposed, with most recent approaches focusing in the use of deep-learning
(DL) assisted models. For example, models for automatic crystal system and
space group classification[8, 9], single phase[10, 11] and possible
multiphase[12, 13] have been proposed. Here we highlight as an example, the
recent work of Syzmanski[13] et. al, where the entire process of simulation
and model building for multiphase classification of the Li-Mn-Ti-O-F chemical
space took approximately 20 hours of computational time, achieving an accuracy
of 92% for multiphase classification. However, these models[9, 10, 11, 12, 13]
usually rely on commercial databases such as the ICSD[14].Evenmore, they
require simulation of large synthetic datasets for model training which are
both time consuming and restricts the model to a certain chemical space. On
the other hand, candidate-based materials search often involves different
materials spaning a large chemical space, where sometimes even the optimal
synthesis route can be unknown. For example, in the search of new
magnetocaloric materials[15, 16, 17], both the work of Bocarsly[16] et. al and
Court[17] et. al, all the identified candidate materials belonged to different
chemical spaces, ie. Mn-Nb-S and Co-Sc-V. Moreover, in the recent work of
Xiong[18] et. al, by high through-put screening of the MP database for
efficient photocatalysts, 11 different materials were identified and
synthesized for experimental evaluation, all belonging to different chemical
spaces, such as Ca-Pb-O, Ba-Nb-Mn-O, Sr-In-O, Ca-In-O, etc. In these
scenarios, the use of DL assisted methods for XRD classification, will
drastically increase the preparation time, since large synthetic datasets
followed by model training would be necessary for each individual chemical
space, slowing down the whole data-driven search process.
Thus, there is a need for an automatic/semi-automatic tool, that allows quick
identification/suggestion of possible synthesis results, which can adapt on-
the-fly to the chemical space. In this work, we introduce an alternative
approach for phase identification, by combining similarity calculation of
X-ray patterns with refinements, which does not require any pre-trained model,
in the open-source python framework Xray Estimation and Refinement Using
Similarity (XERUS) as an attempt to solve this problem.
XERUS takes advantage of open crystal structure databases, not relying on
proprietary data sources, enabling it to adapt to the synthesis chemical space
on-the-fly. In this manner, XERUS can suggest pre-refined multi/single phase
results and provides an interface for optimizing one-shot Rietveld analysis.
It was designed to be used either in Jupyter notebooks[19], as an interactive
tool with researchers, or to be adapted so it can be plugged-in directly into
experimental workflows, providing a quickly analysis of the results. In the
following sections, the core functions of XERUS and the method it uses to
attempt phase identification are described. We also present analysis results
for four different systems and we benchmark the results against the mixture
dataset of Ref. [13].

##  2 XERUS Description

Below we describe the main three main building blocks of XERUS.

  * •

Crystal Structure Management: Responsible for querying the MP, OQMD, AFLOW and
the Crystallographic Open Database (COD), parsing and caching the structures
for a given chemical space.

  * •

Backend: Responsible for pattern simulation, structure refinement,
optimization and visualization. Provides the interface for phase
identification.

  * •

Interface: Allows for ease of visualizing and analyzing the obtained results,
mostly through notebooks or by exporting the data.

These three building blocks are summarized in Figure 1. Figure 2 shows the
workflow of how XERUS works: from receiving the data, elements and directory
to save the results as input, going through phase identification, and finally
reporting the results to the user, while Figure 3 shows all the steps done by
XERUS in the phase identification process. All these steps will be discussed
in the sections below.

Figure 1: Overview of XERUS and its main building blocks, showing how crystal
structure is obtained, how the patterns are simulated and how the data can be
visualized. Figure 2: Workflow diagram of steps taken by XERUS after receiving
the data with given parameters such as elements and directory in which to
save. The preparation blocks present the steps from combining elements, to
querying structures, calculating patterns and a first calculation of
correlation, followed by structure filtering. This is then followed by the
phase identification loop described in Section 2.2 and, finally, the results
are reported to the user.

###  2.1 Preparation

To easily adapt to any given chemical space, XERUS takes advantage of the
computational databases API, such as AFLUX[20], the Materials Project API[21]
(MAPI), interfaced using pymatgen[22], the OQMD API and also the COD[23] REST
API, to obtain the relevant crystal structures for the space. Here, XERUS also
uses a local MongoDB[24] server, to locally cache the structures for future
uses.
Then, for a given list of elements, (eg. Ho, B and O), XERUS creates the
chemical space by first building all possible combinations of each element
[Figure 2 (a-1)], in this case Ho, B, O, Ho-B, Ho-O, B-O and Ho-B-O.
Afterwards, XERUS checks if each combination is already present in the local
database [Figure 2(a-2)]. If they are not found in the local database, XERUS
will automatically query the structures of the missing combination of the
input elements from the databases and obtain the structures in CIF format
[Figure 2 (a-3)], that will be used for simulation of patterns and refinement,
described in more detail the in Section 2.2. The acquired structures are then
parsed by pymatgen to obtain metadata such as composition, space group,
lattice system, and so on, and finally is they are saved into the the local
database. After all structures are obtained, they are finally saved into the
designed working directory [Figure 2 (a-4)].

###  2.2 Phase Identification

To achieve phase identification, XERUS implements a mixed approach of
correlation (similarity) calculations with quick Rietveld refinements,
followed by a pattern reduction, for the case of multiphase identification.
Correlation metrics have been widely studied for use in XRD phase
identification[7, 25, 26], being able to provide an accurate construction of
phase diagrams in combinatorial thin-film synthesis approach. Based on these
previous results, here we use the Pearson correlation metric $C_{ij}$, where i
and j are a given spectra with mean value $\bar{i}$ and $\bar{j}$, is defined
as[25]:

$C_{ij}=\sqrt{\frac{\sum_{k}(i_{k}-\bar{i})(j_{k}-\bar{j})}{\sum_{k}(i_{k}-\bar{i})^{2}\sum_{k}(j_{k}-\bar{j})^{2}}}$
(1)

Then by simulating the patterns of the obtained structures for a given
chemical space using GSASII[27, 28] [Figure 2 (a-5)], the Pearson correlation
between these patterns and the obtained experimental data can be quickly
calculated using the pandas package[29, 30] [Figure 2 (a-6)]. Note, that here
we will always be calculating the correlations between simulation and the
experimental data, therefore the correlation value between a simulated pattern
k and the experimental data will be noted as $C_{k}$.

After the calculation of $C_{k}$ for each simulated pattern, XERUS filter
repeated structures [Figure 2 (a-7)] by taking only one structure for the pair
of a given combination and spacegroup, that has the highest correlation value.
For example, for the structures in the ’Ho-B’ space with spacegroup P6/mmm,
only the structure from this group with highest $C_{k}$ will be used. This
step is necessary, since XERUS will obtain crystal structure from many
different sources, there will be possibly repeated crystal structures for the
same compound and as well off-stoichiometric variations of the compound.
Afterwards, the top g patterns with highest $C_{k}$, where g is a user-defined
parameter (default g= 3), are then refined against the experimental data using
GSAS II. Next, they are re-ranked by the lowest weighted profile residual
factor ($R_{wp}$) value, defined as[31]:

$R_{wp}=\sqrt{\frac{\sum_{i}w_{i}(y_{i}^{obs}-y_{i}^{calc})^{2}}{\sum_{i}w_{i}\left(y_{i}^{obs}\right)^{2}}}$
(2)

Where $w_{i}=\frac{1}{\sigma(y^{obs}_{i})^{2}}$ and $\sigma(y^{obs}_{i})$ is
the standard uncertainty[31]. This is an attempt into taking into account
artifacts that can lead to shift in peak positions such as lattice parameter
changes and zero shift, which are one of the weaknesses of a correlation only
approach[7, 26]. Here, we note that we do not perform an optimal refinement,
rather we just refine background function with fixed set of parameters, zero-
shift, lattice constants and instrument parameters, which is often enough to
give a reasonable first model for the data, in our experience. Then we define
the possible main phase of material as the one who achieved the lowest
$R_{wp}$ [Figure 2(b-1)] and then these g patterns are saved into the first
column of a refinement matrix $M_{g1}$. After defining the main phase, the
whole data (experimental and simulated) is reduced by removing boxes of width
$\delta$, where $\delta$ is a hyperparameter with default value of 1.3 degree
appropriate to our laboratory X-ray source, around the 2$\theta$ corresponding
to each Bragg position of the identified phase [Figure 2(b-2)]. This is
followed by a recalculation of $C_{k}$ for these reduced patterns, and we
again take the top g patterns as the next possible phases [Figure 2(b-3)].
Since the experimental data has been reduced, we do not perform any extra
refinement in this step. Instead, we reduce the data around the Bragg
positions of the highest correlated pattern in this run. These structures are
also saved in the next column of the refinement matrix, that is $M_{g2}$. This
continues until [Figure 2(b-4)]:

  * •

The maximum intensity of the experimental pattern drops below a threshold
(default 10%), this case is defined as the auto mode.

  * •

A number of runs n_runs has been elapsed, defined by the user.

This process is illustrated in Figure 3(a)-(f) for the case of an ideal three
phase search using default threshold of 10%. It is worth mentioning that, in
the case of multiphase identification, the top g structures of the last run
are also saved into the matrix.

Figure 3: Phase detection approach of XERUS. (a): Top g = 3 highest correlated
patterns with the experimental data, (b): Patterns re-ranked by $R_{wp}$ after
a quick Rietveld refinement has been performed, (c): The main phase is
identified and the associated bragg position peaks are removed by a box of
width $\delta$ degrees and these patterns are saved in the refinement matrix
$M_{gn}$, (d): After the removal, the correlations are recalculated and the
new g patterns are selected and saved into the next column of the refinement
matrix. (e): Now the highest correlated pattern is selected as the next
possible phase, and again its associated bragg positions are removed from the
data. (f): After a certain number of runs or the experimental data intensity
has dropped below a threshold, XERUS stops this process. In the case of
multiphase detection, these patterns are also considered as possible secondary
phases. (g): The patterns in the refinement matrix are combined and refined
and ranked by lowest $R_{wp}$. Here, we assume that the lowest $R_{wp}$
combination is the answer for the problem.

Now we turn attention to the refinement matrix $M_{gn}$. It‘s purpose is to
allow the combination of the identified main phase, which by construction lies
in the $M_{11}$ position, with the possible secondary phases, that is the
remaining g(n_runs \- 1) phases in the matrix to be finally refined as
exemplified in Figure 4(a). To do this, ideally one could just then refine the
identified main phase with all possible combinations of the secondary phases.
This would entail, a total of $\frac{k!}{(k-w)!w!}$, where $k=g(n\\_runs-1)$
and $w=n\\_runs-1$, phases for refinement, leading to strong combinatorial
explosion as g and n_runs increases. To soften this explosion, we define a
distance-based method on how the patterns can be combined in the matrix:

  * •

We always start from $M_{11}$, that is phase A in the example matrix of Figure
3

  * •

We allow the first pattern to start making combinations from the first or
second column.

  * •

Runs cannot be skipped. For instance, in Figure 4(b) case, the combination
A-E-F is not allowed.

These rules are illustrated in Figure 4(b), and in Figure 4(c) we show the
comparison of the required number of refinements for a simple combinatorial
method and the distance-based method. It is evident then, that the distance-
based approach effectively suppresses the combinatorial explosion even if
XERUS has to consider the presence of multiple phases, enabling their quick
identification. Here, we stress the benefit of having several runs while using
our similarity based approach. When a sample has multiple phases with minor
fractions in addition to the main phase, the one with highest intensity, those
minor phases usually will not have high similarity score in the first run,
where the main phase is still present. To account for possible extra phases in
the XRD pattern, XERUS reduces the experimental and simulated pattern, with
the expectation that the similarity between the remaining XRD pattern might
show higher score with the reduced simulated pattern for the minor phases in
the sample, in the next runs. Indeed, it is quite often that the identified
secondary phases, or minor phases, are ranked lower than g in the first run,
as it will be shown in the experimental results in Section 3.1.2 (phase D that
appears within g = 3, only in the second run, and it was identified as the
secondary phase) and in Section 3.1.3 (phase F also only appears within g at
the second run).
Here, we also note the case of how XERUS deal with overestimation of number of
phases, that is when the inputted n_runs is higher than the actual number of
phases in the sample. For a given n_runs, XERUS considers all possible number
of total present phases, including lower number of n_runs. Namely, suppose
n_runs = 3, the $R_{wp}$ not only for three-phase refinements but also for
single and two-phase are listed, so XERUS ranks all possible cases of phase
combinations within a given n_runs. In this manner, XERUS avoid the
overestimation of number of phases. Finally, all the considered phases and
their $R_{wp}$ are listed in the result output, where the highest rank would
be considered the the result in the automatic case, while also giving the
option for a final human check, if needed.

Figure 4:  (a) The refinement matrix definition. The possible main phase is
highlighted by the blue circle, while the possible secondary phases are in
red. (b) The overview of the distance-based approach for combination of the
phases. (c) A comparison between this method (solid lines) and simple
combination (dashed lines) of all possible phases for different values of g.

###  2.3 Optimization (Optional)

The main purpose of XERUS is to find out what are the secondary phases other
than the objective material allowing a quick identification of the synthesis
success. Depending on the percentage of such phases, adapting the synthesis
method to minimize such phases, if possible, would be ideal. However, here we
also acknowledge that phase determination is also not the only issue when
dealing with XRD data. Often, it is desired to obtain structural information,
such as lattice parameters, atomic positions, phase fractions (in the case of
multiphase patterns), and so on. This involves fitting structural models such
as Rietveld, which usually requires a initial guess of starting parameters
followed by sequential refinements where the parameters to be refined are
varied[32], to achieve optimal results. This is often done by a human expert,
using a trial and error approach guided by their own experience, which can be
time consuming.
To solve this issue, Ozaki et al. proposed the use of an blackbox optimization
approach (BBO)[33]. In this setup, $R_{wp}(\mathbf{X})$ is treated as the
function to optimize where $\mathbf{X}$ refers to a given input parameters
fora refinement.
XERUS adopts their method, and adapts it to allow for multiphase refinement
case. It can directly use obtained pre-refined patterns after the phase
detection part, without the need of downloading or preparing the search space,
allowing for the obtained results to be quickly optimized [Figure 2 (c-1)]
using the optuna package[34]. Here we note that, as discussed before, the
quick Rietveld refinements used for phase identification are not optimal, and
they are just used to help the phase identification process, and thus need to
be optimized. For more details in this approach, we refer the readers to Ref.
[33], where are in-depth discussion is provided.

###  2.4 User Interface

XERUS provides an low-code, easy-to-use API for exploring the results of phase
identification and optimization (if done). All he results are saved into a
pandas DataFrame, ordered by lowest $R_{wp}$ as discussed in Section 2.2. It
automatically parses GSASII files and provides out-of-the box static or
interactive plots allowing for a quick inspection of the obtained results. On
top of that, it can start the refinement optimization method described in
Section 2.3 of any of the obtained combinations [Figure 2 (c-2)].
XERUS also automatically exports almost all of the data used for phase
identification such as phase identification results, pattern simulations, CIF
structures, and optimization results, which can all be readily visualized with
third-party software if necessary.
The usage of these methods are explained in the example notebook that will be
available at
https://www.github.com/pedrobcst/Xerus/Examples/

##  3 Results and Discussion

To show how XERUS works, we will divide results into experimental data of
multiphase materials into two sections. In the first section, we will describe
the User-guided process where the n_runs (number of phases to search) is given
as an input. Here we will show the role of keeping the human in-the-loop for
the final decision on the results. In the second part, we will show the
straight results from XERUS in automatic mode, where the final result is taken
as the one that achieved the lowest $R_{wp}$ as discussed before. This will be
done in the Li-Mn-Ti-O-F mixture dataset of Ref. [13], where in the end we
will benchmark our results against their deep learning approach trained
specifically for this chemical space. Here we note that, for all cases, we
used the default values of the box width $\delta$ = 1.3, and the number of
patterns considered at each run, g = 3. The use of these parameters, in our
experience, is the one that give the best results in all the patterns we have
tested, and these parameters might require adjustment depending on the data.
We also note that, the compositions listed in this section are not an attempt
at identifying the synthesized material stoichiometry. Rather, they reflect
the composition parsed for a given structure contained at the database of a
given provider. They were also modified, when required, to the nominal
composition of a material for clarity. For example, in the case of the parsed
composition Mo8C4,we used Mo2C instead. The as-parsed compositions are
extensively listed in the Supplementary Information.

###  3.1 User-guided process

To illustrate the use cases of XERUS, we prepared four different samples:

  * •

A simple mixture of Ni and Cr powders

  * •

An alloy of arc-melted Ho and B with 1:2.1 mol proportion, where the target
material is HoB2. The extra amount of boron is to lead to formation of a HoB4
secondary phase.

  * •

A sample of La1-xSrxCuO4 ($x$ = 0.25), as an example of an oxide material.

  * •

An alloy of arc-melted Ho, Mo, C where the target compound is HoMoC2.

####  3.1.1 Ni-Cr powder mixture

The Ni and Cr powder mixture represents the easiest case of two-phase mixture
identification, since both Ni and Cr have non-overlapping peaks, and serves as
an example where the final answer is clearly known, since no reaction of the
elements takes place.
Figure 5 (a)-(b) shows the XERUS process applied for this mixture with n_runs
= 2. The phase cubic Fm$\bar{3}$m Ni phase (A) was identified as the possible
main phase, with the possible secondary phases: cubic Im$\bar{3}$m Cr (B), the
Cr-Ni Cr0.8Ni0.2 alloy (C) with same crystal structure as Cr, and the cubic
Fm$\bar{3}$m theoretical Cr structure[35] (D).

Figure 5: XERUS Results for a simple mixture of Cr-Ni. (a) Identification of
main phase by $R_{wp}$ ranking. Here A refers to the cubic Ni Fm$\bar{3}$m
phase, B is cubic Im$\bar{3}$m Cr phase and C is an alloy of cubic
Im$\bar{3}$m Cr0.8Ni0.2.(b) Second phase identification after the removing the
peaks of A. D refers to the theoretical cubic Fm$\bar{3}$m Cr phase. (c) Final
result after refinement of all combinations between identified phases where
the final shows the mixture of cubic Cr and Ni as expected.

Here we note that all the secondary phases have the same structure as either
Ni or Cr, with exception of the theoretical Cr phase,and thus have the same
pattern as the pure elements. Figure 5(c) shows the final combination chosen
as the answer: cubic Ni and cubic Cr as expected, which achieved the second
lowest $R_{wp}$. The lowest $R_{wp}$ combination was between Ni and the Cr-Ni
alloy. Since this alloy has the same structure as pure Cr, this combination
would also be acceptable as the answer. We note that the difference between
those two combinations in the $R_{wp}$ ranking was less than 1% (See Table SM.
1).

####  3.1.2 Ho-B alloy

To represent a more realistic experimental situation of a binary system, we
prepared a Ho-B alloy by arc melting. The Ho-B phase diagram[36] is simple and
provides a straightforward answer to which phases can exist. To attempt to
obtain a simple two-phase system, the target material was HoB2.1, which should
lead to a mixture of both HoB2 and HoB4 phases.

Figure 6: XERUS Results for the obtained alloy of Ho-B. (a) Identification of
main phase by $R_{wp}$ ranking. A: hexagonal P6/mmm HoB2 , B: cubic
Fm$\bar{3}$m Ho , C: trigonal R$\bar{3}$m Ho. (b) Second run of XERUS process.
D: tetragonal P4/mbm HoB4 and E: cubic Fm$\bar{3}$m HoB12. (c) Final result
after refinement of all combinations between A and the possible secondary
phases, showing that the obtained sample is a mixture of HoB2 and HoB4.

Figure 6 shows the XERUS process for this sample with n_runs = 2. Figure 6(a)
shows the first phase identification done by XERUS, where the main phase is
found as a hexagonal P6/mmm HoB2 (A). The other structures belong to different
metastable structures of Ho: cubic Fm$\bar{3}$m Ho (B) and trigonal
R$\bar{3}$m Ho (C). Figure 6(b) shows the new patterns after the peaks
belonging to the HoB2 phase have been removed. In this case, two new possible
phases appear: tetragonal P4/mbm HoB4 (D), with the highest correlation value,
and cubic Fm$\bar{3}$m HoB12 (E). The final result obtained by XERUS is shown
in Figure 6(c), showing that the sample is indeed a mixture of HoB2 and HoB4,
which is also the combination that achieved the lowest $R_{wp}$ (See Table
SM2).

####  3.1.3 LSCO

As an illustration of the application to an oxide ceramic material, we apply
XERUS with n_runs = 2, to a sample of polycrystalline LSCO synthesized by
solid-state reaction, with the larger chemical space of Sr-La-Cu-O.
Figure 7(a) shows the result for the main phase search. It is identified as
the tetragonal I4/mmm LSCO structure with composition La1.76Sr0.15CuO4 (A).
Here we also note that other phases are also both related to LSCO:
orthorhombic Cmce LSCO La1.87Sr0.13CuO4 (B) and the undoped tetragonal I4/mmm
La2CuO4. (C). The orthorhombic LSCO structure probably corresponds to the low-
temperature, as this material is known to undergo a tetragonal-orthorhombic
structural transition at lower temperatures[37, 38]. Figure 7(b) shows clearly
that the sample has at least a secondary phase. Interestingly, they are all
suggested to be La2O3, one of the starting materials used for sample
synthesis: the highest correlated pattern D trigonal P$\bar{3}$21 La2O3 (D),
hexagonal P63/mmc La2O3 (E) and trigonal P$\bar{3}$m1 La2O3 (F). Note that all
these structures possess extremely similar patterns. Indeed their existence is
due to diverging literature on the correct structure of La2O3 where diverse
authors reported possible structures for this material[39] and thus they all
exist at the COD database.

Figure 7: XERUS Results for the obtained LSCO sample. (a) Identification of
main phase by $R_{wp}$ ranking. A: tetragonal I4/mmm LSCO
(La1.85Sr0.15CuO3.6), B: orthorhombic Cmce LSCO (La1.87Sr0.13CuO4) and C:
tetragonal I4/mmm La2CuO4. (b) Second run of XERUS process. D: trigonal
P$\bar{3}$21 La2O3, E: hexagonal P63/mmc La2O3 and F: trigonal P$\bar{3}$m1
La2O3. (c) Final result after refinement of all combinations between A and the
possible secondary phases, showing that the obtained sample is a mixture of
tetragonal LSCO and trigonal La2O3.

Figure 7(c) shows the final result of the combination with lowest obtained
$R_{wp}$ (See Table SM3.). It shows a combination of I4/mmm tetragonal LSCO as
the main phase and P$\bar{3}$m1 trigonal La2O3 with a large wt%, an indication
that not all the starting materials have reacted. We note that in this
situation, since all the patterns of LSCO and La2O3 are almost identical,
XERUS could possibly find the answer as a combination of any of these
structures. Such a case highlights the importance of a final human
confirmation for an accurate determination of the phases, that is depicted in
Figure 7(c).

####  3.1.4 Ho-Mo-C alloy

As an example of a possible data-driven material, we attempted to synthesize a
sample of HoMoC2.
HoMoC2 is a material whose physical properties have never been reported. Only
its crystal structure is known[40], present in the MP database but not in the
COD. Moreover, the Ho-Mo-C phase diagram, to the best of our knowledge, has
never been reported in the literature, making the prediction of the possible
secondary phases that can appear in this system not straightforward. Figure 8
shows the XERUS process applied to this alloy with n_runs = 3. Figure 8(a)
shows the main phase identification of this alloy. Interestingly, it was
identified as the orthorhombic Pbcn $\alpha$-Mo2C (A). The other two proposed
phases are the orthorhombic Pca21 $\beta$-Mo2C (B) and high temperature
trigonal P$\bar{3}$1m Mo2C phase[41]. In the next run, depicted in Figure
8(b), the highest correlated pattern is found to be P63/mmc MoC (D), followed
by the target material orthorhombic Pnma HoMoC2 (E), and the hexagonal P63/mmc
Ho (F). In the final run, shown in Figure 8(c), two new candidates phases are
found, with the highested correlated with the remaining peaks being the
monoclinic P21/c Ho4C7 (G), and the last one being the I4/mmm HoC2 (H) The
final result, suggested by XERUS is illustrated in Figure 8(d), where the
lowest $R_{wp}$ combination shows that this alloy consists of at least of a
mixture of $\alpha$-Mo2C, HoMoC2 and HoC2 (See Table SM.4), indicating that
the synthesis of the target material was not successful. This case highlights
the importance of the hyperparameter g: both HoMoC2 and HoC2 were not the
highest correlated phase but they lied within the top g. Thus, if g was set to
a lower bound, such as 2, XERUS would miss the HoC2 phase. Also, we stress
here that such impurity phase identification and estimation of the phases
ratios have a massive importance for optimizing the synthesis
process/conditions.

Figure 8: XERUS Results for the obtained alloy of Ho-B. (a) Identification of
main phase by $R_{wp}$ ranking. A: orthorhombic Pbcn $\alpha$-Mo2C, B:
orthorhombic Pca21 $\beta$-Mo2C, C: trigonal P$\bar{3}$1m Mo2C. (b) Second run
of XERUS process. D: hexagonal P63/mmc MoC, E: orthorhombic Pnma HoMoC2, and
F: hexagonal P63/mmc Ho. (c) Third run search. G: monoclinic P21/c Ho4C7 and
H: tetragonal HoC2 I4/mmm (d) Final result after refinement of all
combinations, including two phase mixture and three phase mixture between A
and the possible secondary phases. The final result shows that the sample is
possibly a combination of $\alpha$-Mo2C, HoMoC2 and HoC2.

###  3.2 Automatic Process and Benchmarking

To evaluate better the accuracy of XERUS, a comparison between it and the
previously mentioned DL methods would be interesting. However, training the
models described by previous works would be extremely time consuming and in
most of the cases, they would need to be adapted to accept structures beyond
the ICSD. In addition, most of the experimental dataset of these previous
works are not openly available, with few exceptions. Therefore, here we
benchmark XERUS against the Li-Mn-Ti-O-F experimental mixture dataset of Ref.
[13] for the case of multiphase classification, and we compare the results
obtained with their optimized DL approach trained only for this chemical
space.
This dataset consists of 10 different mixtures of Li2TiO3, MnF2, LiMn2O4,
Li2MnO3, MnO, TiO2, LiMnO2, LiF and Mn2O3, with five three-phase and five
twophases mixtures, summarized in the first column of Table 1.
To compare with the DL method in this dataset, we ran XERUS in using n_runs =
auto with the default values of $\delta=1.3$, g = 3. In this case, a threshold
for the drop of the intensity maximum of the experimental data is necessary.
Here, we use a threshold of 10%, the default value of XERUS, in same manner as
Ref. [13].

The results obtained by XERUS are summarized in the second and third columns
of Table 1. Out of 25 possible phases, XERUS identifies correctly 22 phases
with 2 classification mistakes and one phase not found, achieving an accuracy
of 88% in autonomous mode, a result that is comparable to the DL approach that
scored a 92% accuracy for multiphase classification. To better understand the
possible weaknesses of XERUS an investigation on the source of the errors in
such a controlled dataset is important.
First, the inability of XERUS of finding one of the phases in the
Li2TiO3+MnF2+LiMn2O4 mixture highlights one of its weaknesses: highly
overlapping patterns. Figure 9(a) illustrates this case, where the simulated
patterns of each material in the mixture and the experimental data is shown.
Clearly there is a strong overlap between highest intensity peaks of Li2TiO3
and the other two structures. In such a case, since XERUS will reduce the
patterns based on the Bragg positions of the ”best” candidate of a given run,
thus the highest intensity peaks of the possible secondary phase will all be
removed and thus XERUS will fail to detect them. In this particular case,
LiMn2O4 is found as the main possible phase and the first secondary phase as
MnF2. After the second run, the intensity drops below the threshold and the
iteration stops and only two phases out of three is found. We note that, even
if XERUS is ran specifying the number of iterations, in this case n_runs = 3,
Li2TiO3 is not detected.

Figure 9: Identified pitfalls of XERUS. The black line with open dots refers
to the experimental data while the solid lines are the patterns automatically
simulated by XERUS. (a): A strong overlap between the peaks of the phases can
lead XERUS to fail to find one or more of the highly overlapping phases. (b)
and (c): Chemical compositions with different structures. (b) shows the room
temperature and low-temperature structure of LiMn2O4 and (c) shows diverging
structures reported in the literature regarding the crystal structure of
Li2MnO3.

Now lets turn our attention to the two misclassifications. The first
misclassification, where XERUS finds the tetragonal I41/amd LiMn2O4 phase is
similar to what was observed for the case of LSCO: The LiMn2O4 crystallizes
into a cubic Fd$\bar{3}$m structure at room temperature and upon cooling its
structure changes. Early reports indicated that the cubic structure is
distorted into a tetragonal I41/amd one[42]. However later single crystal
studies revealed that the structure is actually orthorhombic[43].
Nethertheless, similar to LSCO, both these structures exists in the COD. In
particular, the tetragonal distortion pattern is extremely similar to the room
temperature one, as shown in Figure 9(b), with exception of few low intensity
peaks, that will probably be lost in the background. In this situation, it is
natural for XERUS to suggest both phases as a solution. Indeed, the second
best solution found by XERUS is with the cubic structure and the $R_{wp}$
difference between these solutions is as low as 0.02% (See Table SM. 5).
The last misclassification is similar to the case described for La2O3: The
structure of LiMn2O3 was first reported to belong to the C2/c space group.
However, later it was confirmed to actually to have a C2/m symmetry[44].
Similar to La2O3, both of these structures are present in the COD and thus are
considered by XERUS. Figure 9(c) shows the simulated pattern of both
structures. Similar to the other cases, they are almost identical and the
second best combination is indeed with the C2/m structure, with a $R_{wp}$
difference of 0.3% between them (See Table SM. 6).
Here we note an important difference between the inability to find one othe
phases due to overlap and the two misclassifications due to similar
structures: In the case of overlap, the correct phase is never found even when
run manually with a given number of phases to search. In the other cases, the
correct answer usually lies in the top 3 suggested combinations, which can be
quickly checked in the end by a human. In this manner, we identify the major
weaknesses of XERUS as situations where there is strong phase overlap. In the
case of similar structures due to low-temperature phases and conflicting
literature data, all the phases can be easily distinguished by a researcher
with previous knowledge, and in this case XERUS only fails to identify one out
of 25 phases, leading to 96% accuracy. Therefore, XERUS can be exploited by
researchers to obtain a quick identification of their synthesis results,
allowing them to optimize the design of their synthesis method and conditions,
possibly accelerating materials development.

Table 1: Results of XERUS using automatic mode on the experimental data set of Ref. [13].The identified patterns are shown together with their space group. Here we take the result of XERUS as the combination that achieved the lowest $R_{wp}$. To classify as correct or incorrect, we take the answer as the stable structure at room temperature. Even if a structure with matching composition is found, if the structure found is a modification of the room temperature structure, for example low temperature, it is deemed incorrect. Note∗: The original mixture is described in the original data set as LiMn2O4+Li2MnO3+LiF. XERUS analysis indicates that this mixture might has been mislabelled, and LiMnO2 was mixed instead of Li2MnO3. Mixture | Correctly Identified | Mistaken [Not Found]
---|---|---
Li2TiO3+MnF2+LiMn2O4 | LiMn2O4(Fd$\bar{3}$m)+MnF2(P42/mnm) | [Li2TiO3 (C2/c)]
Li2MnO3+MnO+TiO2 | MnO(Fm$\bar{3}$m)+Li2MnO3(C2/m)+TiO2(I41/amd) | -
Mn2O3+TiO2+LiMnO2 | Mn2O3(Ia$\bar{3}$)+LiMnO2(Pmmn)+TiO2(I41/amd) | -
LiMn2O4+LiMnO2*+LiF | LiF(Fm$\bar{3}$m)+LiMnO2(Pmmn) | LiMn2O4(I41/amd)
Mn2O3+Li2TiO3 | Mn2O3(Ia$\bar{3}$)+Li2TiO3(C2/c) | -
Mn2O3+LiF | Mn2O3(Ia$\bar{3}$)+LiF(Fm$\bar{3}$m) | -
TiO2+MnO | MnO(Fm$\bar{3}$m)+TiO2(I41/amd) | -
MnF2+TiO2 | MnF2(P42/mnm)+TiO2(I41/amd) | -
LiMn2O4+Li2MnO3 | LiMn2O4(Fd$\bar{3}$m) | Li2MnO3 (C2/c)
LiMnO2+Li2MnO3 | Li2MnO3(C2/m)+LiMnO2(Pmmn) | -

##  4 Conclusion

In this work, we present the open source XERUS framework for semi-
automatic/automatic analysis of XRD patterns. By taking advantage of materials
simulations databases together with the COD, it can adapt on-the-fly to any
chemical space, obtaining the relevant crystal structures. We show that the
correlation method for pattern matching can be extended to multiphase
classification, by employing a simple iterative peak removal process coupled
with correlations calculations, where finally the found phases can be combined
and refined using GSAS II.
We apply this method first to four different synthesized samples, ranging from
a simple mixture to a complex alloy and show that XERUS can possibly identify
all the phases in these materials. We also benchmark XERUS against a Li-Mn-O-
Ti-F mixture dataset, and we show that XERUS running in automatic mode
achieves an accuracy almost as good as previous deep-learning methods in the
same dataset (88% vs 92%), without the need of simulation of large synthetic
datasets,training of large models or use of proprietary databases. With these
results, we also discuss the possible pitfalls of XERUS.
A tutorial with how to use XERUS, with the steps to reproduce all the results
obtained in this paper will be available at
https://www.github.com/pedrobcst/Xerus

##  5 Experimental Section

###  5.1 Sample Synthesis and Preparation

La2-xSrxCu$O_{4}$(LSCO) with $x=0.25$ was prepared by solid-state reaction.
The starting materials were La2O3, SrCO3(99%) and CuO($\geq 99.9\%)$ powders.
The La2O3 powder was prepared by dehydration of La[OH]3 in a furnace at 1000°C
overnight. The powders were mixed and made into pellets with 10 mm. of
diameter, afterwards placed on an alumina crucible and sintered in a furnace
for 8 hours at 950°C and then cooled slowly. The pellets were made into powder
and pelletized a second time, annealed again at 1000°C for 8 hours and cooled
slowly. The second sample, CrNi, was made by mixing commercial Cr and Ni
powders. The third sample, HoB2, was arc-melted from the elements in an Ar
atmosphere using Ho and B metals with a stoichiometry of 1:2.1. The fourth
sample, HoMoC2, was also arc-melted in an Ar atmosphere with stoichiometry
amounts Ho, Mo and C as starting materials.

###  5.2 X-ray analysis

Powder X-ray diffraction (XRD) analysis was done on a Rigaku MiniFlex 600
diffractometer using Cu-K$\alpha$ radiation, with 2$\theta$ ranging from 10°
to 70° and a 2$\theta$ scanning step of 0.02°.

##  6 Acknowledgements

We are grateful for fruitful discussions with Luca Foppiano. This work was
supported by the JST-Mirai Program (Grant No. JPMJMI18A3), JSPS KAKENHI (Grant
Nos. 19H02177, 20K05070, 20H05644, 20K22420), JSPS Bilateral Program (Grant
No. JPJSBP120214602) and JST CREST (Grant No. JPMJCR20Q4). P.B. Castro
acknowledges the scholarship support from the Ministry of Education, Culture,
Sports, Science and Technology (MEXT), Japan.

## References

  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. a. Persson, The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Kirklin  _et al._ [2015] S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, The open quantum materials database (oqmd): assessing the accuracy of dft formation energies, npj Computational Materials 1, 1 (2015).
  * Curtarolo  _et al._ [2012] S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demchenko, and D. Morgan, Aflow: An automatic framework for high-throughput materials discovery, Computational Materials Science 58, 218 (2012).
  * Balachandran [2019] P. V. Balachandran, Machine learning guided design of functional materials with targeted properties, Computational Materials Science 164, 82 (2019).
  * Olivetti  _et al._ [2020] E. A. Olivetti, J. M. Cole, E. Kim, O. Kononova, G. Ceder, T. Y.-J. Han, and A. M. Hiszpanski, Data-driven materials research enabled by natural language processing and information extraction, Applied Physics Reviews 7, 041317 (2020).
  * Foppiano  _et al._ [2021] L. Foppiano, S. Dieb, A. Suzuki, P. B. de Castro, S. Iwasaki, A. Uzuki, M. G. E. Echevarria, Y. Meng, K. Terashima, L. Romary, Y. Takano, and M. Ishii, Supermat: construction of a linked annotated dataset from superconductors-related publications, Science and Technology of Advanced Materials: Methods 1, 34 (2021).
  * Szymanski  _et al._ [2021a] N. J. Szymanski, Y. Zeng, H. Huo, C. J. Bartel, H. Kim, and G. Ceder, Toward autonomous design and synthesis of novel inorganic materials, Materials Horizons 8, 2169 (2021a).
  * Tiong  _et al._ [2020] L. C. O. Tiong, J. Kim, S. S. Han, and D. Kim, Identification of crystal symmetry from noisy diffraction patterns by a shape analysis and deep learning, npj Computational Materials 6, 1 (2020).
  * Oviedo  _et al._ [2019] F. Oviedo, Z. Ren, S. Sun, C. Settens, Z. Liu, N. T. P. Hartono, S. Ramasamy, B. L. DeCost, S. I. Tian, G. Romano, _et al._ , Fast and interpretable classification of small x-ray diffraction datasets using data augmentation and deep neural networks, npj Computational Materials 5, 1 (2019).
  * Wang  _et al._ [2020] H. Wang, Y. Xie, D. Li, H. Deng, Y. Zhao, M. Xin, and J. Lin, Rapid identification of x-ray diffraction patterns based on very limited data by interpretable convolutional neural networks, Journal of chemical information and modeling 60, 2004 (2020).
  * Lee  _et al._ [2020] J.-W. Lee, W. B. Park, J. H. Lee, S. P. Singh, and K.-S. Sohn, A deep-learning technique for phase identification in multiphase inorganic compounds using synthetic xrd powder patterns, Nature communications 11, 1 (2020).
  * Lee  _et al._ [2021] J.-W. Lee, W. B. Park, M. Kim, S. P. Singh, M. Pyo, and K.-S. Sohn, A data-driven xrd analysis protocol for phase identification and phase-fraction prediction of multiphase inorganic compounds, Inorganic Chemistry Frontiers 8, 2492 (2021).
  * Szymanski  _et al._ [2021b] N. J. Szymanski, C. J. Bartel, Y. Zeng, Q. Tu, and G. Ceder, Probabilistic deep learning approach to automate the interpretation of multi-phase diffraction spectra, Chemistry of Materials 33, 4204 (2021b), https://doi.org/10.1021/acs.chemmater.1c01071 .
  * Hellenbrandt [2004] M. Hellenbrandt, The inorganic crystal structure database (icsd)—present and future, Crystallography Reviews 10, 17 (2004).
  * de Castro  _et al._ [2020] P. B. de Castro, K. Terashima, T. D. Yamamoto, Z. Hou, S. Iwasaki, R. Matsumoto, S. Adachi, Y. Saito, P. Song, H. Takeya, and Y. Takano, Machine-learning-guided discovery of the gigantic magnetocaloric effect in hob sub 2 /sub near the hydrogen liquefaction temperature, NPG Asia Materials 12, 35 (2020).
  * Bocarsly  _et al._ [2017] J. D. Bocarsly, E. E. Levin, C. A. Garcia, K. Schwennicke, S. D. Wilson, and R. Seshadri, A simple computational proxy for screening magnetocaloric compounds, Chemistry of Materials 29, 1613 (2017).
  * Court  _et al._ [2021] C. J. Court, A. Jain, and J. M. Cole, Inverse design of materials that exhibit the magnetocaloric effect by text-mining of the scientific literature and generative deep learning, Chemistry of Materials 33, 7217 (2021).
  * Xiong  _et al._ [2021] Y. Xiong, Q. T. Campbell, J. Fanghanel, C. K. Badding, H. Wang, N. E. Kirchner-Hall, M. J. Theibault, I. Timrov, J. S. Mondschein, K. Seth, R. Katz, A. M. Villarino, B. Pamuk, M. E. Penrod, M. M. Khan, T. Rivera, N. C. Smith, X. Quintana, P. Orbe, C. J. Fennie, S. Asem-Hiablie, J. L. Young, T. G. Deutsch, M. Cococcioni, V. Gopalan, H. D. Abruña, R. E. Schaak, and I. Dabo, Optimizing accuracy and efficacy in data-driven materials discovery for the solar production of hydrogen, Energy & Environmental Science 14, 2335 (2021).
  * Kluyver  _et al._ [2016] T. Kluyver, B. Ragan-Kelley, F. Pérez, B. E. Granger, M. Bussonnier, J. Frederic, K. Kelley, J. B. Hamrick, J. Grout, S. Corlay, _et al._ ,  _Jupyter Notebooks-a publishing format for reproducible computational workflows._ , Vol. 2016 (2016).
  * Rose  _et al._ [2017] F. Rose, C. Toher, E. Gossett, C. Oses, M. B. Nardelli, M. Fornari, and S. Curtarolo, Aflux: The lux materials search api for the aflow data repositories, Computational Materials Science 137, 362 (2017).
  * Ong  _et al._ [2015] S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, The materials application programming interface api : A simple flexible and efficient api for materials data based on representational state transfer rest principles, Computational Materials Science 97, 209 (2015).
  * Ong  _et al._ [2013] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Python materials genomics pymatgen : A robust open-source python library for materials analysis, Computational Materials Science 68, 314 (2013).
  * Graulis  _et al._ [2009] S. Graulis, D. Chateigner, R. T. Downs, A. F. Yokochi, M. Q. s, L. Lutterotti, E. Manakova, J. Butkus, P. Moeck, and A. L. Bail, Crystallography open database - an open-access collection of crystal structures, Journal of Applied Crystallography 42, 726 (2009).
  * Inc. [2009] M. Inc., Mongodb, https://www.mongodb.com/ (2009).
  * Long  _et al._ [2007] C. Long, J. Hattrick-Simpers, M. Murakami, R. Srivastava, I. Takeuchi, V. L. Karen, and X. Li, Rapid structural mapping of ternary metallic alloy systems using the combinatorial approach and cluster analysis, Review of Scientific Instruments 78, 072217 (2007).
  * Iwasaki  _et al._ [2017] Y. Iwasaki, A. G. Kusne, and I. Takeuchi, Comparison of dissimilarity measures for cluster analysis of x-ray diffraction data from combinatorial libraries, npj Computational Materials 3, 4 (2017).
  * Toby and Dreele [2013] B. H. Toby and R. B. V. Dreele, Gsas-ii: The genesis of a modern open-source all purpose crystallography software package, Journal of Applied Crystallography 46, 544 (2013).
  * O’Donnell  _et al._ [2018] J. H. O’Donnell, R. B. Von Dreele, M. K. Chan, and B. H. Toby, A scripting interface for gsas-ii, Journal of Applied Crystallography 51, 1244 (2018).
  * pandas development team [2020] T. pandas development team, pandas-dev/pandas: Pandas (2020).
  * Wes McKinney [2010] Wes McKinney, Data Structures for Statistical Computing in Python, in  _Proceedings of the 9th Python in Science Conference_, edited by Stéfan van der Walt and Jarrod Millman (2010) pp. 56 – 61.
  * Toby [2006] B. H. Toby, R factors in rietveld analysis: How good is good enough?, Powder Diffraction 21, 67–70 (2006).
  * Aimi and Fujimoto [2020] A. Aimi and K. Fujimoto, Development of an automatic, high-throughput structural refinement method using rietveld analysis, ACS Combinatorial Science 22, 35 (2020), pMID: 31840982, https://doi.org/10.1021/acscombsci.9b00158 .
  * Ozaki  _et al._ [2020] Y. Ozaki, Y. Suzuki, T. Hawai, K. Saito, M. Onishi, and K. Ono, Automated crystal structure analysis based on blackbox optimisation, npj Computational Materials 6, 1 (2020).
  * Akiba  _et al._ [2019] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, Optuna: A next-generation hyperparameter optimization framework, in  _Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ (2019).
  * Häglund  _et al._ [1993] J. Häglund, A. Fernández Guillermet, G. Grimvall, and M. Körling, Theory of bonding in transition-metal carbides and nitrides, Phys. Rev. B 48, 11685 (1993).
  * Wei  _et al._ [2020] Q. Wei, M. Rong, S. Li, D. Su, Q. Yao, J. Wang, G. Rao, and H. Zhou, Thermodynamic assessment of the re-b (re=ho, er, tm) binary systems, Calphad 70, 101796 (2020).
  * Day  _et al._ [1987] P. Day, M. Rosseinsky, K. Prassides, W. David, O. Moze, and A. Soper, Temperature dependence of the crystal structure of the ceramic superconducting la1. 85sr0. 15cuo4: a powder neutron diffraction study, Journal of Physics C: Solid State Physics 20, L429 (1987).
  * Radaelli  _et al._ [1994] P. Radaelli, D. Hinks, A. Mitchell, B. Hunter, J. Wagner, B. Dabrowski, K. Vandervoort, H. Viswanathan, and J. Jorgensen, Structural and superconducting properties of la 2- x sr x cuo 4 as a function of sr content, Physical Review B 49, 4163 (1994).
  * Boldish and White [1979] S. I. Boldish and W. B. White, Vibrational spectra of crystals with the a-type rare earth oxide structure—i. la2o3 and nd2o3, Spectrochimica Acta Part A: Molecular Spectroscopy 35, 1235 (1979).
  * Jeitschko and Behrens [1986] W. Jeitschko and R. Behrens, Ternary carbides with ho 2 cr 2 c 3 and umoc 2 type structure, Zeitschrift fuer Metallkunde 77, 788 (1986).
  * Epicier  _et al._ [1988] T. Epicier, J. Dubois, C. Esnouf, G. Fantozzi, and P. Convert, Neutron powder diffraction studies of transition metal hemicarbides m2c1-x—ii. in situ high temperature study on w2c1-x and mo2c1-x, Acta Metallurgica 36, 1903 (1988).
  * Wills  _et al._ [1999] A. Wills, N. Raju, and J. Greedan, Low-temperature structure and magnetic properties of the spinel limn2o4: a frustrated antiferromagnet and cathode material, Chemistry of materials 11, 1510 (1999).
  * Akimoto  _et al._ [2000] J. Akimoto, Y. Takahashi, Y. Gotoh, and S. Mizuta, Single crystal x-ray diffraction study of the spinel-type limn2o4, Chemistry of Materials 12, 3246 (2000), https://doi.org/10.1021/cm0003673 .
  * Boulineau  _et al._ [2009] A. Boulineau, L. Croguennec, C. Delmas, and F. Weill, Reinvestigation of li2mno3 structure: Electron diffraction and high resolution tem, Chemistry of Materials 21, 4216 (2009), https://doi.org/10.1021/cm900998n .


KANAZAWA-21-09

KUNS-2889

Direct detection of pseudo-Nambu-Goldstone dark matter
with light mediator

Yoshihiko<EMAIL_ADDRESS>and Takashi
<EMAIL_ADDRESS>
1Department of Physics, Kyoto University, Kyoto 606-8502, Japan
2Institute of Liberal Arts and Science, Kanazawa University, Kanazawa,
920-1192 Japan
3Institute for Theoretical Physics, Kanazawa University, Kanazawa 920-1192,
Japan

##  1 Introduction

Dark matter direct detection experiments impose strong bounds on thermal dark
matter. The current strongest bound is given by the XENON1T/PandaX-4T [1, 2],
and will be further updated by the future XENONnT experiment [3]. The pseudo-
Nambu-Goldstone boson (pNGB) dark matter has been proposed as a candidate
naturally evading this strong constraint [4]. This is because the amplitude of
the elastic scattering with nuclei is suppressed by the small momentum of dark
matter through the derivative coupling. A non-zero elastic cross section is
induced at loop levels, however it is marginally small compared to future
sensitivities [5, 6, 7]. In addition, a sophisticated global fit of the model
with relevant constraints has also been performed [8].

For the pNGB dark matter model, an explicit soft breaking term is required to
give a mass for the NGB after the spontaneous breaking of the corresponding
global symmetry. It may be natural that the scale of the soft breaking term is
much smaller than the scale of the spontaneous symmetry breaking in the sense
that the global symmetry is approximate. However this is not necessarily true
if a ultra-violet (UV) completion of the pNGB model is considered [9, 10, 11,
12]. In this case, the effective pNGB model derived at low energy can have a
soft breaking term larger than the scale of the spontaneous symmetry breaking.
As a result, the mass of the particle mediating between dark matter and
Standard Model (SM) particles can be much lighter than the dark matter mass.

In this paper, we consider the pNGB dark matter with a light mediator. In this
scenario, we will show that non-zero elastic scattering cross section between
dark matter and nuclei emerges even though the interactions arise from the
derivative coupling. Furthermore, this scenario has a potential to
discriminate the pNGB dark matter and the other dark matter candidates such as
Weakly Interacting Massive Particles (WIMPs) by comparing the recoil energy
spectrum of the event rate at direct detection. We contemplate some relevant
constraints such as thermal production of dark matter, indirect detection
bounds and Higgs invisible decay. Then, we show that some parameter region has
already been excluded by the XENON1T bound, and some other region will be
explored by the future XENONnT experiment.

##  2 The model

We consider the model extended with a complex singlet scalar $S$. The model
has a global $U(1)_{S}$ symmetry under the transformation $S\to e^{i\alpha}S$,
which is broken to the $\mathbb{Z}_{2}$ parity by introducing the soft
breaking term $S^{2}$. The soft breaking term can be derived from a UV
completion of the model [9, 10, 11, 12]. The Lagrangian is given by

$\displaystyle\mathcal{L}=\left|\partial_{\mu}S\right|^{2}-\mathcal{V}(H,S),$
(2.1)

where the scalar potential $\mathcal{V}$ is written down as

$\mathcal{V}=-\frac{\mu_{H}^{2}}{2}|H|^{2}-\frac{\mu_{S}^{2}}{2}|S|^{2}+\frac{\lambda_{H}}{2}|H|^{4}+\lambda_{HS}|H|^{2}|S|^{2}+\frac{\lambda_{S}}{2}|S|^{4}-\frac{m_{\chi}^{2}}{4}\left(S^{2}+{S^{*}}^{2}\right).$
(2.2)

The $\lambda_{HS}$ term is the Higgs portal coupling giving the interaction
between the SM and $S$. The last term is the soft breaking term giving the
mass to the NGB after the spontaneous symmetry breaking. Then, the scalar
fields $H$ and $S$ can be parametrized as

$H=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ v+h\end{pmatrix},\quad
S=\frac{1}{\sqrt{2}}\left(v_{s}+s\right)e^{i\chi/v_{s}},$ (2.3)

where $v,v_{s}\in\mathbb{R}$ are the vacuum expectation values for $H$ and
$S$, respectively. Because of the interests in the low energy dynamics of pNGB
and the simplicity of calculation, we use the non-linear representation for
$S$ in this paper. These expectation values satisfy the following stationary
conditions:

$\displaystyle\mu_{H}^{2}=\lambda_{H}v^{2}+\lambda_{HS}v_{s}^{2},\quad\mu_{S}^{2}+m_{\chi}^{2}=\lambda_{S}v_{s}^{2}+\lambda_{HS}v^{2}.$
(2.4)

The CP-even components $h$ and $s$ mix with each other via the Higgs portal
coupling, and the mass matrix is given by

$\displaystyle\begin{pmatrix}\lambda_{H}v^{2}&\lambda_{HS}vv_{s}\\\
\lambda_{HS}vv_{s}&\lambda_{S}v_{s}^{2}\end{pmatrix}.$ (2.5)

This mass matrix is diagonalized by the unitary matrix, and thus the gauge
eigenstates $h$ and $s$ can be rewritten by the mass eigenstates $h_{1}$ and
$h_{2}$ as

$\displaystyle\begin{pmatrix}h\\\
s\end{pmatrix}=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}h_{1}\\\
h_{2}\end{pmatrix},$ (2.6)

where $h_{1}$ is identified as the SM-like Higgs boson with the mass
$m_{h_{1}}=125~{}\mathrm{GeV}$ and $h_{2}$ is the second Higgs boson whose
mass is assumed to be much lighter than $m_{h_{1}}$. The CP-odd component
$\chi$ is the pNGB with the mass $m_{\chi}$ and can be a stable dark matter
candidate due to a $\mathbb{Z}_{2}$ parity associated with the CP symmetry in
the dark sector. The dark matter mass is given by the soft breaking parameter
$m_{\chi}$.

The quartic couplings $\lambda_{H}$, $\lambda_{HS}$ and $\lambda_{S}$ in the
scalar potential can be rewritten in terms of the physical quantities as

$\displaystyle\lambda_{H}$ $\displaystyle=\frac{\cos^{2}\theta
m_{h_{1}}^{2}+\sin^{2}\theta m_{h_{2}}^{2}}{v^{2}},$ (2.7)
$\displaystyle\lambda_{HS}$
$\displaystyle=\frac{\sin\theta\cos\theta\left(m_{h_{2}}^{2}-m_{h_{1}}^{2}\right)}{vv_{s}},$
(2.8) $\displaystyle\lambda_{S}$ $\displaystyle=\frac{\sin^{2}\theta
m_{h_{1}}^{2}+\cos^{2}\theta m_{h_{2}}^{2}}{v_{s}^{2}},$ (2.9)

and the cubic couplings which are relevant to the subsequent sections are
given by

$\displaystyle\kappa_{111}$
$\displaystyle=3m_{h_{1}}^{2}\left(-\frac{\sin^{3}\theta}{v_{s}}+\frac{\cos^{3}\theta}{v}\right),$
(2.10) $\displaystyle\kappa_{112}$
$\displaystyle=\left(2m_{h_{1}}^{2}+m_{h_{2}}^{2}\right)\sin\theta\cos\theta\left(\frac{\sin\theta}{v_{s}}+\frac{\cos\theta}{v}\right),$
(2.11) $\displaystyle\kappa_{122}$
$\displaystyle=\left(m_{h_{1}}^{2}+2m_{h_{2}}^{2}\right)\sin\theta\cos\theta\left(-\frac{\cos\theta}{v_{s}}+\frac{\sin\theta}{v}\right),$
(2.12) $\displaystyle\kappa_{222}$
$\displaystyle=3m_{h_{2}}^{2}\left(\frac{\cos^{3}\theta}{v_{s}}+\frac{\sin^{3}\theta}{v}\right),$
(2.13)

with the convention

$\displaystyle\mathcal{V}\supset\frac{\kappa_{111}}{3!}h_{1}^{3}+\frac{\kappa_{112}}{2!}h_{1}^{2}h_{2}+\frac{\kappa_{122}}{2!}h_{1}h_{2}^{2}+\frac{\kappa_{222}}{3!}h_{2}^{3}.$
(2.14)

Note that the interactions between the dark matter and the SM particles come
from the kinetic term $|\partial_{\mu}S|^{2}$ and the soft breaking term
$-\frac{m_{\chi}^{2}}{4}(S^{2}+{S^{*}}^{2})$ since the field $S$ is written in
non-linear representation as Eq. (2.3).

It could be natural from the ’t Hooft sense that the soft breaking parameter
$m_{\chi}$ is much smaller than the scale of the spontaneous symmetry breaking
$v_{s}$ because the $\mathbb{Z}_{2}$ parity is enhanced to the global
$U(1)_{S}$ symmetry in $m_{\chi}\to 0$ limit. However note that $m_{\chi}$ and
$v_{s}$ are irrelevant and can be regarded as independent parameters once a UV
completion of the model is considered [9, 11, 10, 12]. Therefore, in the
following we consider the case of $v_{s},m_{h_{2}}\ll m_{\chi}$. From Eq.
(2.4) and $m_{h_{2}}^{2}\approx\lambda_{S}v_{s}^{2}$, it is found that the
quadratic mass parameter $\mu_{S}^{2}$ should be negative and $m_{\chi}^{2}$
is the trigger of the symmetry breaking in this case. Although the breaking
pattern is different, the interactions among the scalars do not change from
the ordinal pNGB dark matter model due to $v_{s}\in\mathbb{R}$ [4]. However,
in the light mediator case ($m_{h_{2}}<m_{\chi}$), we will show that non-zero
elastic cross section with nuclei is induced and can be tested by the current
and future direct detection experiments.

In particular, we focus on the mediator mass range of
$20~{}\mathrm{MeV}\lesssim m_{h_{2}}\lesssim 200~{}\mathrm{MeV}$. This is
because the lighter mass region has already been ruled out by the constraint
of big bang nucleosyntheis (BBN) [13] while the heavier mass region is
insensitive to dark matter direct detection as will be seen below. In the mass
range we focus on, the mixing angle $\sin\theta$ is roughly constrained in the
range of $2\times 10^{-5}\lesssim\sin\theta\lesssim 10^{-4}$ [14]. The upper
and lower bounds come from the constraints of the meson decays and BBN,
respectively.

##  3 The constraints

Before proceeding to the main subject, we review the constraints relevant to
our dark matter scenario.

###  3.1 Relic abundance

In the current scenario, the dark matter annihilation channel $\chi\chi\to
h_{2}h_{2}$ is dominant in the most of the parameter space. The cross section
is calculated from the four diagrams shown in Fig. 1 as

$\displaystyle\sigma_{h_{2}h_{2}}v_{\mathrm{rel}}\approx\frac{\lambda_{S}^{2}}{16\pi
s}\left|1+\frac{m_{h_{1}}^{2}}{m_{h_{2}}^{2}}\frac{\sin^{2}\theta
s}{s-m_{h_{1}}^{2}+im_{h_{1}}\Gamma_{h_{1}}}\right|^{2},$ (3.1)

where $s$ denotes the Mandelstam variable. In the above, $\sin\theta\ll 1$,
$m_{h_{2}}\ll m_{h_{1}},m_{\chi}$ and
$m_{h_{2}}^{2}\approx\lambda_{S}v_{s}^{2}$ are used to simplify the equation.

Figure 1: Diagrams for the main dark matter annihilation $\chi\chi\to
h_{2}h_{2}$.

For additional channels $\chi\chi\to WW,ZZ,f\overline{f}$, the cross sections
are calculated as

$\displaystyle\sigma_{WW}v_{\mathrm{rel}}$
$\displaystyle\approx\frac{\lambda_{S}g_{2}^{2}m_{W}^{2}}{8\pi}\sqrt{1-\frac{4m_{W}^{2}}{s}}\frac{\sin^{2}\theta}{sm_{h_{2}}^{2}}\left[3-\frac{s}{m_{W}^{2}}+\frac{s^{2}}{4m_{W}^{4}}\right]\frac{m_{h_{1}}^{4}}{(s-m_{h_{1}}^{2})^{2}+m_{h_{1}}^{2}\Gamma_{h_{1}}^{2}},$
(3.2) $\displaystyle\sigma_{ZZ}v_{\mathrm{rel}}$
$\displaystyle\approx\frac{\lambda_{S}g_{2}^{2}m_{Z}^{2}}{16\pi\cos^{2}\theta_{W}}\sqrt{1-\frac{4m_{Z}^{2}}{s}}\frac{\sin^{2}\theta}{sm_{h_{2}}^{2}}\left[3-\frac{s}{m_{Z}^{2}}+\frac{s^{2}}{4m_{Z}^{4}}\right]\frac{m_{h_{1}}^{4}}{(s-m_{h_{1}}^{2})^{2}+m_{h_{1}}^{2}\Gamma_{h_{1}}^{2}},$
(3.3) $\displaystyle\sigma_{f\overline{f}}v_{\mathrm{rel}}$
$\displaystyle\approx\frac{\lambda_{S}g_{2}^{2}m_{f}^{2}}{32\pi}\left(1-\frac{4m_{f}^{2}}{s}\right)^{3/2}\frac{\sin^{2}\theta}{m_{W}^{2}m_{h_{2}}^{2}}\frac{m_{h_{1}}^{4}}{(s-m_{h_{1}}^{2})^{2}+m_{h_{1}}^{2}\Gamma_{h_{1}}^{2}},$
(3.4)

where $g_{2}$ is the gauge coupling of the $SU(2)_{L}$ in the SM and
$f\overline{f}$ represents a pair of the SM fermions. Although the other
channels $\chi\chi\to h_{1}h_{1},h_{1}h_{2}$ are also possible, the
expressions are complicated and not shown here. The full expressions of the
annihilation cross sections are given in Appendix A. Note that all the cross
sections except for the main channel $\chi\chi\to h_{2}h_{2}$ are suppressed
by the small mixing angle $\sin^{2}\theta$. Thus it can be seen that the
channel $\chi\chi\to h_{2}h_{2}$ tends to be dominant. The annihilation cross
section for $\chi\chi\to h_{2}h_{2}$ in Eq. (3.1) is proportional to
$\lambda_{S}^{2}$ while those for the other channels in Eqs. (3.2)–(3.4) are
proportional to $\lambda_{S}$.

Then, the thermal averaged cross section can be given by [15]

$\displaystyle\langle\sigma{v}_{\mathrm{rel}}\rangle=\frac{1}{16m_{\chi}^{4}TK_{2}^{2}(m_{\chi}/T)}\int_{4m_{\chi}^{2}}^{\infty}\left(\sigma
v_{\mathrm{rel}}\right)\sqrt{s-4m_{\chi}^{2}}sK_{1}\left(\frac{\sqrt{s}}{T}\right)ds,$
(3.5)

where $T$ is the temperature of the universe and $K_{n}(z)~{}(n=1,2)$ denotes
the second kind modified Bessel function. The Boltzmann equation is
numerically solved by micrOMEGAs [16] with the above analytic formulas for the
annihilation cross sections, and the relic abundance should accommodate the
PLANCK observation $\Omega_{\chi}h^{2}\approx 0.12$ [17].111In the analytic
calculations of the cross sections, a large cancellation between the
corresponding diagrams occurs due to smallness of $v_{s}$. This causes huge
numerical errors if the cross sections are automatically evaluated by using
public codes such as CalcHEP [18]. Thus the above analytic results are used to
evaluate the thermal relic abundance of dark matter in micrOMEGAs [16].  Note
that the early kinetic decoupling effect may change the parameter space which
can reproduce the observed relic abundance when the dark matter mass is close
to the SM-like Higgs resonance [19].

###  3.2 Indirect detection

The main annihilation channel $\chi\chi\to h_{2}h_{2}$ generates cosmic rays
in the galaxy via subsequent $h_{2}$ decays. In the mass range we focus on
($30~{}\mathrm{MeV}\lesssim m_{h_{2}}\lesssim 200~{}\mathrm{MeV}$), the main
decay channel is $h_{2}\to e^{+}e^{-}$, which is constrained by the
$e^{+}e^{-}$ observation of AMS-02 [20]. In addition to the $e^{+}e^{-}$ flux,
gamma rays are also generated via final state radiation from the produced
$e^{+}e^{-}$ and bounded by the observations from dwarf spheroidal galaxies at
Fermi-LAT [21]. However this constraint is weaker than the $e^{+}e^{-}$
production [22].222The other annihilation channels into the SM particles also
produce gamma rays. However the constraint is weaker than that of $e^{+}e^{-}$
due to the small cross sections. The Cosmic Microwave Background (CMB) may
also be distorted by the dark matter annihilation because the produced charged
particles and gamma rays ionize the universe after the recombination era. Thus
the CMB measurement by PLANCK also sets a bound on the model.

We adopt the model independent bounds for multi-step cascade decays derived in
the literature [22] (top left panel in Fig. 11 therein). The extracted bounds
are shown in Fig. 2. As can be seen, the AMS-02 bound is severe for
$m_{\chi}\lesssim 100~{}\mathrm{GeV}$. For $m_{\chi}\gtrsim
300~{}\mathrm{GeV}$, the CMB bound can be stronger than the AMS-02 bound.
However the bound requires
$\langle\sigma{v}_{\mathrm{rel}}\rangle\lesssim\mathcal{O}(10^{-25})~{}\mathrm{cm^{3}/s}$,
which is consistent with thermal relic abundance.

Another comment is that the light mediator $h_{2}$ does not induce a long
range interaction. Thus we do not need to care about the non-perturbative
Sommerfeld effects for the dark matter annihilation [23, 24, 25, 26, 27]. This
is due to the nature of the pNGB dark matter that all the interactions are
written by the derivatives couplings.

Figure 2: Bounds on the annihilation cross section for $\chi\chi\to
h_{2}h_{2}$ from indirect detection extracted from the literature [22].

###  3.3 Higgs decay

Since the mediator $h_{2}$ is light enough, the SM-like Higgs boson $h_{1}$
decays into $h_{2}h_{2}$ whose decay width is calculated as

$\displaystyle\Gamma_{h_{1}\to h_{2}h_{2}}=\frac{\kappa_{122}^{2}}{32\pi
m_{h_{1}}}\sqrt{1-4\frac{m_{h_{2}}^{2}}{m_{h_{1}}^{2}}}\approx\frac{\lambda_{S}m_{h_{1}}^{3}\sin^{2}\theta}{32\pi
m_{h_{2}}^{2}},$ (3.6)

where $\sin\theta\ll 1$ and $m_{h_{2}}\ll m_{h_{1}}$ are assumed and
$m_{h_{2}}^{2}\approx\lambda_{S}v_{s}^{2}$ is used. This decay channel can be
regarded as an invisible decay because the lifetime of $h_{2}$ is long enough
to escape the detector of colliders. In addition to this decay channel,
another channel $h_{1}\to\chi\chi$ is also possible if the dark matter mass is
$m_{\chi}<m_{h_{1}}/2\approx 62.5~{}\mathrm{GeV}$. The decay width is
evaluated as

$\displaystyle\Gamma_{h_{1}\to\chi\chi}=\frac{m_{h_{1}}^{3}\sin^{2}\theta}{32\pi
v_{s}^{2}}\sqrt{1-4\frac{m_{\chi}^{2}}{m_{h_{1}}^{2}}}\approx\frac{\lambda_{S}m_{h_{1}}^{3}\sin^{2}\theta}{32\pi
m_{h_{2}}^{2}}\sqrt{1-4\frac{m_{\chi}^{2}}{m_{h_{1}}^{2}}}.$ (3.7)

The observation of the Higgs signal strength at the LHC is translated into the
constraint of the Higgs invisible decay whose branching fraction should
satisfy $\mathrm{Br}_{\mathrm{inv}}\leq 0.11$ [28, 29].

##  4 Direct detection of pNGB dark matter

It is known that the amplitude for the elastic scattering between the pNGB
dark matter and a nucleus $\chi A\to\chi A$ vanishes in non-relativistic limit
of dark matter. This argument is based on the premise that the mass of the
mediator $h_{2}$ is much heavier than the momentum transfer
$q=\sqrt{2m_{A}E_{R}}$ [4]. However this is not the case in our scenario where
the mass of the mediator $h_{2}$ can be the same order or smaller than the
momentum transfer. Assuming $q^{2},m_{h_{2}}^{2}\ll m_{h_{1}}^{2}$ and
$\sin\theta\ll 1$, the differential cross section for the elastic scattering
with a nucleus $A$ is calculated as

$\displaystyle\frac{d\sigma_{A}}{dE_{R}}\approx\frac{m_{A}\kappa_{A}^{2}}{8\pi
m_{\chi}^{2}v_{\chi}^{2}}\frac{\sin^{2}\theta}{v^{2}v_{s}^{2}}\frac{t^{2}}{\left(t-m_{h_{2}}^{2}\right)^{2}}F^{2}(E_{R}),$
(4.1)

where $t=-q^{2}$ is the Mandelstam variable, $v_{\chi}$ is the velocity of
dark matter and $F(E_{R})$ is the nuclear Helm form factor which is
parametrized by [30]

$\displaystyle F(E_{R})=\frac{3j_{1}(qR)}{qR}e^{-q^{2}\tilde{s}^{2}/2},$ (4.2)

with the spherical Bessel function $j_{1}(qR)$,
$R\approx\sqrt{r^{2}-5\tilde{s}^{2}}~{}\mathrm{fm}$, $r\approx
1.2A^{1/3}~{}\mathrm{fm}$, ($A$ is the mass number of the target nucleus) and
$\tilde{s}\approx 1~{}\mathrm{fm}$. The Helm form factor is normalized as
$F(0)=1$. The coupling $\kappa_{A}$ in Eq. (4.1) is given by

$\displaystyle\kappa_{A}=\kappa_{p}Z+\kappa_{n}(A-Z)=\sum_{q}\Big{[}f_{q}^{p}m_{p}Z+f_{q}^{n}m_{n}(A-Z)\Big{]},$
(4.3)

where $A$ and $Z$ are the mass number and the atomic number of the nucleus,
respectively. The coefficients $f_{q}^{p}$ and $f_{q}^{n}$ are the scalar
quark form factors of a nucleon, which are chosen as the values in the
literature [16]

$\displaystyle f_{u}^{p}=0.0153,\quad f_{d}^{p}=0.0191,\quad
f_{s}^{p}=0.0447,$ $\displaystyle f_{u}^{n}=0.0110,\quad
f_{d}^{n}=0.0273,\quad f_{s}^{n}=0.0447,$ (4.4)

for the light quarks. The form factors of the heavy quarks can be written in
terms of those of the light quarks.

Direct detection experiments provide the limits on the elastic scattering
cross section with a proton at zero momentum transfer, and the current
strongest limit is given by XENON1T/PandaX-4T experiments [1, 2]. Note that
this limit cannot directly be applied to the current scenario because one
cannot simply take zero momentum transfer limit in the differential cross
section in Eq. (4.1) due to a strong dependence on the momentum transfer and
light mediator mass. However, the experimental limit can be translated into
the limit on the total event rate. Then, it is compared with the predicted
event rate at a given parameters in the model.333A similar approach for direct
detection with long range interactions due to a light mediator has been
discussed [31]. The differential event rate is given by

$\displaystyle\frac{dR}{dE_{R}}=\frac{\rho_{\odot}}{m_{\chi}}N_{T}\int_{v_{\chi}>v_{\mathrm{min}}}\frac{d\sigma_{A}}{dE_{R}}v_{\chi}f_{\odot}(\bm{v}_{\chi})d^{3}v_{\chi},$
(4.5)

where $\rho_{\odot}\approx 0.3~{}\mathrm{GeV/cm^{3}}$ is the local dark matter
density, $N_{T}$ is the number of target nucleus, $f_{\odot}(\bm{v}_{\chi})$
is the Maxwell-Boltzmann velocity distribution function at the solar system
[32], and $v_{\mathrm{min}}$ is the minimum velocity of dark matter at the
given recoil energy $E_{R}$

$\displaystyle
v_{\mathrm{min}}=\sqrt{\frac{m_{A}E_{R}}{2}}\frac{1}{\mu_{A\chi}},$ (4.6)

where $\mu_{A\chi}=m_{A}m_{\chi}/(m_{A}+m_{\chi})$ is the reduced mass between
dark matter and nucleus.

Figure 3: (Left): Limit on the total event rate translated from the XENON1T
experiment in the range of the recoil energy
$4.9~{}\mathrm{keV}<E_{R}<40.9~{}\mathrm{keV}$ and $1.0$ t$\times$yr exposure
[1]. The solid (dashed) line represent the limit with (without) the efficiency
of the XENON1T detector. (Right): Comparison of energy spectra of number of
events for pNGB dark matter and singlet scalar dark matter with $1.0$
t$\times$yr exposure where dark matter mass is fixed to be
$m_{\chi}=150~{}\mathrm{GeV}$. The other parameters are taken as
$\lambda_{S}=0.09$, $\sin\theta=3\times 10^{-5}$ and
$m_{h_{2}}=60~{}\mathrm{MeV}$ for the pNGB dark matter, and
$\lambda_{hS}=0.006$ for the singlet scalar dark matter [33]. The solid
(dashed) lines represent the spectra with (without) detector efficiency
(XENON1T) [1].

For the experimental limit, assuming no isospin violation [34], one can
parametrize the differential cross section between dark matter and a nucleus
as

$\displaystyle\frac{d\sigma_{A}^{\mathrm{exp}}}{dE_{R}}=\frac{m_{A}A^{2}}{2\mu_{p\chi}^{2}v_{\chi}^{2}}\sigma_{p}^{\mathrm{exp}}F^{2}(E_{R}),$
(4.7)

where $\sigma_{p}^{\mathrm{exp}}$ is the experimental upper limit on the total
elastic cross section with a proton. Using this parametrization, the
differential event rate with the experimental limit
$\sigma_{p}^{\mathrm{exp}}$ can be rewritten as

$\displaystyle\frac{dR^{\mathrm{exp}}}{dE_{R}}=\frac{\rho_{\odot}}{m_{\chi}}N_{T}\frac{m_{A}A^{2}}{2\mu_{p\chi}^{2}}\sigma_{p}^{\mathrm{exp}}F^{2}(E_{R})\int_{v_{\chi}>v_{\mathrm{min}}}\frac{f_{\odot}(\bm{v}_{\chi})}{v_{\chi}}d^{3}v_{\chi}.$
(4.8)

The left panel of Fig. 3 shows the upper bound on the total event rate
$R^{\mathrm{exp}}$ obtained from the XENON1T limit on the total elastic
scattering cross section $\sigma_{p}^{\mathrm{exp}}$. The interest of region
for the recoil energy is $4.9~{}\mathrm{keV}<E_{R}<40.9~{}\mathrm{keV}$ and
the exposure $1.0~{}\mathrm{t}\times\mathrm{yr}$ is taken into account [1].

Figure 4: Parameter region in the ($m_{\chi}$, $\lambda_{S}$) plane for
$\sin\theta=10^{-4}$ and $m_{h_{2}}=30,60,100$ and $200~{}\mathrm{MeV}$. The
green, orange and violet region are excluded by XENON1T [1], AMS-02 [20] and
Higgs decays [28, 29], respectively. The red line corresponds to the parameter
space which can reproduce the observed relic abundance by PLANCK [17].

Figure 5: Same plots with Fig. 4 for $\sin\theta=3\times 10^{-5}$.

The recoil energy spectrum for the pNGB dark matter could be discriminative
from the other dark matter candidates thanks to its derivative couplings. The
typical energy spectrum for the pNGB dark matter is shown as violet lines in
the right panel of Fig. 3. The spectrum for the singlet scalar dark matter (an
example of usual WIMPs) is also shown as green lines for comparison [33]. The
event rate for the WIMP is enhanced at low energy while that for the pNGB dark
matter is suppressed as in the figure. Note that WIMPs may also be able to
induce energy spectra similar to the pNGB dark matter if the WIMP mass is
heavier and detector efficiency is taken into account because the spectra can
have a longer tail at higher energy for heavier WIMPs.444Similar spectra may
also be induced from the other dark matter scenarios such as WIMPs with
momentum dependent cross sections and inelastic scattering processes [32].
However, the shape of the recoil energy spectrum can be complementarily
utilized to discriminate the pNGB dark matter and the other candidates
combining with the other experiments and observations. Namely, for instance if
the dark matter mass is inferred in narrow range by the other experiments or
observations, one can use the recoil energy spectrum to discriminate the pNGB
dark matter from the usual WIMPs at the end.

Fig. 4 shows the parameter space in the $(m_{\chi},\lambda_{S})$ plane where
the Higgs mixing is fixed to be $\sin\theta=10^{-4}$ and the second Higgs mass
is $m_{h_{2}}=30,60,100$ and $200~{}\mathrm{MeV}$. Fig. 5 shows the same
figure with $\sin\theta=3\times 10^{-5}$. The green, orange and violet region
are excluded by the limit of XENON1T [1], AMS-02 [20] and Higgs decays [28,
29], respectively. The green dashed line represents the future sensitivity of
the XENONnT experiment [3]. The red line can reproduce the thermal relic
abundance consistent with the PLANCK observation $\Omega_{\chi}h^{2}\approx
0.12$ [17]. The constraints of direct detection and Higgs decay tend to be
severe for lighter $h_{2}$ and larger mixing $\sin\theta$ as can be seen from
Eqs. (3.6), (3.7) and (4.1). On the other hand, the AMS-02 bound does not much
depend on $m_{h_{2}}$ except for the region close to the Higgs resonance as
obvious from Eq. (3.1). As discussed in Ref. [35], the perturbative unitarity
constraint imposes the upper bound on the self quartic coupling as
$\lambda_{S}<8\pi/3$, however it is invariably weaker than the constraint from
Higgs decay in our scenario. From these plots, we can see that the thermally
produced pNGB dark matter with the light mediator can be consistent with all
the constraints when the dark matter mass is $m_{\chi}\sim 100~{}\mathrm{GeV}$
for $\sin\theta=3\times 10^{-5}$ and $m_{h_{2}}=60~{}\mathrm{MeV}$, and can be
tested by the future XENONnT experiment.

##  5 Conclusions

In the previous works, the pNGB dark matter has been completely insensitive to
dark matter direct detection experiments because the amplitude for the elastic
scattering vanishes in non-relativistic limit. In this paper, we have focused
on the case that the particle mediating the elastic scattering is light
enough, and we have found that a non-zero contribution to the amplitude
emerges. Together with the relevant constraints such as the thermal relic
abundance, $e^{+}e^{-}$, gamma ray, CMB and Higgs decays, we have shown the
parameter region of this scenario. Some parameter region have already been
excluded by the current XENON1T experiment and some other region can be tested
by the future XENONnT experiment. The pNGB dark matter can be discriminative
from the other dark matter candidates due to complementary study of dark
matter direct detection and cosmic ray observations.

## Acknowledgements

This work was supported by JSPS Grant-in-Aid for Scientific Research KAKENHI
Grant Nos. JP20J11901 (Y.A.), JP20K22349 (T.T.). Numerical computation in this
work was carried out at the Yukawa Institute Computer Facility.

##  Appendix A Annihilation cross sections

In this Appendix, we give analytic expressions for all the annihilation
channels for completeness. The annihilation cross sections for $\chi\chi\to
h_{i}h_{j}$ are given by

$\displaystyle\sigma_{h_{1}h_{1}}v_{\mathrm{rel}}$
$\displaystyle=\frac{1}{4\pi
s}\left(|a_{11}|^{2}+\frac{2}{v_{s}^{4}}\frac{\sin^{4}\theta
m_{h_{1}}^{8}}{c_{11}^{+}c_{11}^{-}}\right)\sqrt{1-\frac{4m_{h_{1}}^{2}}{s}}$
$\displaystyle\hskip 14.22636pt+\frac{\sin^{2}\theta m_{h_{1}}^{4}}{4\pi
s^{3/2}v_{s}^{2}\sqrt{s-4m_{\chi}^{2}}}\left(2\mathop{\mathrm{Re}}(a_{11})+\frac{\sin^{2}\theta}{v_{s}^{2}}\frac{m_{h_{1}}^{4}}{s-2m_{h_{1}}^{2}}\right)\log\left(\frac{c_{11}^{+}}{c_{11}^{-}}\right),$
(A.1) $\displaystyle\sigma_{h_{1}h_{2}}v_{\mathrm{rel}}$
$\displaystyle=\frac{1}{2\pi
s}\left(|a_{12}|^{2}+\frac{2}{v_{s}^{4}}\frac{\sin^{2}\theta\cos^{2}\theta
m_{h_{1}}^{4}m_{h_{2}}^{4}}{c_{12}^{+}c_{12}^{-}}\right)\frac{\sqrt{\lambda(s,m_{h_{1}}^{2},m_{h_{2}}^{2})}}{s}$
$\displaystyle\hskip 14.22636pt+\frac{\sin\theta\cos\theta
m_{h_{1}}^{2}m_{h_{2}}^{2}}{2\pi
s^{3/2}v_{s}^{2}\sqrt{s-4m_{\chi}^{2}}}\left(2\mathop{\mathrm{Re}}(a_{12})+\frac{\sin\theta\cos\theta}{v_{s}^{2}}\frac{m_{h_{1}}^{2}m_{h_{2}}^{2}}{s-m_{h_{1}}^{2}-m_{h_{2}}^{2}}\right)\log\left(\frac{c^{+}_{12}}{c^{-}_{12}}\right),$
(A.2) $\displaystyle\sigma_{h_{2}h_{2}}v_{\mathrm{rel}}$
$\displaystyle=\frac{1}{4\pi
s}\left(|a_{22}|^{2}+\frac{2}{v_{s}^{4}}\frac{\cos^{4}\theta
m_{h_{2}}^{8}}{c_{22}^{+}c_{22}^{-}}\right)\sqrt{1-\frac{4m_{h_{2}}^{2}}{s}}$
$\displaystyle\hskip 14.22636pt+\frac{\cos^{2}\theta m_{h_{2}}^{4}}{4\pi
s^{3/2}v_{s}^{2}\sqrt{s-4m_{\chi}^{2}}}\left(2\mathop{\mathrm{Re}}(a_{22})+\frac{\cos^{2}\theta}{v_{s}^{2}}\frac{m_{h_{2}}^{4}}{s-2m_{h_{2}}^{2}}\right)\log\left(\frac{c_{22}^{+}}{c_{22}^{-}}\right),$
(A.3)

where $s$ is the Mandelstam variable, the coefficients $a_{ij}$ and
$c_{ij}^{\pm}$ are given by

$\displaystyle a_{11}$
$\displaystyle=\frac{\sin^{2}\theta}{v_{s}^{2}}m_{h_{1}}^{2}+\frac{\sin\theta}{2v_{s}}\frac{\kappa_{111}s}{s-m_{h_{1}}^{2}+im_{h_{1}}\Gamma_{h_{1}}}-\frac{\cos\theta}{2v_{s}}\frac{\kappa_{112}s}{s-m_{h_{2}}^{2}+im_{h_{2}}\Gamma_{h_{2}}},$
(A.4) $\displaystyle a_{12}$
$\displaystyle=\frac{\sin\theta\cos\theta}{v_{s}^{2}}\frac{m_{h_{1}}^{2}+m_{h_{2}}^{2}}{2}-\frac{\sin\theta}{2v_{s}}\frac{\kappa_{112}s}{s-m_{h_{1}}^{2}+im_{h_{1}}\Gamma_{h_{1}}}+\frac{\cos\theta}{2v_{s}}\frac{\kappa_{122}s}{s-m_{h_{2}}^{2}+im_{h_{2}}\Gamma_{h_{2}}},$
(A.5) $\displaystyle a_{22}$
$\displaystyle=\frac{\cos^{2}\theta}{v_{s}^{2}}m_{h_{2}}^{2}+\frac{\sin\theta}{2v_{s}}\frac{\kappa_{122}s}{s-m_{h_{1}}^{2}+im_{h_{1}}\Gamma_{h_{1}}}-\frac{\cos\theta}{2v_{s}}\frac{\kappa_{222}s}{s-m_{h_{2}}^{2}+im_{h_{2}}\Gamma_{h_{2}}},$
(A.6) $\displaystyle c_{ij}^{\pm}$
$\displaystyle=s-m_{i}^{2}-m_{j}^{2}\pm\frac{\sqrt{s-4m_{\chi}^{2}}\sqrt{\lambda(s,m_{h_{i}}^{2},m_{h_{j}}^{2})}}{\sqrt{s}},$
(A.7)

with the kinematic function $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2yz-2zx$ and
the cubic couplings $\kappa_{ijk}~{}(i,j,k=1,2)$ are given in Eqs.
(2.10)–(2.13). $\Gamma_{h_{i}}$ denotes the total decay width of $h_{i}$.

For the other channels, the annihilation cross sections are given by

$\displaystyle\sigma_{WW}v_{\mathrm{rel}}$
$\displaystyle=\frac{g_{2}^{2}m_{W}^{2}}{8\pi}\sqrt{s-4m_{W}^{2}}\frac{s^{3/2}}{v_{s}^{2}}\cos^{2}\theta\sin^{2}\theta\left[3-\frac{s}{m_{W}^{2}}+\frac{s^{2}}{4m_{W}^{4}}\right]\left|\tilde{G}_{h_{1}}-\tilde{G}_{h_{2}}\right|^{2},$
(A.8) $\displaystyle\sigma_{ZZ}v_{\mathrm{rel}}$
$\displaystyle=\frac{g_{2}^{2}m_{Z}^{2}}{16\pi\cos^{2}\theta_{W}}\sqrt{s-4m_{Z}^{2}}\frac{s^{3/2}}{v_{s}^{2}}\cos^{2}\theta\sin^{2}\theta\left[3-\frac{s}{m_{Z}^{2}}+\frac{s^{2}}{4m_{Z}^{4}}\right]\left|\tilde{G}_{h_{1}}-\tilde{G}_{h_{2}}\right|^{2},$
(A.9) $\displaystyle\sigma_{f\overline{f}}v_{\mathrm{rel}}$
$\displaystyle=\frac{m_{f}^{2}}{8\pi}\sqrt{s}\left(s-4m_{f}^{2}\right)^{3/2}\frac{\cos^{2}\theta\sin^{2}\theta}{v^{2}v_{s}^{2}}\left|\tilde{G}_{h_{1}}-\tilde{G}_{h_{2}}\right|^{2},$
(A.10)

where $\tilde{G}_{h_{i}}$ denotes the $h_{i}$ propagator and
$\left|\tilde{G}_{h_{1}}-\tilde{G}_{h_{2}}\right|^{2}$ is given by

$\displaystyle\left|\tilde{G}_{h_{1}}-\tilde{G}_{h_{2}}\right|^{2}=\left|\frac{1}{s-m_{h_{1}}^{2}+im_{h_{1}}\Gamma_{h_{1}}}-\frac{1}{s-m_{h_{2}}^{2}+im_{h_{2}}\Gamma_{h_{2}}}\right|^{2}.$
(A.11)

Since the propagator part can be simplified as

$\displaystyle\left|\tilde{G}_{h_{1}}-\tilde{G}_{h_{2}}\right|^{2}\approx\frac{m_{h_{1}}^{4}}{s^{2}}\frac{1}{(s-m_{h_{1}}^{2})^{2}+m_{h_{1}}^{2}\Gamma_{h_{1}}^{2}},$
(A.12)

when $\sqrt{s},m_{h_{1}}\gg m_{h_{2}},\Gamma_{h_{1}},\Gamma_{h_{2}}$, one can
verify the above cross sections do not violate the unitarity at $s\to\infty$.

## References

  * [1] XENON Collaboration, E. Aprile et al., Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121 (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * [2] PandaX-4T Collaboration, Y. Meng et al., Dark Matter Search Results from the PandaX-4T Commissioning Run, arXiv:2107.13438 [hep-ex].
  * [3] XENON Collaboration, E. Aprile et al., Projected WIMP sensitivity of the XENONnT dark matter experiment, JCAP 11 (2020) 031, arXiv:2007.08796 [physics.ins-det].
  * [4] C. Gross, O. Lebedev, and T. Toma, Cancellation Mechanism for Dark-Matter–Nucleon Interaction, Phys. Rev. Lett. 119 (2017) 191801, arXiv:1708.02253 [hep-ph].
  * [5] D. Azevedo, M. Duch, B. Grzadkowski, D. Huang, M. Iglicki, and R. Santos, One-loop contribution to dark-matter-nucleon scattering in the pseudo-scalar dark matter model, JHEP 01 (2019) 138, arXiv:1810.06105 [hep-ph].
  * [6] K. Ishiwata and T. Toma, Probing pseudo Nambu-Goldstone boson dark matter at loop level, JHEP 12 (2018) 089, arXiv:1810.08139 [hep-ph].
  * [7] S. Glaus, M. Mühlleitner, J. Müller, S. Patel, T. Römer, and R. Santos, Electroweak Corrections in a Pseudo-Nambu Goldstone Dark Matter Model Revisited, JHEP 12 (2020) 034, arXiv:2008.12985 [hep-ph].
  * [8] C. Arina, A. Beniwal, C. Degrande, J. Heisig, and A. Scaffidi, Global fit of pseudo-Nambu-Goldstone Dark Matter, JHEP 04 (2020) 015, arXiv:1912.04008 [hep-ph].
  * [9] Y. Abe, T. Toma, and K. Tsumura, Pseudo-Nambu-Goldstone dark matter from gauged $U(1)_{B-L}$ symmetry, JHEP 05 (2020) 057, arXiv:2001.03954 [hep-ph].
  * [10] N. Okada, D. Raut, and Q. Shafi, Pseudo-Goldstone dark matter in a gauged $B-L$ extended standard model, Phys. Rev. D 103 (2021) 055024, arXiv:2001.05910 [hep-ph].
  * [11] Y. Abe, T. Toma, K. Tsumura, and N. Yamatsu, Pseudo-Nambu-Goldstone dark matter model inspired by grand unification, Phys. Rev. D 104 (2021) 035011, arXiv:2104.13523 [hep-ph].
  * [12] N. Okada, D. Raut, Q. Shafi, and A. Thapa, Pseudo-Goldstone Dark Matter in $SO(10)$, arXiv:2105.03419 [hep-ph].
  * [13] A. Fradette and M. Pospelov, BBN for the LHC: constraints on lifetimes of the Higgs portal scalars, Phys. Rev. D 96 (2017) 075033, arXiv:1706.01920 [hep-ph].
  * [14] M. W. Winkler, Decay and detection of a light scalar boson mixing with the Higgs boson, Phys. Rev. D 99 (2019) 015018, arXiv:1809.01876 [hep-ph].
  * [15] P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: Improved analysis, Nucl. Phys. B 360 (1991) 145–179.
  * [16] G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, micrOMEGAs5.0 : Freeze-in, Comput. Phys. Commun. 231 (2018) 173–186, arXiv:1801.03509 [hep-ph].
  * [17] Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [18] A. Belyaev, N. D. Christensen, and A. Pukhov, CalcHEP 3.4 for collider physics within and beyond the Standard Model, Comput. Phys. Commun. 184 (2013) 1729–1769, arXiv:1207.6082 [hep-ph].
  * [19] T. Abe, Early kinetic decoupling and a pseudo-Nambu-Goldstone dark matter model, arXiv:2106.01956 [hep-ph].
  * [20] AMS Collaboration, L. Accardo et al., High Statistics Measurement of the Positron Fraction in Primary Cosmic Rays of 0.5–500 GeV with the Alpha Magnetic Spectrometer on the International Space Station, Phys. Rev. Lett. 113 (2014) 121101.
  * [21] Fermi-LAT, DES Collaboration, A. Albert et al., Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT, Astrophys. J. 834 (2017) 110, arXiv:1611.03184 [astro-ph.HE].
  * [22] G. Elor, N. L. Rodd, T. R. Slatyer, and W. Xue, Model-Independent Indirect Detection Constraints on Hidden Sector Dark Matter, JCAP 06 (2016) 024, arXiv:1511.08787 [hep-ph].
  * [23] J. Hisano, S. Matsumoto, and M. M. Nojiri, Unitarity and higher order corrections in neutralino dark matter annihilation into two photons, Phys. Rev. D 67 (2003) 075014, arXiv:hep-ph/0212022.
  * [24] J. Hisano, S. Matsumoto, and M. M. Nojiri, Explosive dark matter annihilation, Phys. Rev. Lett. 92 (2004) 031303, arXiv:hep-ph/0307216.
  * [25] J. Hisano, S. Matsumoto, M. M. Nojiri, and O. Saito, Non-perturbative effect on dark matter annihilation and gamma ray signature from galactic center, Phys. Rev. D 71 (2005) 063528, arXiv:hep-ph/0412403.
  * [26] J. Hisano, S. Matsumoto, O. Saito, and M. Senami, Heavy wino-like neutralino dark matter annihilation into antiparticles, Phys. Rev. D 73 (2006) 055004, arXiv:hep-ph/0511118.
  * [27] J. Hisano, S. Matsumoto, M. Nagai, O. Saito, and M. Senami, Non-perturbative effect on thermal relic abundance of dark matter, Phys. Lett. B 646 (2007) 34–38, arXiv:hep-ph/0610249.
  * [28] CMS Collaboration, A. M. Sirunyan et al., Search for invisible decays of a Higgs boson produced through vector boson fusion in proton-proton collisions at $\sqrt{s}=$ 13 TeV, Phys. Lett. B 793 (2019) 520–551, arXiv:1809.05937 [hep-ex].
  * [29] ATLAS Collaboration, M. Aaboud et al., Combination of searches for invisible Higgs boson decays with the ATLAS experiment, Phys. Rev. Lett. 122 (2019) 231801, arXiv:1904.05105 [hep-ex].
  * [30] J. D. Lewin and P. F. Smith, Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil, Astropart. Phys. 6 (1996) 87–112.
  * [31] T. Hambye, M. H. G. Tytgat, J. Vandecasteele, and L. Vanderheyden, Dark matter direct detection is testing freeze-in, Phys. Rev. D 98 (2018) 075017, arXiv:1807.05022 [hep-ph].
  * [32] C. McCabe, The Astrophysical Uncertainties Of Dark Matter Direct Detection Experiments, Phys. Rev. D 82 (2010) 023530, arXiv:1005.0579 [hep-ph].
  * [33] GAMBIT Collaboration, P. Athron et al., Status of the scalar singlet dark matter model, Eur. Phys. J. C 77 (2017) 568, arXiv:1705.07931 [hep-ph].
  * [34] J. L. Feng, J. Kumar, D. Marfatia, and D. Sanford, Isospin-Violating Dark Matter, Phys. Lett. B 703 (2011) 124–127, arXiv:1102.4331 [hep-ph].
  * [35] C.-Y. Chen, S. Dawson, and I. M. Lewis, Exploring resonant di-Higgs boson production in the Higgs singlet model, Phys. Rev. D 91 (2015) 035015, arXiv:1410.5488 [hep-ph].