# Magnetic and electronic structure of the topological semimetal YbMnSb2

Jian-Rui Soh  Institute of Physics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), CH-1015 Lausanne, Switzerland  Siobhan M. Tobin  Department of
Physics, University of Oxford, Clarendon Laboratory, Oxford OX1 3PU, United
Kingdom  Hao Su  School of Physical Science and Technology, ShanghaiTech
University, Shanghai 201210, China  Ivica Zivkovic  Institute of Physics,
Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Bachir Ouladdiaf  Institut Laue-Langevin, 6 rue Jules Horowitz, BP 156, 38042
Grenoble Cedex 9, France  Anne Stunault  Institut Laue-Langevin, 6 rue Jules
Horowitz, BP 156, 38042 Grenoble Cedex 9, France  J. Alberto Rodríguez-
Velamazán  Institut Laue-Langevin, 6 rue Jules Horowitz, BP 156, 38042
Grenoble Cedex 9, France  Ketty Beauvois  Institut Laue-Langevin, 6 rue Jules
Horowitz, BP 156, 38042 Grenoble Cedex 9, France  Yanfeng Guo  School of
Physical Science and Technology, ShanghaiTech University, Shanghai 201210,
China  Andrew T. Boothroyd  Department of Physics, University of Oxford,
Clarendon Laboratory, Oxford OX1 3PU, United Kingdom

###### Abstract

The antiferromagnetic (AFM) semimetal YbMnSb2 has recently been identified as
a candidate topological material, driven by time-reversal symmetry breaking.
Depending on the ordered arrangement of Mn spins below the Néel temperature,
$T_{\mathrm{N}}$ = 345 K, the electronic bands near the Fermi energy can ether
have a Dirac node, a Weyl node or a nodal line. We have investigated the
ground state magnetic structure of YbMnSb2 using unpolarized and polarized
single crystal neutron diffraction. We find that the Mn moments lie along the
$c$ axis of the $P4/nmm$ space group and are arranged in a C-type AFM
structure, which implies the existence of gapped Dirac nodes near the Fermi
level. The results highlight how different magnetic structures can critically
affect the topological nature of fermions in semimetals.

PACS numbers

##  I Introduction

Topological semimetals are materials whose electronic bands have a linear
dispersion in the vicinity of the Fermi energy. Examples include crystals
where the valence and conduction bands meet at discrete points (as in Dirac or
Weyl semimetals) or along a one-dimensional curve in k-space (in the case of
nodal line semi-metals) [1, 2, 3]. These systems can host electrons that mimic
the behavior of massless fermions and are robust against perturbations due to
the protection afforded by the topology of the electronic band structure [4,
5].

The exploration of magnetic materials with topologically non-trivial
electronic band structure has become a central topic in quantum materials
physics. A particular interest is to identify materials in which the topology
of the electronic states can be controlled by their magnetic order [6, 7, 8,
9, 10]. Recently, the antiferromagnetic (AFM) metal YbMnSb2 has been proposed
as one such material, with the property that the topology of the bands near
the Fermi energy ($E_{\mathrm{F}}$) depends on the specific spin configuration
on the Mn magnetic sublattice [11, 12, 13, 14, 15, 16].

Previous studies have shown that the crystal structure of YbMnSb2 can be
described by the tetragonal space group $P4/nmm$ (no. 129) [11, 12, 13]. The
unit cell, shown in Fig. 1(a), includes a square sublattice of Sb atoms which
is predicted to host the topological fermions. This layer is in turn
sandwiched between layers of Mn, facilitating the coupling between magnetism
and electronic band topology. Depending on whether YbMnSb2 displays G-type,
C-type (with in-plane Mn moments), or canted C-type AFM order [see Figs.
1(b)–(d)], the electronic bands are predicted to give rise to a nodal-line
dispersion [13], a gapped Dirac crossing [12], or Weyl nodes [11],
respectively. It is essential, therefore, to determine the magnetic order of
the Mn sublattice in order to ascertain the topological nature of the fermions
in YbMnSb2.

Figure 1:  (a) Two unit cells of YbMnSb2 as described by the space group
$P4/nmm$ with cell parameters $a=b=4.31(2)$ Å, $c=10.85(1)$ Å at $T=400$ K.
The square Sb layers (at $z=0$, $c/2$, $c$) are predicted to host different
types of topological fermions, depending on the magnetic structure of the Mn
magnetic sublattice. (b)–(e) show different Mn magnetic structures and the
corresponding predicted topology of the electronic bands: (b) G-type AFM with
in-plane Mn moments stabilizes a line node in the electronic bands; (c) C-type
AFM with in-plane moments gives rise to Dirac nodes with a small gap; (d) a
canted C-type AFM order gives rise to topologically protected Weyl nodes; (e)
C-type AFM with moments along $c$ stabilizes a gapped Dirac state.

In this work we determined the spin configuration of the Mn magnetic
sublattice with high precision using a combination of polarized and
unpolarized single crystal neutron diffraction. We find that below the Néel
temperature the Mn spins point along the crystal $c$ axis and are arranged in
a C-type AFM structure [Fig. 1(e)]. The measurements place an upper limit of
0.01 $\mu_{\rm B}$ on any in-plane FM component (either above or below $T_{\rm
N}$). The particular C-type AFM structure found here is predicted to create
gapped Dirac points in the calculated band structure of YbMnSb2.

##  II Methods

Single crystalline samples of YbMnSb2 were grown by the flux method, as
described in Refs. [11, 13], giving rise to shiny rectangular platelets with
typical dimensions of $\sim\,4\times 4\times 0.5$ mm3. The structure and
quality of the single crystals were checked on a six-circle X-ray
diffractometer (Oxford Diffraction) and Laue diffractometer (Multiwire
Laboratories).

Magnetization measurements were performed on a superconducting quantum
interference device (SQUID) magnetometer (MPMS-3, Quantum Design) with
magnetic field applied parallel and perpendicular to the crystal $a$ axis.
Measurements were performed in the temperature range $T$ = 2 to 400 K and in
fields of up to 3 T. Electrical transport measurements were performed on a
Physical Properties Measurement System (PPMS, Quantum Design) with the
resistivity option, in the temperature range 2 to 400 K in zero field.

Single crystal neutron diffraction with unpolarized neutrons of wavelength
$\lambda=2.36$ Å was performed on the four-circle diffractometer D10 at the
Institut Laue–Langevin (ILL). The scattered neutrons were measured with a
$94\times 94$ mm2 area detector. To determine the magnetic and crystal
structure of YbMnSb2, a total of 382 and 386 $hkl$ reflections were collected
below and above the magnetic ordering temperatures, at $T=2$ K and 400 K,
respectively [Fig. 3(a)]. To ascertain if there is a canting of the Mn moments
away from the crystal $c$ axis [Fig. 4(a),(b)], the temperature dependence
several reflections were measured in the temperature range from 2 to 420 K.

Figure 2: (a) Temperature dependent magnetization curves of YbMnSb2 measured
with a field $\mu_{0}H=3$ T applied parallel and perpendicular to the crystal
$c$ axis. The bifurcation of the curves below $T_{\mathrm{N}}$ = 345(1) K
indicate the onset of AFM order. (b) Integrated intensity of the 100 magnetic
reflection measured by neutron diffraction. (c) Temperature dependence of the
in-plane resistivity of YbMnSb2, displaying metallic conductivity for
$T<T_{\mathrm{N}}$. Figure 3: (a) Reciprocal space representation of the $hkl$
reflections that were measured on the D10 instrument. The blue and red points
correspond to reflections which are allowed and forbidden, respectively, by
the $P4/nmm$ space group. (b) Comparison of measured ($I_{\mathrm{obs}}$) and
calculated ($I_{\mathrm{calc}}$) integrated intensities from the structural
refinement in the paramagnetic phase of YbMnSb2 at $T=400$ K. (c), (d) and (e)
compare $I_{\mathrm{obs}}$ and $I_{\mathrm{calc}}$ for the $T=2$ K data
refined against structural models with G-type AFM, C-AFM (moments in-plane)
and C-AFM (moments along $c$), respectively.

Polarized neutron diffraction measurements ($\lambda=0.832$ Å) were performed
on the single-crystal diffractometer D3 (ILL) with a polarized incident beam
but without polarization analysis (half-polarized setup). The incident neutron
polarization ($P_{\mathrm{i}}$) was oriented either parallel or antiparallel
to $z$ (the normal to the scattering plane) by means of a cryo-flipper, and
maintained by guide fields. The crystal was aligned inside a vertical field
superconducting magnet with the crystal $b$ axis vertical to give access a
series of reflections in the $h0l$ horizontal scattering plane [Fig. 5(b)]. A
magnetic field of $\mu_{0}H=0.4$ T was applied.

For both neutron diffraction experiments on the instruments D3 and D10, the
single crystals were pre-aligned with the neutron Laue diffractometers
OrientExpress [17] and CYCLOPS [18] at ILL. Neutron diffraction data analysis
was performed with the Mag2Pol software [19].

To clarify the topological nature of the electronic band structure in YbMnSb2
we performed density functional theory (DFT) calculations of the electronic
band structure using the Quantum Espresso implementation [20]. Relativistic
pseudo-potentials were used in the calculations to account for the large spin-
orbit coupling (SOC) arising from the heavy Sb ions which might lead to band
inversion. Furthermore, a Hubbard $U=4.1$ eV was used to model the strong
electron-electron correlations due to the Mn $3d$ bands. The Mn spin
configuration used in the calculations was obtained from the analysis of the
neutron scattering data as input parameters. A Monkhorst–Pack k-point sampling
mesh of $4\times 4\times 6$ was used [21].

##  III Results and Discussion

Figure 2(a) plots the magnetization of YbMnSb2 as a function of temperature in
a fixed field of $\mu_{0}H=$ 3 T, applied parallel and perpendicular to the
crystal $a$ axis. On cooling, we observe a bifurcation of the curves below
$T_{\rm N}$ = 345(1) K, indicating the onset of antiferromagnetic order of the
Mn spins (the Yb atoms are divalent and therefore non-magnetic). The
anisotropy of the induced moment in the magnetically ordered phase, with
$M_{\perp c}>M_{\|c}$, suggest that the crystal $c$ axis is an easy axis. At
low temperatures (below $T$ = 50 K), we observe an upturn in the magnetization
in both field configurations, consistent with an earlier report [11]. In a
later section, we will discuss the possible origin of this behavior.

Figure 2(b) displays the temperature dependence of the integrated intensity of
the 100 reflection measured by neutron diffraction. The 100 reflection is
structurally forbidden by the $n$ glide, and its intensity has an order
parameter-like dependence. Its appearance coincides with the bifurcation in
the magnetic susceptibility at $T_{\mathrm{N}}$. These observations confirm
the existence of antiferromagnetism in YbMnSb2.

In Fig. 2(c) we plot the in-plane longitudinal resistivity ($\rho_{xx}$) as a
function of temperature. On cooling, we observe that $\rho_{xx}$ increases to
a maximum at $T_{\mathrm{N}}$ then decreases with temperature down to $T$ = 2
K (the lowest temperature reached) in good agreement with earlier studies [12,
11, 13].

The broad resistivity peak in the vicinity of $T_{\mathrm{N}}$ could indicate
a coupling between charge transport and fluctuations of the Mn moments, which
diverge at $T_{\mathrm{N}}$. Ab inito calculations of the electronic band
structure of YbMnSb2 demonstrate that the Sb square layer, which is
structurally sandwiched between two Mn layers [Fig. 1(a)], is mainly
responsible for charge transport [13]. Coupling between Mn $3d$ and Sb $4p$
states makes it possible for different Mn magnetic structures to cause
different Sb $4p$ electronic topologies [see Figs. 1(b)–(d)], and this is why
it is important to understand the details of the AFM structure observed below
$T_{\rm N}$.

All of the magnetic structures depicted in Fig. 1 can be described with a
magnetic propagation vector $\textbf{q}=(0,0,0)$. For $T<T_{\rm N}$,
therefore, magnetic Bragg peaks observed by neutron diffraction may coincide
with structural Bragg peaks. Figure 3(a) illustrates the set of reflections
measured in our diffraction study. These reflections are mostly structurally
allowed in the $P4/nmm$ space group (indicated by the blue circles), but there
are also a small number of reflections of the form $h+k=2m+1$ ($m$ integer)
which are forbidden by the $n$ glide (red squares). The latter reflections, if
present, are therefore purely magnetic in origin.

To identify the magnetic contribution to the structural Bragg peaks when
$T<T_{\rm N}$ it is important to obtain a good model for the crystal structure
in the paramagnetic phase. Figure 3(b) shows the results of our refinement in
the paramagnetic phase at $T$ = 400 K. The observed and calculated peak
intensities are seen to agree well. The refined structural parameters are
given in the Supplemental Material [22].

For completeness, we mention that we also observed some very weak reflections
at half-integer $l$ positions which are structurally forbidden in $P4/nmm$ but
allowed in the closely related $Pnma$ space group. The latter space group was
reported to give an equally good description of x-ray diffraction data from
polycrystalline YbMnSb2 [11] and is known in the wider family of $A$MnSb2 ($A$
= Eu, Sr, Ca) compounds [23, 24, 25, 26, 27, 28]. These half-integer
reflections can arise from a subtle structural distortion which doubles the
unit cell in the out-of-plane direction, but due to limited data we have not
attempted to refine the structure further. The central conclusions of the
study, which concern the magnetic structure, are not affected by the existence
of such a distortion.

Figures 3(c)–(e) compare the measured ($I_{\mathrm{obs}}$) and predicted
($I_{\mathrm{calc}}$) integrated intensities based on the magnetic structures
shown in Figs. 1(b), (c) and (e). There are large discrepancies between
$I_{\mathrm{obs}}$ and $I_{\mathrm{calc}}$ for the G-type AFM model, Fig.
3(b), especially for the purely magnetic reflections which have zero intensity
for this magnetic structure. Hence, the G-type AFM order, shown in Fig. 1(b),
can safely be excluded. Similarly, the $I_{\mathrm{calc}}$ values for the in-
plane C-type AFM model significantly underestimate the intensities of the
purely magnetic reflections [Fig. 3(c)]. Therefore, the in-plane C-type AFM
model shown in Fig. 1(c) can also be excluded.

Figure 4: (a) Unpolarized neutrons were used on the D10 instrument (ILL), with
the scattering vector Q along $h00$. (b) Since magnetic neutron diffraction is
sensitive to the component of M perpendicular to Q, namely
$\textbf{M}_{\perp}$, a canting of the Mn moments away from the crystal $c$
axis will produce a magnetic contribution to the nuclear Bragg reflections
along $00l$. (c) The temperature dependence of the integrated intensities of
the 001 and the 002 reflections. The line corresponds to the fit to the Debye-
Waller factor.

On the other hand, there is good agreement [see Fig. 3(d)] between the
predicted and measured integrated intensities for the C-type AFM model with Mn
moments along the crystal $c$ axis, shown in Fig. 1(e). The size of the Mn
moments obtained from the refinement is 3.48(9) $\mu_{\mathrm{B}}$, which is
somewhat less than the value of 5 $\mu_{\mathrm{B}}$ expected for divalent
manganese with localized electrons ($3d^{5}$, $S=5/2$), but consistent with
the Mn ordered moments reported in related $A$Mn$X_{2}$ intermetallic
compounds ($A$ = Ca, Sr, Ba, Yb; $X$ = Sb, Bi) [26, 29, 30, 31, 32].

Let us now consider the structure shown in Fig. 1(d), which involves a canting
of the C-type AFM structure so as to produce a net in-plane ferromagnetic (FM)
moment. The possibility of a small FM canting was suggested by Wang et al.
[11] to account for the upturn in the magnetization at low temperatures, and
was reported in neutron scattering measurements on the sister compound SrMnSb2
[26]. In that study, Liu et al. reported an increase in the integrated
intensities of the $200$ reflection in the vicinity of $T_{\mathrm{N}}$, and
attributed the increase to an in-plane ferromagnetic order which sets in at
$T_{\mathrm{C}}\sim$ 565 K (well above $T_{\mathrm{N}}$). This in-plane
ferromagnetic order in SrMnSb2 is thought persist into the AFM ordered phase
below $T_{\mathrm{N}}$, giving rise to the canted AFM order shown in Fig.
1(d).

To search for a small in-plane FM component we studied the temperature
dependence of several $00l$ reflections. Magnetic neutron diffraction is
caused by the component of the magnetic moment perpendicular to the scattering
vector Q, so if the moments are aligned along the $c$ axis then there no
magnetic contribution to any $00l$ reflections. On the other hand, any canting
of the moments will produce a magnetic component perpendicular to $c$ which
will give magnetic scattering at $00l$ reflections. Figures 4(a) and (b)
illustrate the arrangement just described.

Figure 4(c) plots the temperature dependence of the integrated intensities of
the $001$ and $002$ reflections. At low temperatures (below $T$ = 50 K), we do
not observe an upturn in the integrated intensities which would indicate the
onset of in-plane ferromagnetic order. Therefore, the increase in magnetic
susceptibility at low temperatures, shown in Fig. 2(a), cannot be attributed
to the development of spontaneous ferromagnetic component in the crystal
$a$-$b$ plane. This conclusion is consistent with behavior of the
magnetization as a function of field at low $T$ which also does not indicate
spontaneous in-plane ferromagnetism [see Supplemental Materials [22]].
Instead, we attribute the upturn in the magnetic susceptibility at low
temperatures as due to a very small concentration of paramagnetic impurities,
as reported in other members of the $A$Mn$X_{2}$ family [33, 30]. Hence, the
canted AFM order, shown in Fig. 1(d), can be excluded at low temperatures.

However, we observed a gradual increase in the integrated intensity of the 001
and 002 reflections above $T\simeq$ 200 K [Fig. 4(c)]. This increase, which is
too large to be accounted for by paramagnetic diffuse scattering, could
potentially arise from an in-plane ferromagnetic order of the Mn moments,
giving rise to a canted magnetic structure [Fig. 1(d)] when $T<T_{\rm N}$, as
was reported in the neutron diffraction study of SrMnSb2 [26]. To test for the
presence of in-plane ferromagnetism in YbMnSb2 we used polarized neutron
diffraction to measure the flipping ratio ($R$) of a series of reflections
within the $h0l$ scattering plane at $T=400\,\textrm{K}>T_{\mathrm{N}}$ [see
Fig. 5(a) and (b)].

Figure 5: (a) To access the $h0l$ reflections in the half-polarized
experimental setup on D3, a YbMnSb2 single crystal was aligned with the
crystal $b$ axis vertical. The incident neutrons are polarized along $z$
($P_{\mathrm{i}}\parallel z$) and the scattered neutron polarization
$P_{\mathrm{f}}$ is not analyzed. (b) Reciprocal space plot of the measured
reflections (blue spheres) in the $h0l$ scattering plane. (c) Measured
flipping ratios (blue circles) along with the calculated $R$ for several in-
plane magnetic moments from 0 $\mu_{\mathrm{B}}$ (dashed line) to 0.03
$\mu_{\mathrm{B}}$.

For centrosymmetric crystals and a magnetic field that is perpendicular to the
scattering plane the flipping ratio is given by [34],

$R=\frac{F_{\mathrm{N}}^{2}+F_{\mathrm{M}}^{2}+2F_{\mathrm{N}}F_{\mathrm{M}}}{F_{\mathrm{N}}^{2}+F_{\mathrm{M}}^{2}-2F_{\mathrm{N}}F_{\mathrm{M}}},$
(1)

where $F_{\mathrm{M}}$ and $F_{\mathrm{N}}$ are the magnetic and nuclear
structure factors. When there is no magnetic contribution to the measured
reflections ($F_{\mathrm{M}}=0$) the flipping ratio is seen from Eq. (1) to be
$R=1$. Conversely, a flipping ratio that deviates from unity implies a non-
zero $F_{\mathrm{M}}$. As the applied field was small ($\mu_{0}H=0.4$ T), the
induced moment on paramagnetic Mn is negligible. Hence, if $R\neq 1$ we can
infer a small FM contribution.

Figure 5(c) plots the measured flipping ratios $R$ of the reflections that
were studied. The plot also includes the calculated $R$ for several values of
the in-plane ferromagnetic moment up to 0.03 $\mu_{\mathrm{B}}$. In
particular, the dashed horizontal line denotes $R=1$, which corresponds to
zero ferromagnetic moment.

It can be seen from Fig. 5(c) that the measured $R$ is consistent with an in-
plane ferromagnetic component of less than 0.01$\,\mu_{\mathrm{B}}$. This is
far too small to explain the observed increase in integrated intensity of the
$00l$ reflections above $T\sim$ 200 K, which we estimate would require an in-
plane moment of about 0.8 $\mu_{\mathrm{B}}$. Hence, the polarized neutron
measurements exclude spontaneous in-plane magnetic order above $T\sim$ 200 K.
Furthermore, the magnetization as a function of field at $T$ = 400 K also
indicates that any in-plane ferromagnetic component is negligible ($M\lesssim
0.001\,\mu_{\mathrm{B}}$) (see Supplemental Materials [22]). Therefore,
although it is unclear what causes the observed rise in the integrated
intensity of $00l$ reflections above $T\sim$ 200 K [Fig. 4(c)], we can safely
exclude a ferromagnetic origin.

The evidence presented here excludes all of the proposed magnetic structures
presented in Figs. 1(b)–(d) as acceptable models for the magnetic order of
YbMnSb2, and instead we find C-type AFM order with Mn moments along the
crystal $c$ axis [Fig. 1(e)]. We note that some of the predicted magnetic
structures in Figs. 1(b)–(d) were based on ab initio calculations, and the
difference between the self-consistent energies obtained for the various
magnetic structures can sometimes be of the same order of magnitude as the
computational uncertainty. Hence, our work highlights the importance of
identifying the ground state magnetic structure experimentally before
proceeding to calculate the electronic band structure.

Figure 6:  Ab initio electronic dispersion of YbMnSb2 along in-plane high
symmetry directions. The calculation assumes the magnetic structure determined
in this study.

Finally, Fig. 6 shows the in-plane electronic dispersion of YbMnSb2 calculated
with the magnetic structure determined in this study. The band structure
contains two avoided Dirac crossings close the Fermi energy, $\Delta
E\sim-0.2$ eV along $\Gamma$–$X$ high symmetry line and $\Delta E\sim 0.15$ eV
between the $M$ and $\Gamma$ point. These gaps arises from the hybridization
of the Sb $5p$ bands [14, 35]. The gapped Dirac dispersion of YbMnSb2 is
similar to that reported in other compounds exhibiting C-type AFM magnetic
order in the $A$MnSb2 family [27, 36, 37, 38].

##  IV Conclusion

In this study we tested experimentally several magnetic structures which have
been predicted for YbMnSb2. Magnetic structure refinements made from
unpolarized neutron diffraction data provide conclusive evidence for C-type
AFM order with Mn spins along the $c$ axis. Polarized neutron diffraction
measurements could not detect any in-plane FM order above or below $T_{\rm
N}$. The results do not support predictions that YbMnSb2 is either a nodal
line or a Weyl semimetal, and instead imply that the band structure contains
gapped Dirac nodes near the Fermi level. More generally, this work highlights
how different magnetic structures can profoundly affect the topological
behavior of electrons in metals.

###### Acknowledgements.

The authors wish to thank H. Rønnow for helpful discussions, and D.
Prabhakaran for assistance with sample preparation. The proposal numbers for
the neutron scattering experiments on YbMnSb2 are 5-41-1049 (D10), 5-53-305
(D3 [39]), EASY-671 (CYCLOPS) and EASY-667 (OrientExpress). Y-F.G. and A.T.B.
acknowledge support from the Oxford–ShanghaiTech collaboration project. S. M.
T. was supported by a scholarship from the Rhodes Trust. J.-R. S. acknowledges
support from the Singapore National Science Scholarship, Agency for Science
Technology and Research.

## References

  * Yang  _et al._ [2018] S.-Y. Yang, H. Yang, E. Derunova, S. S. P. Parkin, B. Yan, and M. N. Ali, Symmetry demanded topological nodal-line materials, Advances in Physics: X 3, 1414631 (2018).
  * Fang  _et al._ [2016] C. Fang, H. Weng, X. Dai, and Z. Fang, Topological nodal line semimetals, Chin. Phys. B 25, 117106 (2016).
  * Kumar  _et al._ [2021] N. Kumar, S. N. Guin, K. Manna, C. Shekhar, and C. Felser, Topological quantum materials from the viewpoint of chemistry, Chem. Rev. 121, 2780 (2021).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Nagaosa  _et al._ [2020] N. Nagaosa, T. Morimoto, and Y. Tokura, Transport, magnetic and optical properties of Weyl materials, Nat. Rev. Mater. 5, 015001 (2020).
  * Xu  _et al._ [2020] Y. Xu, L. Elcoro, Z.-D. Song, B. J. Wieder, M. G. Vergniory, N. Regnault, Y. Chen, C. Felser, and B. A. Bernevig, High-throughput calculations of magnetic topological materials, Nature 586, 702 (2020).
  * Zhang  _et al._ [2021] Z. Zhang, Z.-M. Yu, and S. A. Yang, Magnetic higher-order nodal lines, Phys. Rev. B 103, 115112 (2021).
  * Gui  _et al._ [2019] X. Gui, I. Pletikosic, H. Cao, H.-J. Tien, X. Xu, R. Zhong, G. Wang, T.-R. Chang, S. Jia, T. Valla, W. Xie, and R. J. Cava, A new magnetic topological quantum material candidate by design, ACS Cent. Sci. 5, 900 (2019).
  * Wang  _et al._ [2019] L.-L. Wang, N. H. Jo, B. Kuthanazhi, Y. Wu, R. J. McQueeney, A. Kaminski, and P. C. Canfield, Single pair of weyl fermions in the half-metallic semimetal $\mathrm{EuC}{\mathrm{d}}_{2}\mathrm{A}{\mathrm{s}}_{2}$, Phys. Rev. B 99, 245147 (2019).
  * Soh  _et al._ [2019a] J.-R. Soh, F. de Juan, M. G. Vergniory, N. B. M. Schröter, M. C. Rahn, D. Y. Yan, J. Jiang, M. Bristow, P. Reiss, J. N. Blandy, Y. F. Guo, Y. G. Shi, T. K. Kim, A. McCollam, S. H. Simon, Y. Chen, A. I. Coldea, and A. T. Boothroyd, Ideal Weyl semimetal induced by magnetic exchange, Phys. Rev. B 100, 201102(R) (2019a).
  * Wang  _et al._ [2018] Y.-Y. Wang, S. Xu, L.-L. Sun, and T.-L. Xia, Quantum oscillations and coherent interlayer transport in a new topological Dirac semimetal candidate YbMnSb2, Phys. Rev. Mater. 2, 021201(R) (2018).
  * Kealhofer  _et al._ [2018] R. Kealhofer, S. Jang, S. M. Griffin, C. John, K. A. Benavides, S. Doyle, T. Helm, P. J. W. Moll, J. B. Neaton, J. Y. Chan, J. D. Denlinger, and J. G. Analytis, Observation of a two-dimensional Fermi surface and Dirac dispersion in YbMnSb2, Phys. Rev. B 97, 045109 (2018).
  * Qiu  _et al._ [2019] Z. Qiu, C. Le, Z. Liao, B. Xu, R. Yang, J. Hu, Y. Dai, and X. Qiu, Observation of a topological nodal-line semimetal in YbMnSb2 through optical spectroscopy, Phys. Rev. B 100, 125136 (2019).
  * Klemenz  _et al._ [2019] S. Klemenz, S. Lei, and L. M. Schoop, Topological semimetals in square-net materials, Ann. Rev. Mater. Res. 49, 185 (2019).
  * Klemenz  _et al._ [2020] S. Klemenz, A. K. Hay, S. M. L. Teicher, A. Topp, J. Cano, and L. M. Schoop, The role of delocalized chemical bonding in square-net-based topological semimetals, J. Amer. Chem. Soc. 142, 6350 (2020), pMID: 32142261\.
  * Pan  _et al._ [2021] Y. Pan, F.-R. Fan, X. Hong, B. He, C. Le, W. Schnelle, Y. He, K. Imasato, H. Borrmann, C. Hess, B. Büchner, Y. Sun, C. Fu, G. J. Snyder, and C. Felser, Thermoelectric properties of novel semimetals: A case study of YbMnSb2, Adv. Mater. 33, 2003168 (2021).
  * Ouladdiaf  _et al._ [2006] B. Ouladdiaf, J. Archer, G. McIntyre, A. Hewat, D. Brau, and S. York, OrientExpress: A new system for Laue neutron diffraction, Physica B: Condensed Matter 385-386, 1052 (2006).
  * Ouladdiaf  _et al._ [2011] B. Ouladdiaf, J. Archer, J. R. Allibon, P. Decarpentrie, M.-H. Lemée-Cailleau, J. Rodríguez-Carvajal, A. W. Hewat, S. York, D. Brau, and G. J. McIntyre, CYCLOPS — a reciprocal-space explorer based on CCD neutron detectors, J. Appl. Cryst. 44, 392 (2011).
  * Qureshi [2019] N. Qureshi, Mag2Pol: a program for the analysis of spherical neutron polarimetry, flipping ratio and integrated intensity data, J. Appl. Cryst. 52, 175 (2019).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Mat. 21, 395502 (2009).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * [22] See Supplemental Material at http://link.aps.org/supplemental/… magnetometry measurements, and the structural refinement.
  * Soh  _et al._ [2019b] J.-R. Soh, P. Manuel, N. M. B. Schröter, C. J. Yi, F. Orlandi, Y. G. Shi, D. Prabhakaran, and A. T. Boothroyd, Magnetic and electronic structure of Dirac semimetal candidate EuMnSb2, Phys. Rev. B 100, 174406 (2019b).
  * Gong  _et al._ [2020] D. Gong, S. Huang, F. Ye, X. Gui, J. Zhang, W. Xie, and R. Jin, Canted Eu magnetic structure in EuMnSb2, Phys. Rev. B 101, 224422 (2020).
  * Liu [2017] J. Liu, Quantum transport and magnetic properties of topological semimetals in AMnSb2 (A = Sr, Ba, and Yb), PhD Thesis, Tulane University (2017).
  * Liu  _et al._ [2017a] J. Y. Liu, J. Hu, Q. Zhang, D. Graf, H. B. Cao, S. M. A. Radmanesh, D. J. Adams, Y. L. Zhu, G. F. Cheng, X. Liu, W. A. Phelan, J. Wei, M. Jaime, F. Balakirev, D. A. Tennant, J. F. DiTusa, I. Chiorescu, L. Spinu, and Z. Q. Mao, A magnetic topological semimetal Sr1-yMn1-zSb2 ($y$, $z<0.1$), Nat. Mater. 16, 905 (2017a).
  * He  _et al._ [2017] J. B. He, Y. Fu, L. X. Zhao, H. Liang, D. Chen, Y. M. Leng, X. M. Wang, J. Li, S. Zhang, M. Q. Xue, C. H. Li, P. Zhang, Z. A. Ren, and G. F. Chen, Quasi-two-dimensional massless Dirac fermions in CaMnSb2, Phys. Rev. B 95, 045128 (2017).
  * [28] Q. Zhang, J. Liu, H. Cao, W. A. Phelan, J. F. DiTusa, D. A. Tennant, and Z. Mao, Tuning quantum transport by controlling spin reorientations in Dirac semimetal candidates Eu1-xSrxMnSb2, arXiv2010.10405 .
  * Liu  _et al._ [2016a] J. Liu, J. Hu, H. Cao, Y. Zhu, A. Chuang, D. Graf, D. J. Adams, S. M. A. Radmanesh, L. Spinu, I. Chiorescu, and Z. Mao, Nearly massless Dirac fermions hosted by Sb square net in BaMnSb2, Sci. Rep. 6, 10.1038/srep30525 (2016a).
  * Guo  _et al._ [2014] Y. F. Guo, A. J. Princep, X. Zhang, P. Manuel, D. Khalyavin, I. I. Mazin, Y. G. Shi, and A. T. Boothroyd, Coupling of magnetic order to planar Bi electrons in the anisotropic Dirac metals $A$MnBi2 ($A$ = Sr, Ca), Phys. Rev. B 90, 075120 (2014).
  * Wang  _et al._ [2016] A. Wang, I. Zaliznyak, W. Ren, L. Wu, D. Graf, V. O. Garlea, J. B. Warren, E. Bozin, Y. Zhu, and C. Petrovic, Magnetotransport study of Dirac fermions in YbMnBi2 antiferromagnet, Phys. Rev. B 94, 165161 (2016).
  * Liu  _et al._ [2017b] J. Y. Liu, J. Hu, D. Graf, T. Zou, M. Zhu, Y. Shi, S. Che, S. M. A. Radmanesh, C. N. Lau, L. Spinu, H. B. Cao, X. Ke, and Z. Q. Mao, Unusual interlayer quantum transport behavior caused by the zeroth Landau level in YbMnBi2, Nat. Commun. 8, 646 (2017b).
  * Soh  _et al._ [2019c] J.-R. Soh, H. Jacobsen, B. Ouladdiaf, A. Ivanov, A. Piovano, T. Tejsner, Z. Feng, H. Wang, H. Su, Y. Guo, Y. Shi, and A. T. Boothroyd, Magnetic structure and excitations of the topological semimetal YbMnBi2, Phys. Rev. B 100, 144431(R) (2019c).
  * Boothroyd [2020] A. T. Boothroyd,  _Principles of Neutron Scattering from Condensed Matter_ (Oxford University Press, 2020).
  * Park  _et al._ [2011] J. Park, G. Lee, F. Wolff-Fabris, Y. Y. Koh, M. J. Eom, Y. K. Kim, M. A. Farhan, Y. J. Jo, C. Kim, J. H. Shim, and J. S. Kim, Anisotropic Dirac fermions in a Bi square net of SrMnBi2, Phys. Rev. Lett. 107, 126402 (2011).
  * Ramankutty  _et al._ [2018] S. V. Ramankutty, J. Henke, A. Schiphorst, R. Nutakki, S. Bron, G. Araizi-Kanoutas, S. K. Mishra, L. Li, Y. K. Huang, T. K. Kim, M. Hoesch, C. Schlueter, T. L. Lee, A. de Visser, Z. Zhong, J. van Wezel, E. van Heumen, and M. S. Golden, Electronic structure of the candidate 2D Dirac semimetal SrMnSb2: a combined experimental and theoretical study, SciPost Phys. 4, 010 (2018).
  * Huang  _et al._ [2017] S. Huang, J. Kim, W. A. Shelton, E. W. Plummer, and R. Jin, Nontrivial berry phase in magnetic bamnsb2 semimetal, Proceedings of the National Academy of Sciences 114, 6256 (2017), https://www.pnas.org/content/114/24/6256.full.pdf .
  * Liu  _et al._ [2016b] J. Liu, J. Hu, H. Cao, Y. Zhu, A. Chuang, D. Graf, D. J. Adams, S. M. A. Radmanesh, L. Spinu, I. Chiorescu, and Z. Mao, Nearly massless dirac fermions hosted by sb square net in bamnsb2, Scientific Reports 6, 30525 (2016b).
  * [39] The neutron scattering data recorded on the D3 diffractometer at the Institut Laue–Langevin is available from https://doi.ill.fr/10.5291/ILL-DATA.5-53-305.


# Synthetic Aperture Radar Image Change Detection via Layer Attention-Based
Noise-Tolerant Network

Desen Meng, Feng Gao, Junyu Dong, Qian Du, and Heng-Chao Li This work was
supported in part by the National Key Research and Development Program of
China under Grant 2018AAA0100602, and in part by the Natural Science
Foundation of Shandong Province under Grant ZR2019QD011. Desen Meng, Feng Gao,
and Junyu Dong are with the School of Information Science and Engineering,
Ocean University of China, Qingdao 266100, China. _(Corresponding author: Feng
Gao.)_ Qian Du is with the Department of Electrical and Computer Engineering,
Mississippi State University, Starkville, MS 39762 USA. H. -C. Li is with the
Sichuan Provincial Key Laboratory of Information Coding and Transmission,
Southwest Jiaotong University, Chengdu 610031, China.

###### Abstract

Recently, change detection methods for synthetic aperture radar (SAR) images
based on convolutional neural networks (CNN) have gained increasing research
attention. However, existing CNN-based methods neglect the interactions among
multilayer convolutions, and errors involved in the preclassification restrict
the network optimization. To this end, we proposed a layer attention-based
noise-tolerant network, termed LANTNet. In particular, we design a layer
attention module that adaptively weights the feature of different convolution
layers. In addition, we design a noise-tolerant loss function that effectively
suppresses the impact of noisy labels. Therefore, the model is insensitive to
noisy labels in the preclassification results. The experimental results on
three SAR datasets show that the proposed LANTNet performs better compared to
several state-of-the-art methods. The source codes are available at
https://github.com/summitgao/LANTNet.

###### Index Terms:

Change detection; Layer attention; Synthetic aperture radar; Noise-tolerant
loss.

##  I Introduction

Remote sensing image change detection aims to identify regions that changed in
the same area using multitemporal images. In fact, it facilitates various
applications, such as disaster assessment, deforestation monitoring, and
object detection [1] [2] [3]. The synthetic aperture radar (SAR) is an ideal
source for change detection, since the SAR sensor is independent of sunlight
illumination and weather conditions. Therefore, SAR change detection has
received extensive research attention recently.

Figure 1: Comparisons of the traditional method with the proposed LANTNet. (a)
The traditional CNN-based feature extraction model for change detection. They
overlook the interactions among multi-layer convolutions. (b) The proposed
LANTNet. The layer attention module adaptively computes weights for features
from different layers and thus improves the feature representations.

Many supervised and unsupervised methods have been devoted to detect changes
in multitemporal SAR images. In fact, it commonly costs a lot of effort to
collect high-quality labeled samples that reflect the changed information, and
learning the changed information in an unsupervised manner is more efficient.
Therefore, in this letter, we focus on unsupervised SAR change detection.

Traditionally, the task of SAR image change detection are solved by clustering
methods, such as $k$-means [4], fuzzy $c$-means (FCM) [5] and hierarchical
clustering [6]. Recently, leveraging the representation power of convolutional
neural networks (CNNs) and the attention mechanism to capture changed
information has become a common solution for change detection. Gao et al. [7]
designed an adaptive fusion convolution module to encode the input, allowing
the network to focus on the more informative channels. Zhang et al [8]
proposed a sparse feature clustering network for SAR image change detection.
Geng et al. [9] presented a saliency-guided network to suppress the influence
of speckle noise and improve the SAR change detection performance. In [10],
transfer learning is employed for the parameter initialization for change
detection deep networks. Furthermore, wavelet-CNN [11], discrete cosine
transform [12], and ensemble learning [13] are also used to enhance the
representation power of deep networks.

Figure 2: The pipeline of LANTNet. It is comprised of the preclassification,
convolution stem, layer attention module and the FC layer. The convolution
stem consists of four convolutional layers. Correspondingly, the layer
attention module uses four feature groups to exploit the correlations among
the four-layer features generated by the convolution stem.

Indeed, the CNN-based models have shown impressive performance in change
detection task. However, they still show limitations in the following aspects:
1) Feature interaction among convolution layers. Research in human vision
system finds that visual representations in the cortex are activated in
parallel fashion [14]. Existing self-attention-enhanced CNNs commonly
incorporate an attention block after a convolution layer, which overlooks the
interactions among multi-layer convolutions. In light of this, it is essential
to explore the information flow among the convolutional layers. 2) Errors
involved in the pseudo-labels. Existing unsupervised change detection methods
generally use clustering to generate pseudo-labels with high confidence for
training. However, the errors involved in the pseudo-labels restrict network
optimization [15]. Therefore, it is vital to build a noise-tolerant model for
unsupervised change detection task.

To address the aforementioned issues, we propose a Layer Attention-based
Noise-Tolerant Network, dubbed as LANTNet, which interconnects convolutional
layers and is less sensitive to noisy labels. Specifically, we design a layer
attention module for change detection, which adaptively weights features from
different layers, as illustrated in Fig. 1. In addition, we design a noise-
tolerant loss function that can effectively suppress the impact of noisy
labels. It combines the Cross Entropy (CE) loss and the Mean Absolute Error
(MAE). Hence, the network is less sensitive to noisy labels and has a fast
convergence rate. We conducted extensive experiments on three SAR datasets,
and the results demonstrate that the proposed LANTNet yields better
performance as compared to state-of-the-art methods.

In summary, the contributions of this letter are threefold:

  * •

We propose a layer attention module, which exploits the correlations of multi-
layer convolutions. Features from different convolution layers cooperate with
each other, effectively improving the representation power of the network.

  * •

We introduce a noise-tolerant loss function to alleviate the influence of
noisy labels in pseudo-labeled samples. It makes the model less sensitive to
noisy labels and has fast convergence rate.

  * •

Experimental results confirm the effectiveness of the proposed LANTNet.
Furthermore, we have made our data and codes publicly available to benefit
other research of the remote sensing community.

##  II Methodology

Given the multitemporal SAR images $I_{1}$ and $I_{2}$ acquired in the same
geographical zone, the goal of change detection is to generate a change map
where each spatial location is assigned to one type of change. In this work,
we consider the generation of change map as a binary classification task,
which means that the label is either 1 (changed) or 0 (not changed).

As shown in Fig. 2, the proposed LANTNet comprises three components:
preclassification module, convolutional stem and layer attention module.
First, multitemporal coregistered SAR images $I_{1}$ and $I_{2}$ are fed into
the preclassification module to obtain samples that are highly likely to be
changed or unchanged. New image patches will be generated and then fed into
the convolutional stem to compute feature groups $X$. The feature groups are
then updated to attention features $Y$ by the layer attention module. Finally,
$Y$ and the features of the convolutional stem are combined for classification
to calculate the final change map.

###  II-A Preclassification and Convolutional Stem

In the preclassification module, the logarithmic ratio operator and the
hierarchical Fuzzy C-Means Algorithm are employed, and details of the module
can be found in [6]. After preclassification, the image patches ($R\times R$
pixels) with pseudo-labels are extracted from $I_{1}$, $I_{2}$, and the
difference image, respectively. They are concatenated to form new image
patches ($3\times R\times R$ pixels) as training samples.

The convolutional stem consists of four convolution layers. 1$\times$1
convolution is first used to extract the shallow feature $F_{0}$. Then, three
3$\times$3 convolution layers are employed to extract intermediate features
$F_{i}$ as

$F_{i}=\textrm{Conv}(F_{i-1}),~{}~{}~{}i=1,2,3.$ (1)

For the shallow feature $F_{0}$, we set the channel numbers to 16, and for
intermediate features $F_{1}$, $F_{2}$, and $F_{3}$, we set the channel
numbers to 32. Then, $1\times 1$ convolution is employed to increase the
channel dimensions of $F_{0}$ to 32. After that, $F_{0}$, $F_{1}$, $F_{2}$,
and $F_{3}$ are combined to form a feature group, and fed into the layer
attention module.

###  II-B Layer Attention Module

The attention mechanism helps the network focus on the most important part and
neglect irrelevant information. However, the correlations among multi-layer
convolutions are rarely taken into account. To this end, we design a layer
attention module to weight the features $F_{0}$, $F_{1}$, $F_{2}$, and $F_{3}$
from multiple layers by learning the correlation matrix. Therefore, layer
attention effectively fuses the spatial information of low-level layers and
the semantic information from high-level layers. The informative layers are
emphasized, and the redundant ones are suppressed. Fig. 2 shows the details of
the proposed layer attention module. The input of the module is the feature
group $X\in\mathbb{R}^{N\times C\times R\times R}$ from the convolution stem.
In our implementation, $N$ is set to 4. Then we reshape $X$ to a feature
matrix with the dimension of $N\times CRR$. Subsequently, we apply matrix
multiplication between the new feature matrix and a diagonal matrix $W$ to
weight the features from different layers. Note that $W$ is initialized to
identity matrix and only the diagonal elements of $W$ are updated in each
epoch. The weight matrix adaptively assign weights to the input feature
groups. After that, a weighted feature matrix $\hat{X}\in\mathbb{R}^{N\times
CRR}$ can be obtained.

Next, we calculate an attention matrix $A\in\mathbb{R}^{N\times N}$ as the
correlation matrix as follows:

$A=\textrm{Softmax}(\hat{X}\otimes\hat{X}^{T}),$ (2)

where the element $A_{i,j}$ in the attention matrix denotes the correlation
index between the $i$-th and $j$-th layer features, and $\otimes$ denotes the
matrix multiplication. Matrix multiplication is performed between the weighted
feature matrix $\hat{X}$ and its transpose matrix $\hat{X}^{T}$. Softmax
operation is applied to generate the attention matrix $A$.

Finally, we multiply the weighted feature matrix $\hat{X}$ by the attention
matrix $A$, and reshape the output to the dimension of $N\times C\times
R\times R$. The input $X$ are combined to form the final output $Y$ as:

$Y=\phi(A\otimes\hat{X})\oplus X,$ (3)

where $\phi$ denotes the reshape operation.

To compute the final change map, $Y$ is reshaped to the dimension of $NC\times
R\times R$, and then fed into a $1\times 1$ convolution layer for channel
reduction, as illustrated in Fig. 2. After that, the combined features are
handled by a fully connected (FC) layer. The softmax function is used to
calculate the possibility of changed or unchanged. After training, the final
change map can be obtained by classifying all pixels from the multitemporal
SAR images.

###  II-C Noise-tolerant loss function

Like most existing unsupervised methods, we employ clustering to generate
pseudo-labels with high probability for network training. However, errors are
commonly involved in these pesudo-labels, and restrict the network
optimization [15]. Therefore, it is critical to build a noise-tolerant model.

Toward this end, we design a noise-tolerant loss which combines the Cross
Entropy (CE) loss and the Mean Absolute Error (MAE). Suppose the dataset for
change detection is represensted as $\\{(x_{i},y_{i})|1\leq i\leq m\\}$, where
$x_{i}$ is the input image patch, $y_{i}$ is the label associated with
$x_{i}$, and $m$ is the total number of samples. The MAE loss is defined as:

$\mathcal{L}_{\textrm{MAE}}(f(x),y)=\|e_{y}-f(x)\|_{1},$ (4)

where $f(\cdot)$ represents the network and $e_{y}$ is a one-hot vector with
$e_{yj}=1$ if $j=y$, otherwise $e_{yj}=0$.

The CE loss is computed as:

$\mathcal{L}_{\textrm{CE}}(f(x),y)=-e_{y}\log f(x).$ (5)

Ghosh et al. [16] empirically demonstrated that MAE is robust to noisy labels,
but it often costs more time to converge. To the contrary, CE tends to overfit
noisy labels, but it performs better than MAE when training with clean data.
To leverage the good convergence performance of CE, and exploit the noise
robustness of provided by MAE, we design a robust loss function which combines
the CE and MAE. Formally, the robust loss function is defined as:

$\mathcal{L}=\alpha\mathcal{L}_{\textrm{CE}}+\beta\mathcal{L}_{\textrm{MAE}}$
(6)

where $\alpha$ and $\beta$ are two decoupled hyperparameters. In our
implementation, we set $\alpha$ to 0.1 and $\beta$ to 0.9 for a trade-off
between convergence efficiency and the noise robustness.

##  III Experiments and Analysis

###  III-A Datasets and Evaluation Metrics

To evaluate the proposed LANTNet, we conducted extensive experiments on three
real SAR datasets captured by different sensors. The first dataset is Chao
Lake dataset with the size of $384\times 384$ pixels. It was captured by the
Sentinel-1 sensor in May and July 2020 in China, respectively. The second
dataset is the Ottawa dataset with the size of $290\times 350$ pixels, which
was acquired by the Radarsat sensor in May and August 1997, respectively. The
Ottawa dataset shows the changed areas afflicted by floods. The third dataset
is the Sulzberger dataset which was captured by the Envisat satellite of the
European Space Agency on March 11 and 16, 2011. Both images show a breakup of
the ice shelf caused by a tsunami.

Five commonly-used evaluation metrics, i.e., false positive (FP), false
negative (FN), overall error (OE), percentage of correct classification (PCC)
and Kappa coefficient (KC) are used to evaluate the change detection
performance.

###  III-B Analysis of the Patch Size

The image patch size $R$ is a critical parameter. To explore the relationship
between the patch size $R$ and PCC, we set $R$ from 5 to 15, and the
experimental results are shown in Fig. 3. It can be seen that all curves tend
to increase first and then decrease. When $R$ is small, the contextual
information around each pixel is insufficient. However, if $R$ becomes too
large, irrelevant neighborhood information would be taken into account.
Therefore, we set $R=9$ for the Chao Lake dataset, and $R=7$ for the Ottawa
and Sulzberger datasets.

Figure 3: Relationship between image patch size and the PCC value. Figure 4: Visualization of the feature representations on the Ottawa dataset. (a) Features before the layer attention module. (b) Features after the layer attention module. TABLE I: Ablation Studies of the Proposed LANTNet Method | PCC on different datasets ($\%$)
---|---
Chao Lake | Ottawa | Sulzberger
Basic CNN | 97.49 | 98.22 | 96.36
w/o layer attention | 97.82 | 98.32 | 96.64
w/o noise-tolerant loss | 97.98 | 98.42 | 96.69
Proposed LANTNet | 98.49 | 98.47 | 96.93
Figure 5: Visualized results of different change detection methods on the Chao
Lake dataset (first row), Ottawa dataset (second row) and Sulzberger dataset
(third row) : (a) Image captured at $t_{1}$. (b) Image captured at $t_{2}$.
(c) Ground truth image. (d) Result by PCAKM [4]. (e) Result by NR-ELM [17].
(f) Result by GaborPCANet [6]. (g) Result by CWNN [11]. (h) Result by DDNet
[12]. (i) Result by the proposed LANTNet.

###  III-C Ablation Study

To verify the effectiveness of the proposed layer attention module and the
noise-tolerant loss, we conduct ablation experiments on three datasets. Three
variants are designed for comparison. First, we design a basic CNN with the
same structure as LANTNet, which removes the layer attention module and uses
cross entropy as loss function instead of the noise-tolerant loss. Then we run
our model without layer attention module (_w/o layer attention_) and without
the noise-tolerant loss(_w/o noise-tolerant loss_). The results are shown in
Table I. It can be observed that both the layer attention module and the
noise-tolerant loss are effective for change detection performance
improvement.

Furthermore, to gain deeper insights toward the layer attention, we visualized
the features before and after the layer attention module with the tool of
t-SNE. As shown in Fig. 4, the distance between features clusters derived by
the layer attention is larger, which reflects the feature representations by
our model are indeed more discriminative.

TABLE II: Change Detection Results of Different Methods on Three Datasets Method | Results on the Chao Lake dataset
---|---
FP | FN | OE | PCC | KC
PCAKM [4] | 8521 | 2248 | 10769 | 92.70 | 65.58
NR-ELM [17] | 595 | 3836 | 4431 | 97.00 | 81.27
GaborPCANet [6] | 2946 | 1771 | 4717 | 96.80 | 82.66
CWNN [11] | 959 | 2397 | 3356 | 97.72 | 86.63
DDNet [12] | 2097 | 1068 | 3165 | 97.85 | 88.31
Proposed LANTNet | 1358 | 867 | 2225 | 98.49 | 91.65
Method | Results on the Ottawa dataset
FP | FN | OE | PCC | KC
PCAKM [4] | 729 | 2458 | 3187 | 96.86 | 87.67
NR-ELM [17] | 517 | 1310 | 1827 | 98.20 | 93.10
GaborPCANet [6] | 755 | 1080 | 1835 | 98.19 | 93.15
CWNN [11] | 1291 | 434 | 1725 | 98.30 | 93.75
DDNet [12] | 504 | 1191 | 1695 | 98.33 | 93.62
Proposed LANTNet | 638 | 911 | 1549 | 98.47 | 94.23
Method | Results on the Sulzberger dataset
FP | FN | OE | PCC | KC
PCAKM [4] | 3308 | 701 | 4009 | 93.88 | 84.49
NR-ELM [17] | 2386 | 646 | 3032 | 95.37 | 88.07
GaborPCANet [6] | 2485 | 494 | 2979 | 95.45 | 88.34
CWNN [11] | 1598 | 1132 | 2730 | 95.83 | 88.98
DDNet [12] | 1710 | 526 | 2236 | 96.59 | 91.10
Proposed LANTNet | 1229 | 784 | 2013 | 96.93 | 91.87

###  III-D Experimental Results and Comparison

To verify the effectiveness of the proposed LANTNet, five closely related
change detection methods are used for comparison, including PCAKM [4], NR-ELM
[17], GaborPCANet [6], CWNN [11], and DDNet [12]. Note that the methods
mentioned above are implemented with the default parameters provided in their
work. All visual results are illustrated in Fig. 5, and the corresponding
quantitative results are shown in Table II.

On the Chao Lake dataset (the first row of Fig. 5), PCAKM has the highest FP
value among all methods and the FP value affects the overall performance of
PCAKM. At the same time, GaborPCANet and DDNet also suffer from high FP
values, while NR-ELM and CWNN suffer from high FN values. CWNN and DDNet are
two CNN-based methods and have strong feature representation power, as shown
in Table II. However, the pseudo-labels utilized by CWNN and DDNet are not
accurate enough, and their models may overfit the noisy labels. The KC value
of the proposed LANTNet is impoved by 26.07$\%$, 10.38$\%$, 8.99$\%$, 5.02$\%$
and 3.34$\%$ over PCAKM, NR-ELM, GaborPCANet, CWNN and DDNet, respectively. It
is evident that the proposed method achieves the best performance on the Chao
Lake dataset.

On the Ottawa dataset (the second row of Fig. 5), PCAKM, NR-ELM, GaborPCANet,
and DDNet suffer from high FN values and many change pixels are not correctly
classified. CWNN suffers from a high FP value and generates many noisy spots
in the change map. The proposed LANTNet achieves the best PCC value, and it
demonstrates that layer attention effectively improves the feature
representation. The comparisons demonstrate the good performance of the
proposed LANTNet on the Ottawa dataset.

On the Sulzberger dataset (the third row of Fig. 5), all the methods suffer
from high FP values, and the proposed LANTNet achieves the lowest FP value. It
can be seen from Fig. 5 that the proposed LANTNet is effective in suppressing
the label noise and reduces the FP values. The KC value of the proposed
LANTNet is impoved by 7.38$\%$, 3.8$\%$, 3.53$\%$, 2.89$\%$ and 0.77$\%$ over
PCAKM, NR-ELM, GaborPCANet, CWNN and DDNet, respectively. The comparison shows
that the LANTNet is powerful in noise suppression on the Sulzberger dataset.

##  IV Conclusions

In this letter, we propose a layer attention-based noise-tolerent neural
network (LANTNet) for SAR change detection. Compared with existing CNN-based
methods, our LANTNet improves in two aspects. First, we design a layer
attention module that adaptively weights the feature of different convolution
layers. In addition, since the pseudo-labels obtained by preclassification are
not accurate enough, we design a noise-tolerant loss to help the network avoid
overfitting to the noisy labels. Experimental results on three multi-temporal
SAR datasets demonstrate that the proposed LANTNet is superior to several
state-of-the-art change detection methods.

## References

  * [1] Z. Lv, X. Yang, X. Zhang, and J. A. Benediktsson, “Object-based sorted-histogram similarity measurement for detecting land cover change with VHR remote sensing images,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 19, pp. 1–5, 2022.
  * [2] J. Wang, F. Gao, J. Dong, S. Zhang, and Q. Du, “Change detection from synthetic aperture radar images via graph-based knowledge supplement network,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 15, pp. 1823–1836, 2022.
  * [3] J. Ai, R. Tian, Q. Luo, J. Jin, and B. Tang, “Multi-scale rotation-invariant Haar-like feature integrated CNN-based ship detection algorithm of multiple-target environment in SAR imagery,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 57, no. 12, pp. 10 070–10 087, 2019.
  * [4] T. Celik, “Unsupervised change detection in satellite images using principal component analysis and $k$-means clustering,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 6, no. 4, pp. 772–776, 2009.
  * [5] M. Gong, L. Su, M. Jia, and W. Chen, “Fuzzy clustering with a modified MRF energy function for change detection in synthetic aperture radar images,” _IEEE Transactions on Fuzzy Systems_ , vol. 22, no. 1, pp. 98–109, 2014.
  * [6] F. Gao, J. Dong, B. Li, and Q. Xu, “Automatic change detection in synthetic aperture radar images based on PCANet,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 13, no. 12, pp. 1792–1796, 2016.
  * [7] Y. Gao, F. Gao, J. Dong, and H.-C. Li, “SAR image change detection based on multiscale capsule network,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 18, no. 3, pp. 484–488, 2021.
  * [8] W. Zhang, L. Jiao, F. Liu, S. Yang, W. Song, and J. Liu, “Sparse feature clustering network for unsupervised SAR image change detection,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–13, 2022.
  * [9] J. Geng, X. Ma, X. Zhou, and H. Wang, “Saliency-guided deep neural networks for SAR image change detection,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 57, no. 10, pp. 7365–7377, 2019.
  * [10] Y. Gao, F. Gao, J. Dong, and S. Wang, “Transferred deep learning for sea ice change detection from synthetic aperture radar images,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 16, no. 10, pp. 1655–1659, 2019\.
  * [11] F. Gao, X. Wang, Y. Gao, J. Dong, and S. Wang, “Sea ice change detection in SAR images based on convolutional-wavelet neural networks,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 16, no. 8, pp. 1240–1244, 2019.
  * [12] X. Qu, F. Gao, J. Dong, Q. Du, and H.-C. Li, “Change detection in synthetic aperture radar images using a dual-domain network,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 19, pp. 1–5, 2022.
  * [13] J. Wang, Y. Wang, B. Chen, and H. Liu, “LCS-EnsemNet: A semisupervised deep neural network for SAR image change detection with dual feature extraction and label-consistent self-ensemble,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 11 903–11 925, 2021.
  * [14] L. Chelazzi, J. Duncan, E. Miller, and R. Desinone, “Responses of neurons in inferior temporal cortex during memory-guided visual search,” _Journal of Neurophysiology_ , vol. 80, no. 6, pp. 2918–2940, 1998.
  * [15] Z. Zhang and M. Sabuncu, “Generalized cross entropy loss for training deep neural networks with noisy labels,” in _32nd Conference on Neural Information Processing Systems_ , 2018, pp. 1–14.
  * [16] A. Ghosh, H. Kumar, and P. Sastry, “Robust loss functions under label noise for deep neural networks,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2017, pp. 1919–1925.
  * [17] F. Gao, J. Dong, B. Li, Q. Xu, and C. Xie, “Change detection from synthetic aperture radar images based on neighborhood-based ratio and extreme learning machine,” _Journal of Applied Remote Sensing_ , vol. 10, no. 4, pp. 1–14, 2016.


SI.pdf

# On the dynamical stability of copper-doped lead apatite

Sun-Woo Kim<EMAIL_ADDRESS>Department of Materials Science and Metallurgy,
University of Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, United
Kingdom  Kang Wang  Department of Materials Science and Metallurgy, University
of Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom  Siyu
Chen  Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue,
Cambridge CB3 0HE, United Kingdom  Lewis J. Conway  Department of Materials
Science and Metallurgy, University of Cambridge, 27 Charles Babbage Road,
Cambridge CB3 0FS, United Kingdom  Advanced Institute for Materials Research,
Tohoku University, 2-1-1 Katahira, Aoba, Sendai 980-8577, Japan  G. Lucian
Pascut  MANSiD Research Center and Faculty of Forestry, Stefan Cel Mare
University (USV), Suceava 720229, Romania  Ion Errea  Fisika Aplikatua Saila,
Gipuzkoako Ingeniaritza Eskola, University of the Basque Country (UPV/EHU),
Europa Plaza 1, 20018 Donostia/San Sebastián, Spain  Centro de Física de
Materiales (CSIC-UPV/EHU), Manuel de Lardizabal Pasealekua 5, 20018
Donostia/San Sebastián, Spain  Donostia International Physics Center (DIPC),
Manuel de Lardizabal Pasealekua 4, 20018 Donostia/San Sebastián, Spain  Chris
J. Pickard  Department of Materials Science and Metallurgy, University of
Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom
Advanced Institute for Materials Research, Tohoku University, 2-1-1 Katahira,
Aoba, Sendai 980-8577, Japan  Bartomeu Monserrat<EMAIL_ADDRESS>Department
of Materials Science and Metallurgy, University of Cambridge, 27 Charles
Babbage Road, Cambridge CB3 0FS, United Kingdom  Cavendish Laboratory,
University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, United
Kingdom

###### Abstract

The recent claim of room temperature superconductivity in a copper-doped lead
apatite compound, called LK-99, has sparked remarkable interest and
controversy. Subsequent experiments have largely failed to reproduce the
claimed superconductivity, while theoretical works have identified multiple
key features including strong electronic correlation, structural
instabilities, and dopability constraints. A puzzling claim of several recent
theoretical studies is that both parent and copper-doped lead apatite
structures are dynamically unstable at the harmonic level, questioning decades
of experimental reports of the parent compound structures and the recently
proposed copper-doped structures. In this work, we demonstrate that both
parent and copper-doped lead apatite structures are dynamically stable at room
temperature. Anharmonic phonon-phonon interactions play a key role in
stabilizing some copper-doped phases, while most phases are largely stable
even at the harmonic level. We also show that dynamical stability depends on
both volume and correlation strength, suggesting controllable ways of
exploring the copper-doped lead apatite structural phase diagram. Our results
fully reconcile the theoretical description of the structures of both parent
and copper-doped lead apatite with experiment.

##  I Introduction

Copper-doped lead apatite Pb10-xCux(PO4)6O with $0.9<x<1.1$, known as LK-99,
has been recently claimed to exhibit superconductivity above room temperature
and at ambient pressure LK99-1 ; LK99-2 . This remarkable claim is backed by
magnetic (half-)levitation on a permanent magnet and by a sudden drop in
resistivity at the claimed superconducting transition temperature. However,
subsequent extensive experimental efforts by other groups have failed to
confirm the superconductivity Kumar_2023_synthesis ; Liu_2023_semiconducting ;
wu2023successful ; Guo2023 ; Wang2023 ; zhang2023ferromagnetism ; zhu2308first
; timokhin2023synthesis ; kumar2023absence ; Liu_2023_phases ;
puphal2023single . The magnetic half-levitation is reproduced in some
insulating samples, where it is attributed to soft ferromagnetism Guo2023 ;
Wang2023 ; zhang2023ferromagnetism . A plausible explanation for the sudden
resistivity drop is provided by a first order phase transition of Cu2S
impurities zhu2308first ; Jain2023 , which is further supported by the highly
insulating nature of a single crystalline sample without Cu2S impurities
puphal2023single .

On the theoretical front, initial density functional theory calculations
reported an electronic structure exhibiting relatively flat bands near the
Fermi level for a simple model of copper-doped lead apatite griffin2023origin
; si2023electronic ; LAI202466 ; kurleto2023pb . However, subsequent
calculations showed that the inclusion of spin-orbit coupling or non-local
correlations lead to an insulating electronic structure bai2023ferromagnetic ;
swift2023comment ; pashov2023multiple , a conclusion that is also reached with
the inclusion of local correlations using dynamical mean-field theory
korotin2023electronic ; si2308pb10 ; yue2023correlated . Different estimates
of critical superconducting temperatures have so far delivered values
significantly lower than room temperature oh2023s ; witt2023no ;
paudyal2023implications .

These state-of-the-art electronic structure calculations all assume a specific
structural model as a starting point, often suggested by experiment. However,
other theoretical works have questioned the suitability of these structural
models both in terms of the thermodynamic feasibility of copper doping
shen2023phase  or the dynamical stability of the experimentally proposed
structures Bernevig_phonon ; shen2023phase ; hao2023first ; liu2023symmetry ;
liu2023different ; cabezas2023theoretical . Indeed, one of the most basic
quantities used to characterize a material is its dynamical stability. A
dynamically stable structure corresponds to a local minimum of the potential
(free) energy surface, and its phonon frequencies are real. A dynamically
unstable structure corresponds to a saddle point of the potential (free)
energy surface, and some of its phonon frequencies are imaginary with
associated eigenvectors that encode atomic displacement patterns that lower
the energy of the system. Only dynamically stable structures can represent
real materials. Puzzlingly, recent computational works have claimed that the
experimentally reported structures of the parent lead apatite Bernevig_phonon
; shen2023phase ; hao2023first ; liu2023symmetry  and of the copper-doped lead
apatite compounds Bernevig_phonon ; shen2023phase ; hao2023first ;
liu2023different ; cabezas2023theoretical  are dynamically unstable at the
harmonic level, which would imply that they cannot be the true structures of
the materials underpinning LK-99, and would question the validity of most
electronic structure calculations to date.

In this work, we demonstrate that both parent lead apatite and copper-doped
lead apatite compounds are dynamically stable at room temperature. The parent
compounds are largely stable at the harmonic level, with some exhibiting very
slight instabilities which are suppressed by quartic anharmonic terms. For the
copper-doped compounds, dynamical stability at the harmonic level depends on
the doping site, but even those that are dynamically unstable at the harmonic
level are overall stable at room temperature with the inclusion of anharmonic
phonon-phonon interactions.

##  II Lead apatite

Lead apatite is a compound that was first experimentally reported over 70
years ago rooksby1952identification . Figure 1 depicts an example of lead
apatite, with a hexagonal lattice (space group $P6_{3}/m$) and general formula
Pb10(PO4)6X2, where X is either a halide atom or an OH group. The variant
Pb10(PO4)6O, which is claimed to be the parent structure of LK-99 LK99-1 ;
LK99-2 , has also been known experimentally for decades
rooksby1952identification ; merker1960oxypromorphit ; merker1970lead ;
krivovichev2003crystal . The X site corresponds to Wyckoff position $4e$,
giving a multiplicity of four in the unit cell, but these sites are only
partially filled with an occupation of $\frac{1}{2}$ for halide atoms and the
OH group, and an occupation of $\frac{1}{4}$ for O. We consider two
representative cases, Pb10(PO4)6O and Pb10(PO4)6(OH)2, where the specific
distribution of species on the X site results in space groups $P3$ and
$P6_{3}$, respectively.

Figure 1:  Crystal structure and phonon band structures of parent lead
apatite. a. Crystal structure of lead apatite Pb10(PO4)6O or Pb10(PO4)6(OH)2.
b. Harmonic and anharmonic (50 K) phonon band structures of Pb10(PO4)6O. c.
Harmonic phonon band structure of Pb10(PO4)6(OH)2.

The phonon dispersion of the parent Pb10(PO4)6O and Pb10(PO4)6(OH)2 compounds
is shown in Fig. 1. At the harmonic level, Pb10(PO4)6O exhibits imaginary
phonon frequencies at the zone boundary points M and K of the $k_{z}=0$ plane,
and at the zone boundary point H of the $k_{z}=\frac{\pi}{c}$ plane. However,
the absolute values of these imaginary frequencies are less than $2.3$ meV and
the resulting anharmonic potentials have a dominant quartic term that strongly
suppresses the instability to about $0.2$ meV/f.u. (see Supplementary Fig.
S4). As a result, the calculation of self-consistent phonons including
anharmonic interactions fully stabilizes the structure at the relatively low
temperature of $50$ K, and potentially lower. Earlier works reported that
Pb10(PO4)6O is dynamically unstable at the harmonic level Bernevig_phonon ;
shen2023phase ; hao2023first , a result we confirm, but our work further
demonstrates that Pb10(PO4)6O is overall dynamically stable at room
temperature driven by higher-order anharmonic terms.

Pb10(PO4)6(OH)2 is dynamically stable at the harmonic level. Earlier works
reported that Pb10(PO4)6(OH)2 is dynamically unstable Bernevig_phonon ;
shen2023phase , an opposite conclusion that we attribute to unconverged
harmonic calculations (see Supplementary Fig. S3). We note that to fully
converge the harmonic calculations of Pb10(PO4)6(OH)2, the required coarse
$\mathbf{q}$-point grid includes all of $\Gamma$, M, K, and A points, which
can only be accomplished with a regular grid of minimum size $6\times 6\times
2$ or alternatively a non-uniform Farey grid farey_grid  of minimum size
$(2\times 2\times 2)\cup(3\times 3\times 1)$. In our calculations, we use the
latter as it is computationally more efficient.

Overall, we find that both lead apatite compounds Pb10(PO4)6O and
Pb10(PO4)6(OH)2 are dynamically stable, a conclusion that is in full agreement
with multiple experimental reports of the structure of lead apatite over the
past 70 years rooksby1952identification ; merker1960oxypromorphit ;
merker1970lead ; krivovichev2003crystal ; BRUCKNER1995209 .

##  III Copper doped lead apatite

The claim of room temperature superconductivity in LK-99 is based on copper
doping of lead apatite, with copper replacing about 1 in 10 lead atoms leading
to a Pb9Cu(PO4)6O stoichiometry. There are two symmetrically distinct lead
sites, labelled Pb(1) and Pb(2) in the literature (see Fig. 1a), and doping at
these sites results in structures with the space groups $P3$ and $P1$,
respectively. The original LK-99 work suggested that the doping site is Pb(1)
LK99-1 ; LK99-2 , but subsequent experimental works have suggested that both
Pb(1) and Pb(2) sites can be doped puphal2023single ; Bernevig_phonon .
Computational works find that the relative energy between the two doping sites
depends on the exchange-correlation functional and the magnitude of the
Hubbard $U$ parameter used on the copper atom, with most choices favouring
doping at the Pb(2) site, a prediction we confirm with our own calculations.
For completeness, in this work we explore doping at both sites.

Figure 2:  Electron and phonon band structures of Pb9Cu(PO4)6O. a,b Electronic
band structures of Pb9Cu(PO4)6O for copper doping at the a Pb(1) and b Pb(2)
sites. NM, FM, SOC represent nonmagnetic, ferromagnetic, and spin-orbit
coupling, respectively. For doping at the Pb(2) site, the initial nonmagnetic
configuration converged to a ferromagnetic configuration in the presence of
spin-orbit coupling. c Harmonic and anharmonic (300 K) phonon band structures
of Pb9Cu(PO4)6O structure with $P3$ symmetry for doping at the Pb(1) site. d
Harmonic phonon band structure of Pb9Cu(PO4)6O structure with $P1$ symmetry
for doping at the Pb(2) site. The data are obtained with $U=3$ eV on the
copper $3d$ orbital.

The electronic structures of Pb9Cu(PO4)6O with copper on the Pb(1) and Pb(2)
sites are shown in Fig. 2. For doping at the Pb(1) site, a non-magnetic
calculation leads to a metallic state in which the Fermi energy crosses four
relatively flat bands (a pair of doubly-degenerate bands). Inclusion of spin-
orbit coupling while maintaining the non-magnetic configuration leads to a
splitting of the pair of doubly-degenerate bands, and the Fermi energy crosses
a pair of singly-degenerate relatively flat bands. A calculation including
spin-orbit coupling and allowing a non-zero magnetic moment leads to a
ferromagnetic configuration in which the system is gapped. The latter
ferromagnetic configuration is the most energetically favourable, but may not
be directly relevant for room temperature experiments as single crystal
measurements suggest the material is a non-magnetic insulator exhibiting a
diamagnetic response with potentially a small ferromagnetic component
puphal2023single . Additionally, the ferromagnetic ordering may be an artifact
of the DFT calculations, as dynamical mean-field theory calculations
korotin2023electronic ; si2308pb10 ; yue2023correlated  suggest a gap opens
due to a Mott-like band splitting without the need of ferromagnetic ordering.
For doping at the Pb(2) site, we also find a metallic state with a single band
crossing the Fermi level in non-magnetic calculations, and a gapped state in
ferromagnetic calculations. For both doping sites, we find that the phonon
dispersion is only weakly affected by the level of electronic structure theory
used (see Supplmenetary Note 3.1), so the discussion below should be largely
independent of the precise electronic structure of the system.

The phonon dispersions of Pb9Cu(PO4)6O with copper on the Pb(1) and Pb(2)
sites are shown in Fig. 2. Doping at the Pb(2) site leads to a dynamically
stable structure at the harmonic level of theory. By contrast, doping at the
Pb(1) site leads to a dynamically unstable structure at the harmonic level
that exhibits two imaginary phonon branches of frequencies about $15i$ meV
across the entire Brillouin zone. This harmonic instability is present
irrespective of the level of theory used, including a Hubbard $U$ parameter on
the copper $d$ orbitals, spin-orbit coupling, and ferromagnetic ordering (see
Supplementary Fig. S5). Importantly, anharmonic phonon-phonon interactions
strongly suppress the instability and the structure becomes dynamically stable
at $300$ K. We reach similar conclusions for copper doping of Pb10(PO4)6(OH)2
(see Supplementary Fig. S6). Overall, copper-doped lead apatite is dynamically
stable at room temperature for doping at either site.

Figure 3:  $\Gamma$-distorted $P1$ structure and its electron and phonon band
structures. a. Potential energy surface along the imaginary phonon modes at
high symmetry points in the Brillouin zone of the $P3$ structure for doping at
the Pb(1) site (see its harmonic phonon dispersion in Fig. 2c). b,c. b
Harmonic phonon and c electron band structures of the $P1$ structure distorted
along the $\Gamma$ mode of the $P3$ structure. d. Nonmagnetic density of
states (DOS) and Cu partial DOS of both the $P3$ and $\Gamma$-distorted $P1$
structures. See also partial DOS of other atoms in Supplementary Fig. S7. The
data are obtained with $U=3$ eV on the copper $3d$ orbital.

The original paper claiming superconductivity in LK-99 suggested that copper
doping of lead apatite occurs on the Pb(1) site LK99-1 ; LK99-2 . As the
associated structure exhibits a dynamical instability at the harmonic level,
we further explore its properties by considering the potential energy surface
along the imaginary phonon modes at high symmetry points in the Brillouin zone
(Fig. 3a). The dominant instability is driven by a $\Gamma$ point phonon mode,
and fully relaxing the structure along this instability leads to a new
structure of $P1$ symmetry, which is dynamically stable at the harmonic level
of theory (Fig. 3b). We ascribe the harmonic stability of the $P1$ structure
to a downward shift of the occupied part of the density of states compared to
the $P3$ structure, dominated by copper-derived orbitals (Fig. 3d). The four
bands (a pair of doubly-degenerate bands) that cross the Fermi level in the
$P3$ structure (Fig. 2a) split under the distortion, such that the resultant
$P1$ structure has a metallic state with a single doubly-degenerate band
crossing the Fermi level in the nonmagnetic configuration (Fig. 3c). The
distorted $P1$ structure becomes an insulator in the presence of ferromagnetic
ordering, similar to lead apatite with copper doping at the Pb(2) site (Fig.
2b).

Figure 4:  Volume-Hubbard $U$ phase diagram. a. Representative harmonic phonon
band structures of the $P3$ structure doped at the Pb(1) site for different
$U$ values and volumes. b. Volume-Hubbard $U$ phase diagram of Pb9Cu(PO4)6O
for doping at the Pb(1) site. The volume of the $P3$ structure is presented on
the $x$-axis as a reference, corresponding to a pressure range of $0-25$ GPa.
We find slightly larger volume changes in the $\Gamma$-distorted $P1$
structure. The color bar indicates the enthalpy difference between the $P3$
and $P1$ structures (in meV/atom).

Interestingly, the relative energy of the $P3$ structure compared to the
$\Gamma$-distorted $P1$ structure is strongly dependent on both the volume and
the electronic correlation strength as measured by the Hubbard $U$ parameter.
Specifically, we find that harmonic instabilities favouring the $P1$ phase
occur for large values of $U$ and large volumes, while the harmonic
instability of the $P3$ phase completely disappears for small values of $U$
and small volumes. This is evident from the phonon dispersion of the $P3$
phase for different $U$ values and volumes (Fig. 4a) and from the enthalpy
difference between $P3$ and $P1$ phases indicated by the colour bar in the
phase diagram in Fig. 4b. These observations suggest that controlling volume,
for example through hydrostatic pressure or strain, and controlling the degree
of electronic correlation, for example by applying a gate voltage or doping,
can be used to navigate the structural phase diagram of compounds based on
lead apatite. Specifically, it may be possible to observe a temperature-driven
structural phase transition between a low temperature $P1$ phase and high
temperature $P3$ phase in a regime with a large harmonic dynamical
instability. Finally, we note that electronic correlation beyond the static
description provided by a Hubbard $U$ correction may play an important role on
this phase diagram korotin2023electronic ; si2308pb10 ; yue2023correlated , so
further work is required to fully characterise it.

##  IV Discussion

Since the original claim of room temperature superconductivity in LK-99, seven
phonon dispersion calculations have been reported in the literature. Of these,
parent Bernevig_phonon ; shen2023phase ; hao2023first ; liu2023symmetry  and
copper-doped lead apatite with a $P3$ space group Bernevig_phonon ;
shen2023phase ; hao2023first ; liu2023different ; cabezas2023theoretical  are
claimed to be dynamically unstable at the harmonic level, while another work
claims that the copper-doped lead apatite is dynamically stable at the
harmonic level paudyal2023implications .

We attribute these puzzling and contradictory conclusions about the dynamical
stability of lead apatite to the complexity of harmonic phonon calculations in
this system, with a unit cell containing at least $41$ atoms, and to the
subtle interplay between volume, electronic correlation strength, and phonons.
First, we find that fully converged phonon calculations for the parent
compounds require relatively large coarse $\mathbf{q}$-point grids. For
example, a converged calculation for Pb10(PO4)6(OH)2 requires the inclusion of
the $\Gamma$, M, K, and A points in the coarse grid, but none of the
previously reported calculations include all these points, and as a result
they incorrectly conclude that this compound is dynamically unstable at the
harmonic level. Second, we find that dynamical stability for the copper-doped
compounds at the harmonic level depends on both the value of the Hubbard $U$
parameter and the volume of the system (see Fig. 4). In this context, we
rationalise the seemingly contradictory conclusions about dynamical stability
of the copper-doped compounds by suggesting that different works use different
volumes and different choices for the Hubbard $U$ parameter.

Beyond clarifying the dynamical stability of lead apatite at the harmonic
level, we have shown that anharmonic phonon-phonon interactions play a key
role in stabilising multiple lead apatite compounds. Overall, our calculations
indicate that both parent and copper-doped lead apatite compounds are
dynamically stable at room temperature.

We believe that lead apatite is a nice example to illustrate the ability of
state-of-the-art first principles methods to fully characterise a complex and
experimentally relevant system. However, our work also demonstrates that
reliable results and conclusions can only be reached with a careful
consideration of convergence parameters, such as the size of the
$\mathbf{q}$-point grid, and physical models, such as the inclusion of
anharmonic phonon-phonon interactions.

##  V Conclusions

We show that the experimentally suggested structures of lead apatite and
copper-doped lead apatite are dynamically stable at room temperature. Most
structures are dynamically stable at the harmonic level, but some key
structures, including the structure claimed to be responsible for
superconductivity at ambient conditions, only becomes dynamically stable with
the inclusion of anharmonic phonon-phonon interactions. Our results resolve a
puzzling suggestion by multiple earlier computational works that claimed that
the experimentally reported structures of both parent and copper-doped lead
apatite compounds were dynamically unstable, and fully reconcile the current
experimental and theoretical description of the structure of lead apatite.

##  VI Methods

Electronic structure calculations. - We perform density functional theory
(DFT) calculations using the Vienna ab initio simulation package vasp VASP1 ;
VASP2 , which implements the projector-augmented wave method PAW . We employ
PAW pseudopotentials with valence configurations $5d^{10}6s^{2}6p^{2}$ for
lead, $3d^{10}4s^{1}$ for copper, $3s^{2}3p^{3}$ for phosphorus,
$2s^{2}2p^{4}$ for oxygen, and $1s^{1}$ for hydrogen. For the exchange-
correlation energy, we use both the generalized-gradient approximation
functional of Perdew-Burke-Ernzerhof (PBE) PBE  and its modified version for
solids (PBEsol) PBEsol . We find that experimental lattice parameters agree
well with those predicted by PBEsol, and the data presented in the main text
has been obtained using PBEsol (see comparison between PBE and PBEsol in
Supplemenatry Fig. S5a). An on-site Hubbard interaction $U$ is applied to the
copper $3d$ orbitals based on the simplified rotationally invariant DFT+$U$
method by Dudarev and co-workers DFT+U . We have checked that DFT+$U$ gives
almost identical lattice parameters to DFT. Converged results are obtained
with a kinetic energy cutoff for the plane wave basis of $600$ eV and a
$\mathbf{k}$-point grid of size $4\times 4\times 5$ and $6\times 6\times 8$
for the primitive cell of the parent and copper-doped lead apatite,
respectively (see convergence test results in Supplementary Note 1). The
geometry of the structures is optimised until all forces are below $0.01$ eV/Å
and the pressure is below $1$ kbar.

Harmonic phonon calculations. - We perform harmonic phonon calculations using
the finite displacement method in conjunction with nondiagonal supercells
nondiagonal_supercells ; phonon_review . A converged calculation requires a
minimum coarse $\mathbf{q}$-point grid including the high symmetry points
$\Gamma$, M, K, and A, which we accomplish by means of a Farey nonuniform grid
farey_grid  of size $(2\times 2\times 2)\cup(3\times 3\times 1)$. To evaluate
the force derivatives, we use a three-point central formula with a finite
displacement of $0.02$ bohr. The underlying electronic structure calculations
are performed using the same parameters as those described above. We have also
cross-checked the phonon band structures with calculations using castep CASTEP
and Quantum Espresso QE  (see Supplementary Fig. S2).

Anharmonic phonon calculations. - We perform anharmonic phonon calculations
using the stochastic self-consistent harmonic approximation (SSCHA) SSCHA1 ;
SSCHA2 ; SSCHA3 , which accounts for anharmonic effects at both zero and
finite temperature. The self-consistent harmonic approximation SSHA  is a
quantum variational method on the free energy, and the variational
minimization is performed with respect to a trial harmonic system. In its
stochastic implementation, the forces on atoms are calculated in an ensemble
of configurations drawn from the trial harmonic system. We use vasp to perform
electronic structure calculations using the same parameters as those described
above, and consider configurations commensurate with a $2\times 2\times 2$
supercell. The number of configurations needed to converge the free energy
Hessian is of the order of $4{,}000$ for the parent lead apatite structure and
of the order of $8{,}000$ configurations for the copper-doped structure.

##  VII Acknowledgments

###### Acknowledgements.

S.-W.K., K.W., and B.M. are supported by a UKRI Future Leaders Fellowship
[MR/V023926/1]. B.M. also acknowledges support from the Gianna Angelopoulos
Programme for Science, Technology, and Innovation, and from the Winton
Programme for the Physics of Sustainability. S.C. acknowledges financial
support from the Cambridge Trust and from the Winton Programme for the Physics
of Sustainability. The computational resources were provided by the Cambridge
Tier-2 system operated by the University of Cambridge Research Computing
Service and funded by EPSRC [EP/P020259/1] and by the UK National
Supercomputing Service ARCHER2, for which access was obtained via the UKCP
consortium and funded by EPSRC [EP/X035891/1].

##  VIII References

## References

  * (1) Lee, S., Kim, J.-H. & Kwon, Y.-W.  The firs room-temperature ambient-pressure superconductor.  _arXiv preprint arXiv:2307.12008_ (2023).
  * (2) Lee, S. _et al._ Superconductor Pb10-xCux(PO4)6O showing levitation at room temperature and atmospheric pressure and mechanism.  _arXiv preprint arXiv:2307.12037_ (2023).
  * (3) Kumar, K., Karn, N. K. & Awana, V. P. S.  Synthesis of possible room temperature superconductor LK-99: Pb9Cu(PO4)6O.  _Superconductor Science and Technology_ 36, 10LT02 (2023).  URL https://dx.doi.org/10.1088/1361-6668/acf002.
  * (4) Liu, L. _et al._ Semiconducting transport in Pb10-xCux(PO4)6O sintered from Pb2SO5 and Cu3P.  _Advanced Functional Materials_ n/a, 2308938.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adfm.202308938.
  * (5) Wu, H., Yang, L., Xiao, B. & Chang, H.  Successful growth and room temperature ambient-pressure magnetic levitation of LK-99.  _arXiv preprint arXiv:2308.01516_ (2023).
  * (6) Guo, K., Li, Y. & Jia, S.  Ferromagnetic half levitation of LK-99-like synthetic samples.  _Science China Physics, Mechanics & Astronomy_ 66, 107411 (2023).  URL https://doi.org/10.1007/s11433-023-2201-9.
  * (7) Wang, P. _et al._ Ferromagnetic and insulating behavior in both half magnetic levitation and non-levitation LK-99 like samples.  _Quantum Frontiers_ 2, 10 (2023).  URL https://doi.org/10.1007/s44214-023-00035-z.
  * (8) Zhang, Y., Liu, C., Zhu, X. & Wen, H.-H.  Ferromagnetism and insulating behavior with a logarithmic temperature dependence of resistivity in Pb10-xCux(PO4)6O.  _arXiv preprint arXiv:2308.05786_ (2023).
  * (9) Zhu, S., Wu, W., Li, Z. & Luo, J.  First order transition in Pb10-xCux(PO4)6O ($0.9<x<1.1$) containing Cu2S.  _arXiv preprint arXiv:2308.04353_ (2023).
  * (10) Timokhin, I., Chen, C., Yang, Q. & Mishchenko, A.  Synthesis and characterisation of LK-99.  _arXiv preprint arXiv:2308.03823_ (2023).
  * (11) Kumar, K., Karn, N., Kumar, Y. & Awana, V.  Absence of superconductivity in LK-99 at ambient conditions.  _arXiv preprint arXiv:2308.03544_ (2023).
  * (12) Liu, C. _et al._ Phases and magnetism at microscale in compounds containing nominal Pb10-xCux(PO4)6O.  _Phys. Rev. Mater._ 7, 084804 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevMaterials.7.084804.
  * (13) Puphal, P. _et al._ Single crystal synthesis, structure, and magnetism of Pb10-xCux(PO4)6O.  _arXiv preprint arXiv:2308.06256_ (2023).
  * (14) Jain, P. K.  Superionic phase transition of copper (I) sulfide and its implication for purported superconductivity of LK-99.  _The Journal of Physical Chemistry C_ (2023).  URL https://doi.org/10.1021/acs.jpcc.3c05684.
  * (15) Griffin, S. M.  Origin of correlated isolated flat bands in copper-substituted lead phosphate apatite.  _arXiv preprint arXiv:2307.16892_ (2023).
  * (16) Si, L. & Held, K.  Electronic structure of the putative room-temperature superconductor Pb9Cu(PO4)6O.  _arXiv preprint arXiv:2308.00676_ (2023).
  * (17) Lai, J., Li, J., Liu, P., Sun, Y. & Chen, X.-Q.  First-principles study on the electronic structure of Pb10-xCux(PO4)6O ($x=0,1$).  _Journal of Materials Science & Technology_ 171, 66–70 (2024).  URL https://www.sciencedirect.com/science/article/pii/S1005030223006291.
  * (18) Kurleto, R. _et al._ Pb-apatite framework as a generator of novel flat-band CuO based physics, including possible room temperature superconductivity.  _arXiv preprint arXiv:2308.00698_ (2023).
  * (19) Bai, H., Gao, L. & Zeng, C.  Ferromagnetic ground state and spin-orbit coupling induced bandgap open in LK99.  _arXiv preprint arXiv:2308.05134_ (2023).
  * (20) Swift, M. W. & Lyons, J. L.  Comment on ”Origin of correlated isolated flat bands in copper-substituted lead phosphate apatite”.  _arXiv preprint arXiv:2308.08458_ (2023).
  * (21) Pashov, D., Acharya, S., Lany, S., Dessau, D. S. & van Schilfgaarde, M.  Multiple slater determinants and strong spin-fluctuations as key ingredients of the electronic structure of electron-and hole-doped Pb10-xCux(PO4)6O.  _arXiv preprint arXiv:2308.09900_ (2023).
  * (22) Korotin, D. M., Novoselov, D. Y., Shorikov, A. O., Anisimov, V. I. & Oganov, A. R.  Electronic correlations in promising room-temperature superconductor Pb9Cu(PO4)6O: a DFT+ DMFT study.  _arXiv preprint arXiv:2308.04301_ (2023).
  * (23) Si, L. _et al._ Pb10-xCux(PO4)6O: a Mott or charge transfer insulator in need of further doping for (super) conductivity.  _arXiv preprint arXiv:2308.04427_ (2023).
  * (24) Yue, C., Christiansson, V. & Werner, P.  Correlated electronic structure of Pb10-xCux(PO4)6O.  _arXiv preprint arXiv:2308.04976_ (2023).
  * (25) Oh, H. & Zhang, Y.-H.  S-wave pairing in a two-orbital t-J model on triangular lattice: possible application to Pb10-xCux(PO4)6O.  _arXiv preprint arXiv:2308.02469_ (2023).
  * (26) Witt, N., Si, L., Tomczak, J. M., Held, K. & Wehling, T.  No superconductivity in Pb9Cu1(PO4)6O found in orbital and spin fluctuation exchange calculations.  _arXiv preprint arXiv:2308.07261_ (2023).
  * (27) Paudyal, H., Flatté, M. E. & Paudyal, D.  Implications of the electron-phonon coupling in CuPb9(PO4)6O for high-temperature superconductivity: an ab initio study.  _arXiv preprint arXiv:2308.14294_ (2023).
  * (28) Shen, J. _et al._ Phase stability of lead phosphate apatite Pb10-xCux(PO4)6O, Pb10-xCux(PO4)6(OH)2, and Pb8Cu2(PO4)6.  _arXiv preprint arXiv:2308.07941_ (2023).
  * (29) Jiang, Y. _et al._ Pb9Cu(PO4)6(OH)2: Phonon bands, localized flat band magnetism, models, and chemical analysis.  _arXiv preprint arXiv:2308.05143_ (2023).
  * (30) Hao, L. & Fu, E.  First-principles calculation on the electronic structures, phonon dynamics, and electrical conductivities of Pb10(PO4)6O and Pb9Cu(PO4)6O compounds.  _arXiv preprint arXiv:2308.05618_ (2023).
  * (31) Liu, J. _et al._ Symmetry breaking induced insulating electronic state in Pb9Cu(PO4)6O.  _arXiv preprint arXiv:2308.11766_ (2023).
  * (32) Liu, R., Guo, T., Lu, J., Ren, J. & Ma, T.  Different phase leads to different transport behavior in Pb9Cu(PO4)6O compounds.  _arXiv preprint arXiv:2308.08454_ (2023).
  * (33) Cabezas-Escares, J., Barrera, N., Cardenas, C. & Munoz, F.  Theoretical insight on the LK-99 material.  _arXiv preprint arXiv:2308.01135_ (2023).
  * (34) Rooksby, H.  Identification by X-ray diffraction of crystalline inclusions in glass.  _Analyst_ 77, 759–765 (1952).
  * (35) Merker, L. & Wondratschek, H.  Der oxypromorphit Pb10(PO4)6O und der ausschnitt Pb4P2O9-Pb3(PO4)2 des systems PbO-P2O5.  _Zeitschrift für anorganische und allgemeine Chemie_ 306, 25–29 (1960).  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/zaac.19603060105.
  * (36) Merker, L., Engel, G., Wondratschek, H. & Ito, J.  Lead ions and empty halide sites in apatites.  _American Mineralogist: Journal of Earth and Planetary Materials_ 55, 1435–1437 (1970).
  * (37) Krivovichev, S. V. & Burns, P. C.  Crystal chemistry of lead oxide phosphates: crystal structures of Pb4O(PO4)2, Pb8O5(PO4)2 and Pb10(PO4)6O.  _Zeitschrift für Kristallographie-Crystalline Materials_ 218, 357–365 (2003).
  * (38) Chen, S., Salzbrenner, P. T. & Monserrat, B.  Nonuniform grids for Brillouin zone integration and interpolation.  _Phys. Rev. B_ 106, 155102 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.106.155102.
  * (39) Brückner, S., Lusvardi, G., Menabue, L. & Saladini, M.  Crystal structure of lead hydroxyapatite from powder X-ray diffraction data.  _Inorganica Chimica Acta_ 236, 209–212 (1995).  URL https://www.sciencedirect.com/science/article/pii/002016939504636N.
  * (40) Kresse, G. & Furthmüller, J.  Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set.  _Computational Materials Science_ 6, 15–50 (1996).  URL https://www.sciencedirect.com/science/article/pii/0927025696000080.
  * (41) Kresse, G. & Furthmüller, J.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  _Phys. Rev. B_ 54, 11169–11186 (1996).  URL https://link.aps.org/doi/10.1103/PhysRevB.54.11169.
  * (42) Blöchl, P. E.  Projector augmented-wave method.  _Phys. Rev. B_ 50, 17953–17979 (1994).  URL https://link.aps.org/doi/10.1103/PhysRevB.50.17953.
  * (43) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Phys. Rev. Lett._ 77, 3865–3868 (1996).  URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3865.
  * (44) Perdew, J. P. _et al._ Restoring the density-gradient expansion for exchange in solids and surfaces.  _Phys. Rev. Lett._ 100, 136406 (2008).  URL https://link.aps.org/doi/10.1103/PhysRevLett.100.136406.
  * (45) Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P.  Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study.  _Phys. Rev. B_ 57, 1505–1509 (1998).  URL https://link.aps.org/doi/10.1103/PhysRevB.57.1505.
  * (46) Lloyd-Williams, J. H. & Monserrat, B.  Lattice dynamics and electron-phonon coupling calculations using nondiagonal supercells.  _Phys. Rev. B_ 92, 184301 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevB.92.184301.
  * (47) Monserrat, B.  Electron-phonon coupling from finite differences.  _J. Phys. Condens. Matter_ 30, 083001 (2018).  URL https://dx.doi.org/10.1088/1361-648X/aaa737.
  * (48) Clark, S. J. _et al._ First principles methods using CASTEP.  _Zeitschrift für Kristallographie - Crystalline Materials_ 220, 567–570 (2005).  URL https://doi.org/10.1524/zkri.220.5.567.65075.
  * (49) Giannozzi, P. _et al._ QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials.  _Journal of Physics: Condensed Matter_ 21, 395502 (2009).  URL https://dx.doi.org/10.1088/0953-8984/21/39/395502.
  * (50) Errea, I., Calandra, M. & Mauri, F.  Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides.  _Phys. Rev. B_ 89, 064302 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.89.064302.
  * (51) Bianco, R., Errea, I., Paulatto, L., Calandra, M. & Mauri, F.  Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation.  _Phys. Rev. B_ 96, 014111 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevB.96.014111.
  * (52) Monacelli, L. _et al._ The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects.  _Journal of Physics: Condensed Matter_ 33, 363001 (2021).  URL https://dx.doi.org/10.1088/1361-648X/ac066b.
  * (53) Hooton, D.  LI. a new treatment of anharmonicity in lattice thermodynamics: I.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ 46, 422–432 (1955).  URL https://doi.org/10.1080/14786440408520575.

in 1,…,0 See pages \x of SI.pdf