eqs

$\begin{split}\BODY\end{split}$ (1)

††thanks: These authors contributed equally.††thanks: These authors
contributed equally.††thanks: SBL grew crystals and solved the SCXD solutions
reported in this paper.††thanks: These authors contributed equally.††thanks:
These authors contributed equally.††thanks: These authors contributed
equally.††thanks: These authors contributed equally.††thanks: These authors
contributed equally.††thanks: These authors contributed equally.††thanks:
These authors contributed equally.††thanks: These authors contributed
equally.††thanks: These authors contributed equally.

# Pb9Cu(PO4)6(OH)2: Phonon bands, Localized Flat Band Magnetism, Models, and
Chemical Analysis

Yi Jiang  Donostia International Physics Center (DIPC), P. Manuel de
Lardizabal 4, 20018 Donostia-San Sebastian, Spain  Scott B. Lee  Department of
Chemistry, Princeton University, Princeton, New Jersey 08544, USA  Jonah
Herzog-Arbeitman  Department of Physics, Princeton University, Princeton, New
Jersey 08544, USA  Jiabin Yu  Department of Physics, Princeton University,
Princeton, New Jersey 08544, USA  Xiaolong Feng  Max Planck Institute for
Chemical Physics of Solids, 01187 Dresden, Germany  Haoyu Hu  Donostia
International Physics Center (DIPC), P. Manuel de Lardizabal 4, 20018
Donostia-San Sebastian, Spain  Dumitru Călugăru  Department of Physics,
Princeton University, Princeton, New Jersey 08544, USA  Parker S. Brodale
Department of Chemistry and Biochemistry, University of Oregon, Eugene, OR,
97403  Eoghan L. Gormley  Department of Chemistry and Biochemistry, University
of Oregon, Eugene, OR, 97403  Maia G. Vergniory  Donostia International
Physics Center (DIPC), P. Manuel de Lardizabal 4, 20018 Donostia-San
Sebastian, Spain  Max Planck Institute for Chemical Physics of Solids, 01187
Dresden, Germany  Claudia Felser  Max Planck Institute for Chemical Physics of
Solids, 01187 Dresden, Germany  S. Blanco-Canosa  Donostia International
Physics Center (DIPC), P. Manuel de Lardizabal 4, 20018 Donostia-San
Sebastian, Spain  IKERBASQUE, Basque Foundation for Science, 48013 Bilbao,
Spain  Christopher H. Hendon  Department of Chemistry and Biochemistry,
University of Oregon, Eugene, OR, 97403  Leslie M. Schoop  Department of
Chemistry, Princeton University, Princeton, New Jersey 08544, USA  B. Andrei
Bernevig  Department of Physics, Princeton University, Princeton, New Jersey
08544, USA  Donostia International Physics Center (DIPC), P. Manuel de
Lardizabal 4, 20018 Donostia-San Sebastian, Spain  IKERBASQUE, Basque
Foundation for Science, 48013 Bilbao, Spain

###### Abstract

In a series of recent reports, doped lead apatite (LK-99) has been proposed as
a candidate ambient temperature and pressure superconductor. However, from
both an experimental and theoretical perspective, these claims are largely
unsubstantiated. To this end, our synthesis and subsequent analysis of an
LK-99 sample reveals a multiphase material that does not exhibit high-
temperature superconductivity. We study the structure of this phase with
single-crystal X-ray diffraction (SXRD) and find a structure consistent with
doped Pb10(PO4)_6(OH)2. However, the material is transparent which rules out a
superconducting nature. From _ab initio_ defect formation energy calculations,
we find that the material likely hosts OH- anions, rather than divalent
$O^{2-}$ anions, within the hexagonal channels and that Cu substitution is
highly thermodynamically disfavored. Phonon spectra on the equilibrium
structures reveal numerous unstable phonon modes. Together, these calculations
suggest it is doubtful that Cu enters the structure in meaningful
concentrations, despite initial attempts to model LK-99 in this way. However
for the sake of completeness, we perform _ab initio_ calculations of the
topology, quantum geometry, and Wannier function localization in the Cu-
dominated flat bands of four separate doped structures. In all cases, we find
they are atomically localized by irreps, Wilson loops, and the Fubini-Study
metric. It is unlikely that such bands can support strong superfluidity, and
instead are susceptible to ferromagnetism (or out-of-plane antiferromagnetism)
at low temperatures, which we find in _ab initio_ studies. In sum,
Pb9Cu(PO4)6(OH)2 could more likely be a magnet, rather than an ambient
temperature and pressure superconductor.

##  I Introduction

A wave of scientific and social interest has followed a recent claim that
LK-99[1, 2], with the proposed composition of Pb10-xCux(PO4)6O, exhibits
ambient temperature and pressure superconductivity [3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. Although the data presented in
the original reports is insufficient to support such a remarkable claim, a
large body of immediate work has followed. Experimentally, it is unclear what
the structure and composition of the material really are; most likely it is a
multi-phase sample. Experiments have suggested diamagnetic behavior, and
levitation experiments have discussed that it may arise from either
diamagnetism or small ferromagnetic impurities [8]. As the sample is likely to
contain multiple phases, it is possible that two different compounds
contribute to each property, i.e. one part is diamagnetic, and the other
metallic. Thus, clarification of the composition of the material is necessary
_before_ trusting the models based on postulated structures. Assuming some of
the numerous possible compositions, a number of ab initio band structures have
been produced. Claims of flat bands – and their positive influence on
superconductivity – have been made.

However, given the stakes, computation and prediction of physical properties
requires an elevated level of accuracy. While flat bands provide a theoretical
platform for high-temperature superconductivity, non-trivial quantum geometry
is a compulsory for superfluidity, i.e. the Meissner effect[22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32]. This is because the mass of the condensing Cooper
pair is inversely proportional to the minimal Fubini-Study metric [33, 34], a
rigorous measure of quantum geometry. A variety of nonzero lower bounds on the
minimal Fubini-Study metric exist for non-atomic bands [22, 35, 36, 37, 38].
However, flat bands are in fact detrimental to superfluidity if they are
atomically localized, i.e. if their narrow dispersion comes from Wannier
localization as opposed to destructive interference [39, 40, 41]. The
inescapability of this conclusion is seen in the limiting case of a tight-
binding model with all hoppings vanishing. The perfectly flat band that
results cannot support transport of any kind, much less superconductivity,
even in the presence of attractive Hubbard interactions. The key physics of
flat band superconductivity thus lies in deviations from this limit, measured
by quantum geometry [42]. Furthermore, if Coulomb repulsion overwhelms the
would-be attractive interaction, flat band ferromagnetism, rather than
superconductivity, would be favored [43, 44]. Thus a comprehensive
understanding of the active bands at the Fermi level is required for
predictions of the many-body state.

Figure 1: The crystal structure for (a) Pb10(PO4)6O[45] and (b)
Pb10(PO4)6(OH)2[46]. The Wyckoff positions of atoms in space group 176
$P6_{3}/m$ are labelled, with their coordinates given in Table. 1. The O atoms
surrounding P atoms that form PO4 are at $6h$ or generic positions and are not
labeled for simplicity. We also mark two possible Cu doping positions on Pb2
with yellow and blue circles, which are $1b$ and $1c$ Wyckoff positions in
space group 143 $P3$ and will be called Cu1\- and Cu2-doping, respectively.
The H atoms are not shown in (b) for simplicity, which are close to the
trigonal O atoms at $4e$ Wyckoff position.

Due to new advances in understanding and classification of band structures
[47, 48], we can analyze their symmetry, localization, and topology with
extreme accuracy that only depends on the accuracy of the DFT calculation. In
this work, we study two lead apatites (see Fig. 1) available in materials
databases, Pb10(PO4)6O[45] and Pb10(PO4)6(OH)2[46] – which we emphasize may
not be the ultimate material structure. Making assumptions about the location
of the fractionally occupied O anion and the Cu dopant, we obtain several
microscopic models for the electron and phonon bands in this system. Our
results can be summarized as follows.

Firstly in the Pb10(PO4)6O initial structure (Fig. 1(a)), we consider Cu
replacing Pb at two possible positions, either the 1b or 1c Wyckoff positions,
referred to as the Cu1- and Cu2-dopings respectively. Both possibilities
result in a set of two bands at the Fermi level dominated by Cu (see Sec.
III.1). They are narrow with a $\sim 100$-meV bandwidth, and form an
elementary band representation [47] of the Cu $d$-orbitals. We compute their
Fubini-Study metric and non-abelian Wilson loops, which show strong
localization, although the the Cu2 structure does have a significantly reduced
gap to the O bands below. For both possibilities, we obtain 4-band, short-
ranged, symmetric tight-binding models which demonstrate that the weak
dispersion of the Cu orbitals arises primarily from hybridization with nearby
O bands. Secondly we consider the Pb10(PO4)6(OH)2 structure (Fig. 1(b)). For
both locations of the Cu dopant, we again find a set of Cu bands in an
elementary band representation with $\sim 100$-meV bandwidth. However, the gap
to the nearby O bands is much larger, and we provide a 2-band model built
entirely from Cu Wannier functions. Again, the Fubini-Study metric indicates
atomic localization. See Sec. III and Sec. IV for more details.

In all cases, the two bands at the Fermi level lack strong quantum geometry.
Due to their flat, localized nature, ferromagnetism seems to be the preferred
configuration of these states in ab initio studies (see Sec. III.1). The
absence of extended states in these bands does not support a theory of high-
temperature superconductivity based on the flat bands we obtain in these
structures. However, preliminary calculations of the phonon spectrum show that
more careful relaxation of the doped compound is required to fully stabilize
their structures, which may result in changes to the band geometry. A phonon-
driven mechanism for superconductivity must also compete with the strong
Hubbard repulsive interaction, which we also construct using ab initio
Hubbard-Kanamori parameters.

In Refs. [1, 2], LK-99 shows a sharp drop of the resistivity around 400 K,
towards a state claimed to be superconducting. Nevertheless, the reported
value of resistivity is 2-3 orders of magnitude higher than that of good
metals; for instance, Cu presents a resistivity value of 10-6 $\Omega$cm.
Moreover, the analysis of the reported specific heat shows no transition up to
400 K, in principle, at odds with a jump expected due to the release of
entropy of a second order phase transition. Indeed, the specific heat seems to
drop with temperature, which adds more controversy to the claim of room-
temperature superconductivity. Ref. [49] has recently found evidence for the
Cu2S present in the multi-phase compound sourcing the resistivity transition.
Moreover, the presence of a diamagnetic signal does not necessarily imply SC.
In Fig. 2(A), we show the magnetic response of the diamagnetic compound HOPG
(Highly Oriented Pyrolityc Graphite). As we can see, the diamagnetism of HOPG
is clearly different from any superconductor [50].

Figure 2: (A) Magnetization vs field for (HOPG) Highly Oriented Pyrolityc
Graphite, showing a diamagnetic response. (B) Hystereses in the curves of the
torque ($\tau$) vs H for NbSe2.

On the other hand, besides the typical hysteric behavior of the magnetization
expected in a SC, in Fig. 2B, we show the magnetic torque response
$\boldsymbol{\tau}=\mathbf{M}\times\mathbf{B}$ of the type II SC NbSe2. The
hysteric behavior between the sweep-up and -down curves defines the different
field regions of a type II SC (vortex solid, vortex liquid, critical field,
etc) [51].

##  II Chemical Structure

Pb-apatite is structurally related to the parent compound, Ca-apatite, and
much can be learned through their comparison. This crystal type features
1-dimensional channels filled with requisite charge balancing anions. For
example, the calcium-oxo analogue has been shown to only exist as Ca10(PO4)6O
at temperatures exceeding 1000${}^{o}C$,[52] below of which it rapidly
includes water into the lattice, forming the neutral hydrate Ca10(PO4)6(OH)2.
While the Pb-based system may have different dehydration temperatures to
interconvert between Pb10(PO4)6O[45] and Pb10(PO4)6(OH)2[46] (see Fig. 1),
historical data would indicate that the hydrated material is thermodynamically
preferred. To this end, we performed a number of DFT calculations to assess
the formation enthalpy of the hydroxy and oxo Pb apatites considered here. Our
calculations reveal that the inclusion of H2O, forming Pb10(PO4)6(OH)2, is
exothermic with its inclusion favored by 38.5 kcal/mol. These data suggest
while Pb10(PO4)6O may be a metastable form, direct syntheses of that material
will tend to form the hydrate if sufficient protons are available. And even if
Pb10(PO4)6O is formed it will interconvert to the hydrate upon exposure to
air.

Structurally, the oxo and hydroxy apatites are better distinguished by their
lattice parameters, and less so by direct crystallographic measurements. Ca-
apatite exhibits a dramatic lattice contraction in both the a and c
parameters[53], proportional to the extent of hydration, i.e. the conversion
from O2- to OH- contracts the lattice.[54]. Such contraction has been
attributed to Cu-inclusion within the Pb10(PO4)6O lattice[2, 7], but a similar
effect is predicted through simple hydration.

While theoretical studies thus far have focused on the pure Pb10-xCux(PO4)6O
or Pb10-xCux(PO4)6(OH)2 phase, the reported synthesis methods cannot result in
single phase samples of either of them. In the final step of the reaction,
Lanarkite (Pb2(SO4)O) is mixed in a 1:1 molar ratio with Cu3P, but this 2:1
Pb/P ratio is stoichiometrically inconsistent with the desired product. Even
loosening the ratio of the reactants would create a significant amount of Cu
impurities: the final product has 6 P, so there will be 18 Cu per formula unit
of Pb9Cu1(PO4)6O or Pb9Cu1(PO4)6(OH)2. The powder X-ray diffraction (PXRD)
pattern presented in the original reports[1, 2] suggests that the majority of
the sample exists in a structure related to Pb10(PO4)6(OH)2, but also features
many crystalline impurities. Due to the severe off-stoichiometry of the
reaction, additional amorphous phases are also to be expected. Only few
Rietveld analyses have thus far been performed[12], and it is difficult to
separate the structural effects of Cu inclusion versus hydration without high-
quality single crystal data. To this end, is unclear if any Cu atoms were
actually incorporated into the structure.

To address this, we are able to compute the formation enthalpies of Cu
substitutions at Pb lattice sites. Following the standard procedure[55], the
Cu defects were computed using a 2$\times$2$\times$2 super cell containing
1831 electrons for the neutral substitution. Computations were referenced to
bulk Cu/Pb (rich potentials) and CuO/PbO (poor potentials) and the formation
enthalpy was computed using a 2$\times$2$\times$2 k-mesh, with the FNV
correctional scheme employed for charged defects and averaged diagonalized
dielectric tensors[56, 57], and Hubbard model (U = 4.0 eV for Cu). The
formation enthalpies for the defects are depicted in 3. Notably, the
experimentally observed Cu2+ (i.e. charge-neutral) substitution is predicted
to form in p-type conditions, at minimum occurring with a 1.2 eV penalty
depending on the reference potentials and amounting to many orders of
magnitude sub-stoichiometric Cu concentration, essentially forbidding its
inclusion within the material. However, the charge transition from inclusion
of Cu2+ to Cu1+ occurs very near the Fermi level for the undoped Pb-apatite
material, and Cu1+ inclusions becomes more favored in n-type conditions. The
experimental absence of appreciable Cu1+ suggests that if Cu2+ is
incorporating, it is being aided by other correlated defects not considered
here. Additionally, other mechanisms of Cu inclusion are possible (e.g.
interstitial formation, as observed in Pb2+-containing in lead halide
perovskites[58]). Pb-apatite may host other dopants e.g. S, which is also
present during the reaction[59]. These may be studied in future work.

Figure 3: The defect formation enthalpies for Cu substitution in Pb lattice
sites at both rich and poor potentials show that CuPb is unfavorably
incorporated into either Pb site within the apatite lattice.

To address whether Cu and OH- are included, and to generally deduce the
structure, we performed a modified synthesis of the putative LK-99 compound
described in Ref.[1, 2]. Sample purity of all intermediate reagents was
confirmed via powder X-ray diffraction using a STOE Stadi P powder X-ray
diffractometer equipped with a Mo K$\alpha$ ($\lambda$=0.71073 Å) sealed-tube
x-ray source and graphite monochromator at room temperature in either Debye-
Scherrer or transmission geometry Fig. 25, Fig. 26. In the final step, Cu3P
and Pb2(SO4)O in a 1:1 molar ratio were ground into a fine powder in mortar
and pestle. The powders were loaded into an alumina crucible, placed in a
quartz tube, and sealed under dynamic vacuum without any argon back-filling.
The tube was then loaded in a furnace and heated to 950 ∘C over 4 hours, kept
at this temperature for 6 hours, and then shut off to cool quickly. (Upon
pulling the sample out of the furnace we observed that our sample, due to its
loading in an alumina crucible, did not attack our quartz tube in contrast to
Refs. [1, 2]. This suggests the possibility of further chemical reactions in
the protocol of Refs. [1, 2].) The resulting product inside the crucible
contained mostly white and orange colored powder in addition to metallic gray
lumps as seen in Fig. 28. A PXRD pattern of this sample can be seen in Fig.
27, showing that our sample has more impurity phases than the original LK-99
sample, but that the characteristic peaks are present and thus it seems
possible to isolate the main phase reported in the original sample from our
sample as well. Energy-dispersive X-ray spectroscopy (EDS) was also utilized
to examine these impurity phases. In Fig. 29, we show that the transparent-
orange phase likely gets its color from pockets of metallic Cu. Likewise, in
Fig. 30 we observe islands of Cu2S embedded in a matrix of the transparent-
white phase. As noted previously[60, 61], these Cu2S islands are likely the
cause of the tenfold drop in resistivity observed at 104.8 ∘C in the original
preprint.

Under magnification, the white powder was observed to be transparent single
crystals with a clear hexagonal rod habit (Fig. 4). Knowing that Pb10(PO4)6O
was reported to crystallize in a hexagonal space group, a transparent single
crystal of dimensions (0.55 $\times$ 0.092 $\times$ 0.176 mm) was picked for
single crystal X-ray diffraction (SXRD) analysis using an APEX2 CCD
diffractometer equipped with a Mo K$\alpha$ ($\lambda$=0.71073 Å) sealed-tube
X-ray source and graphite monochromator at room temperature. Initial unit cell
refinement obtained a unit cell in a hexagonal setting with lattice parameters
a=9.8508(1) Åand c=7.4395(2) Å, agreeing with other previously reported
structures for lead apatite [46, 45]. Indexation and integration proceeded
smoothly for a full hemisphere collection out to a resolution of 0.5 Å. Run
list generation and frame data processing were done in APEX 2 [62]. An
analytical absorption correction was used to scale the data before importing
the peak list into JANA2020 [63].

Figure 4: A single crystal similar in transparency, shape, and size to what we
picked for SXRD.

The initial structure solution was obtained in space group P63/$m$, consistent
with observed systematic absences Fig. 31, from the charge-flipping algorithm
as implemented in SUPERFLIP [64]. An initial refinement, done on $F^{2}$, was
needed to obtain a valid initial structure with 6 distinct sites (2 Pb, 1 P,
and 3 O). Notably, during this initial refinement Oxygen atoms along the
c-axis were removed. Running a Fourier transform of residual electron density,
visualized with the Vesta software package [65], reveals 2 unique
crystallographic pockets of significant electron density around points [0 0 0]
(Wyckoff position 2$b$) and [0 0 $1/4$] (2$a$) in P63/$m$, shown in Fig. 5,
suggestive of a Pb10(PO4)6(OH)2 structural solution. After Oxygen atoms were
added to these points, freely refining their occupancies results in a site
occupancy factor of the 2$b$ Oxygen centered at [0 0 0] of 1.163, an
unphysical value for an OH- molecular unit. Clearly, an atom with more
electron density is needed. Replacing the oxygen instead with sulfur, a
chemically similar element that is also present in the reaction, seems like a
likely candidate, although further elemental analysis methods will be needed
to distinguish between an SH-, PH${}_{2}^{-}$ or other possible dopants at
this site.

Figure 5: Residual electron density located at the 2a and 2b positions along
the c axis. Positive isosurfaces are drawn at 7e Å-3.

Freely refining occupancies for both an O at [0 0 $1/4$] and S located at [0 0
0] results in site occupancy factors of 0.54 and 0.56, respectively. The
freely refined occupancies having a summation very close to 1 seems suggestive
that this could be a physically real interpretation of the structural
solution. After refining the isotropic parameters anisotropically, a final
refinement converges with a goodness of fit (GoF) parameter of 1.22 and R
value of 3.94 compared to all reflections. Here, our site occupancies for O
(1.08) and S (1.12) sum to above 1. If we decide to restrict the occupancies
of the 2$b$ S and 2$a$ O to sum to a site occupancy factor of 1, we end
refinement with a final composition Pb10(PO4)6(OH)0.94(SH)1.06. This
refinement has nominally the same refinement statistics (GoF(all) $=1.22$,
R(all) $=3.95$).

We also investigated two ways of possible Cu doping within our refinement: one
in which the Cu atom substitutes on the Pb lattice site, and another which the
Cu atom is inserted along the chain. Attempts to dope both Pb sites with Cu
results in a refined composition of Pb9.55Cu0.45(PO4)6(OH)0.94(SH)1.06 with
similar statistics (GoF(all)$=1.21$, R(all) $=3.93$). We found the Cu has to
be added in as a split site after Pb position has been refined completely. The
last two steps of occupancy and anisotropic parameter refinement had to be
done by restraining the Pb/Cu split site position with automatic refining keys
switched off.

Figure 6: Energy dispersive X-ray Spectroscopy map of a translucent crystal
coated in carbon.

The compositional refinement of this structure certainly warrants skepticism.
Unfortunately for us, this system has complications for EDS, namely, the Pb
Lα1 edge lies at the same energy as the S Kα1 edge. This, along with neutron
diffraction being needed to confirm H incorporation means future studies are
needed. However, in search of confirming the refinements Cu incorporation, we
carbon coated a small transparent needle and ran EDS on a Quanta environmental
scanning electron microscope. Using the elemental mapping technique, we were
surprised to find that we were able to detect Cu, and it is homogeneously
distributed within the crystal. Furthermore, in Fig. 6 the Pb:Cu ratio we
detected over a roughly 15 minute period was 58.3:2.8, which roughly
translates to a ratio of 9.55:0.46, an indicator that the Cu incorporation in
our refined composition of Pb9.55Cu0.45(PO4)6(OH)0.94(SH)1.06 could be
correct. However, we would like to acknowledge the fact that 2.8 atomic
percent is extremely low for EDS characterization and future studies utilizing
synchrotron radiation are needed to confirm incorporation. In addition, as
there is a lot of Cu in the reaction it cannot be ruled out that Cu is at the
surface of the crystal.

We can also obtain a reasonable refinement to the data if we place Cu into the
channel with composition Pb10(PO4)6(OH) 1.11 Cu0.49 (GoF $=1.22$, R $=3.93$).
Still, as the crystals are transparent, charge balance needs to be maintained
and thus this last solution is not chemically reasonable. Full collection,
integration, and refinement statistics can be found in App. C. We would to
stress here again that the transparent nature of the crystal should rule out
superconductivity as a property, it rather indicates a wide band gap.

Next, we compare simulated powder patterns of our structural solution, the
ICSD reported structures, and the relaxed DFT structures of Cu doped variants
with the published PXRD pattern from [1, 2]. To do this, structural position
files were loaded into VESTA, and patterns simulated for using a Cu K$\alpha$
($\lambda=1.5406\text{\AA}$) wavelength. Simulated patterns were then overlaid
atop the experimental pattern, extracted using an in-house Mathematica code.
For ease of visualization, we adjust the simulated patterns via a zero point
shift to match the peaks expected around 18 degrees in the experimental data.
This zero-point correction is not uncommon for Rietveld refinements, and is
needed if the diffractometer used in the experimental pattern is misaligned.
Most zero-point corrections were minimal, i.e. $0.2^{\circ}$. We find a good
agreement (see Fig. 32 and Fig. 33) of our SXRD solution with the reported
data. Other structures fit the data less well and a discussion is given in the
SI. This analysis is preliminary as the cropped data from the original
preprint are not of high enough quality to perform a Rietveld analysis. It
does however show that our structure obtained from SRXD, measured in a
_transparent_ crystal, agrees with the powder pattern published in the
original LK-99 paper. Further analysis of our samples will follow.

Due to the numerous uncertainties regarding the ultimate structural
composition of the LK-99 material, we will investigate a variety of scenarios
using ab initio density functional theory. In our ab initio studies, we still
focus on structures on which Cu substitutes for Pb, as we come to different
conclusions as previous theoretical works assuming the same substitutions.

Figure 7: Electronic band structure for (a) Pb10(PO4)6O and (b)
Pb10(PO4)6(OH)2, and zoom in plots in (c) for Pb10(PO4)6O and (d) for
Pb10(PO4)6(OH)2. The quasi-flat bands near $E_{f}$ are mainly formed by the
$(p_{x},p_{y})$ orbitals of the trigonal-O atoms, which are shown with blue
weights.

##  III ab initio Results

Compound | Atom | Wyckoff position
---|---|---
Pb10(PO4)6O | Pb1 | $6h$
Pb2 |  $4f$, $z=0.004$
P | $6h$
tri-O |  $4e$, $z=0.134$, $\frac{1}{4}$-occu
Pb10(PO4)6(OH)2 | Pb1 | $6h$
Pb2 |  $4f$, $z=0.994$
P | $6h$
tri-(OH)2 |  $4e$, $z=0.040$, $\frac{1}{2}$-occu
Table 1: The atomic positions in Pb10(PO4)6O[45] and Pb10(PO4)6(OH)2[46]. The
lattice constants for Pb10(PO4)6O are $a=9.865,c=7.431$ Å, and for
Pb10(PO4)6(OH)2 are $a=9.866,c=7.426$ Å. The Wyckoff positions in SG 176 have
the following coordinates:
$4f=(\frac{1}{3},\frac{2}{3},z),(\frac{2}{3},\frac{1}{3},z+1/2),(\frac{2}{3},\frac{1}{3},-z),(\frac{1}{3},\frac{2}{3},-z+1/2)$,
$4e=(0,0,\pm 1),(0,0,\pm 1+\frac{1}{2})$, and
$6h=(x,y,\frac{1}{4}),(-y,x-y,\frac{1}{4}),(-x+y,-x,\frac{1}{4}),(-x,-y,\frac{3}{4}),(y,-x+y,\frac{3}{4}),(x-y,x,\frac{3}{4})$.
For both compounds, the Pb2 atoms at $4f$ approximately lie on the honeycomb
lattices on $z=0,\frac{1}{2}$ planes. The ‘tri’ in the table stands for
trigonal lattice. The tri-O and tri-OH in the two compounds have fractional
occupancies, i.e., $\frac{1}{4}$ for tri-O and $\frac{1}{2}$ for tri-(OH)2.

We consider two different experimental structures of lead apatite for ab-
initio calculations in this work, i.e., Pb10(PO4)6O[66, 45] and
Pb10(PO4)6(OH)2[67, 46], both have space group (SG) 176 $P6_{3}/m$ symmetry.
Their crystal structures are shown in Fig. 1. The O atoms at $4e$ Wyckoff
position in Pb10(PO4)6O have $\frac{1}{4}$ occupancy, while (OH)2 atoms at
$4e$ in Pb10(PO4)6(OH)2 have $\frac{1}{2}$ occupancy. The atomic positions are
summarized in Table. 1. We remark that lead apatite has many experimentally
reported structures, with the position of trigonal-O or (OH)2 being slightly
different. For example, the structures in Ref.[68, 69, 70, 71] has
trigonal-(OH)2 located at Wyckoff position
$2a=(0,0,\frac{1}{4}),(0,0,\frac{3}{4})$ without fractional occupancy. The
LK-99[1, 2] Pb10-xCux(PO4)6O ($0.9<x<1.1$) is hypothesized to be synthesized
by doping Pb atoms with Cu at $4f$ Wyckoff positions.

Figure 8: Electronic band structure for Pb9Cu1(PO4)6O, where (a) is the PM
phase and (b) is the FM phase of Cu1-doping. (c) and (d) are similar for
Cu2-doping. The atomic positions of two Cu-doping can be found in Table. 5.

###  III.1 Electronic Structure

In this section, we discuss the electronic structures for Pb10(PO4)6O and
Pb10(PO4)6(OH)2 in both undoped and Cu-doped phases.

We first consider the undoped phase. As reported in literature[45, 46], the
trigonal-O in Pb10(PO4)6O has $\frac{1}{4}$ occupancy, while the
trigonal-(OH)2 in Pb10(PO4)6(OH)2 has $\frac{1}{2}$ occupancy. Such fractional
occupancy is difficult to treat in DFT. For simplicity, we fix their positions
to remove the fractional occupancy, i.e., fix O at $(0,0,0.634)$, and fix O in
OH2 at $(0,0,0.04),(0,0,0.54)$ and H at $(0,0,-0.10),(0,0,0.40)$, s.t the
total number of electrons in the unit cell is the same as the fractionally
occupied structure. Mind that the positions of H atoms are not given in the
original experimental structure[46] we use and are manually added using the
$O$-$H$ bond length in H2O molecular, i.e., about 1Å, which also agrees with
the $O$-$H$ length reported in another experiment structure[72]. Notice that
after fixing the trigonal-O and (OH)2 positions, the original SG 176
$P6_{3}/m$ symmetry is lowered to SG 143 $P3$ and SG 173 $P6_{3}$ for two
structures, respectively. Remark that fixing trigonal-O at any of the four
$4e$ positions are equivalent, as they are related by
$\\{C_{6z|00\frac{1}{2}}\\}$ and $\\{M_{z}|000\\}$ in SG 176. For trigonal
(OH)2 the scenario is similar. Thus we will focus on the aforementioned
positions of trigonal-O and (OH)2 in the following.

We relax the structure without fractional occupancy and obtain the relaxed
lattice constants and atomic positions summarized in Appendix. A.1 (Table. 4,
5). The symmetry is maintained during the relaxation, i.e., SG 143 for
Pb10(PO4)6O and SG 173 for Pb10(PO4)6(OH)2. The band structures for two
relaxed structures are shown in Fig. 7, which host a large band gap of 2.761
and 2.635 eV, respectively. The highest occupied bands are quasi-flat and
mainly come from the $(p_{x},p_{y})$ orbitals of the trigonal-O atoms.

We then consider the Cu-doped phases. As hypothesized in Ref.[1, 2], the Cu-
doped LK-99 Pb10-xCux(PO4)6O ($0.9<x<1.1$) has Cu doping the Pb atoms at $4f$
position. By fixing $x=1$ for simplicity, there exist four possible Cu-doping
structures by placing the Cu atom at one of the $4f$ positions. The four $4f$
positions are equivalent in SG 176. However, after fixing (OH)2 in
Pb10(PO4)6(OH)2, the SG is lowered to SG 173 $P6_{3}$, which gives two
inequivalent Cu doping positions, which we call the one at
$(\frac{1}{3},\frac{2}{3},z)$ Cu1-doping and $(\frac{2}{3},\frac{1}{3},-z)$
Cu2-doping, as marked using yellow and blue circles in Fig. 1. For
Pb10(PO4)6O, after fixing the position of trigonal-O, the SG is lower to 143
$P3$ which makes all 4 positions inequivalent. For simplicity, we only
consider the Cu1-doping and Cu2-doping, as the other two possible dopings have
similar band structures as reported in Ref.[14].

Figure 9: Electronic band structure for Pb9Cu1(PO4)6(OH)2, where (a) is the PM
phase and (b) is the FM phase of Cu1-doping. (c) and (d) are similar for
Cu2-doping. The atomic positions of two Cu-doping can be found in Table. 5.

We relax the structure of Cu1\- and Cu2-doping for Pb10(PO4)6O and
Pb10(PO4)6(OH)2, both having SG 143 $P3$ symmetry, with lattice constants and
atomic positions summarized in Appendix. A.1 (Table. 4, 5). The four relaxed
structures are used to perform DFT calculations. We consider the paramagnetic
(PM) phase and ferromagnetic (FM) for each structure, with their band
structures shown in Fig. 8, 9, and the orbital projections given in Appendix.
A.2 (Fig. 13, 14).

In the PM phase, for all four structures, there exist two quasi-flat bands
with $\frac{3}{4}$ filling at the Fermi level $E_{f}$, contributed mainly by
the $(d_{xz},d_{yz})$ orbitals of Cu (with Cu $(d_{xz},d_{yz})$ weight about
$50\%$, $(d_{xy},d_{x^{2}-y^{2}})$ about $20\%$, and $p$ orbital of O atoms
that close to Cu weight about $30\%$), forming elementary band representation
(EBR) ${}^{1}E^{2}E@1b$ in SG P3 (notice in the PM phase the time-reversal
symmetry exists and enforces ${}^{1}E$ and ${}^{2}E$ to be degenerate) for Cu1
doping and ${}^{1}E^{2}E@1c$ for Cu2 doping. For two Cu dopings of
Pb10(PO4)6O, there are two bands close the to Cu bands with a small band gap,
mainly formed by the bands of $(p_{x},p_{y})$ of the trigonal-O and form EBR
${}^{1}E^{2}E@1a$, as shown in Fig. 8, 13. Although here these four bands
formed by $(d_{xz},d_{yz})$ of Cu and $(p_{x},p_{y})$ of the trigonal-O do not
have band inversions and are topologically trivial, we can adjust the hoppings
and make them topological, which we show in Sec. IV using minimal tight-
binding models. For two Cu dopings of Pb10(PO4)6(OH)2, however, the bands
below the Cu bands have a larger band gap and weak hybridization with Cu
bands.

The FM phases of the four structures all have a lower total energy compared
with the PM phase, as shown in Table. 2, suggesting that the FM phase is
energetically more favored than the PM phase. Evidence for soft
ferromagnetism, consistent with the well-localized bands and strong Coulomb
repulsion we find, has recently appeared in experiments [8]. In the spin-
polarized orbital projections shown in Fig. 13, 14, the quasi-flat bands of Cu
near $E_{f}$ have a large spin-splitting of about $0.6$ eV, while the bands of
trigonal-O have negligible spin-splitting. The calculated magnetic moment is
about 0.6 $\mu_{B}$ for Cu and small moments for O atoms surrounding Cu which
sum to about 0.3 $\mu_{B}$ in total, in agreement with the fact that the two
flat bands near $E_{f}$ in the PM phase also have about $30\%$ weight from O.
We also calculate the total energy of the A-type antiferromagnetic phase
(AFM), i.e., FM inplane and AFM out of the plane, as shown in Table. 2. The FM
and AFM phases have almost the same energy per unit cell, both being lower
than the PM phase. We leave for future studies to investigate the magnetic
ground state.

Compound | Phase | PM | FM | AFM
---|---|---|---|---
Pb10(PO4)6O | Cu1-doping | 0 | -0.129 | -0.130
Cu2-doping | 0.195 | 0.089 | 0.092
Pb10(PO4)6(OH)2 | Cu1-doping | 0 | -0.136 | -0.136
Cu2-doping | -0.177 | -0.332 | -0.332
Table 2: Comparison between total energy of PM, FM, and A-type AFM phases per
unit cell (in eV) calculated in DFT. The total energies in the Cu1-doping
phase in Pb10(PO4)6O and Pb10(PO4)6(OH)2 are used as zero. It can be seen that
in all four cases, the FM and AFM phases have very close total energy, which
are both lower than the PM phase. In Pb10(PO4)6O the Cu1 doping has lower
energy while Pb10(PO4)6)(OH)2 the Cu2 phase is lower.

###  III.2 Phonon Spectrum

(a) Pb-(OH)2 at low-T.

(b) Pb-(OH)2 at high-T.

(c) Pb-(OH)2 at high-T in 112 supercell.

Figure 10: Phonon spectrum for relaxed Pb10(PO4)6(OH)2 (without Cu doping)
structure at (a) low and (b) high temperature calculated using $1\times
1\times 1$ supercell, which shows imaginary, unstable phonons in the
$k_{3}=\pi$ plane. (c) Phonon spectrum calculated using a $1\times 1\times 2$
supercell at high temperature, which stabilizes the $k_{3}=\pi$ modes but
softens a single mode at the $\Gamma$ point. In addition, we observe that the
soft phonon modes shown in panels (a) and (b) are flat modes in the
$k_{3}=\pi$ plane. A series of flat phonon modes are also presented at the
finite frequency in the high-$T$ calculations, such as flat modes around 1 THz
in panels (c).

In this section, we perform the phonon calculations to check the stability of
proposed structures in the literature, including Pb10(PO4)6O, Pb10(PO4)6(OH)2,
Pb9Cu1(PO4)6O and Pb9Cu1(PO4)6(OH)2. We find that, in their nominal 111 unit
cells, _all_ structures display imaginary phonon frequencies. We focus our
discussion below on Pb9Cu1(PO4)6(OH)2 since powder X-ray diffraction shows a
strong similarity between it and LK-99. Supplemental figures for other
structures may be found in App. A.4.

For the undoped structures, both the Pb-O and Pb-(OH)2 present
negative/imaginary phonon at low-T and high-T with $1\times 1\times 1$ cell,
where the phonon of Pb-(OH)2 is shown in Fig. 10(a,b). Here, the temperature
effect is included via electronic smearing in the harmonic approximation
level. As the atomic projection presents, the phonon instability is mainly
contributed - as theoretically expected - by the heaviest Pb atoms at both Pb1
and Pb2 sites. As the temperature goes higher, the imaginary phonons harden,
especially in the $k_{z}=0$ plane where imaginary modes disappear at high-T.
Since the calculation is performed in $1\times 1\times 1$ cell with a short
$c$ compared to $a$, the negative phonon in $k_{z}=\pi$ plane may be caused by
the short cutoff in the $c$ direction. To verify this scenario, we also
perform phonon calculation in $1\times 1\times 2$ supercell at high-T. As
shown in Fig. 10(c), the negative branches on the $k_{z}=\pi$ plane in the
$1\times 1\times 1$ cell calculation become positive, as the cell is enlarged
to include the force constants between atoms with longer distance. We expect
that the residual imaginary phonon at the $\Gamma$ point can be eliminated by
performing phonon calculations in a larger supercell after relaxation.

We also observe that the phonon spectrum presents a good separation of
frequency based on the mass of elements: the heaviest Pb dominates the lowest
frequency and gives the imaginary modes, H phonons lie at a much higher
frequency, which is not shown in the plot, and P and O phonons lie in the
middle.

For the Cu-doped structure, previous DFT calculations suggest a
(ferro/antiferro)magnetic ground state. Therefore, the phonon calculations are
performed in both paramagnetic and ferromagnetic phases as presented in App.
A.4 with a $1\times 1\times 1$ cell. Similar to the undoped case, the doped
structures show negative phonon modes, but tend to harden in the FM phase.
Both O and the Cu dopant contribute to the imaginary phonon at low $T$ (see
App. A.4 ) and one should relax the structure with much care for the doping
effect to obtain a stable phonon spectrum. The difficulty in obtaining a
stable structure even for the undoped parent compounds, which exist in nature,
underscores an important challenge for first principles studies of the doped
compound, whatever its nature.

##  IV Tight-Binding Models

In this section, we construct short-range tight-binding models for both
Pb9Cu1(PO4)6O and Pb9Cu1(PO4)6(OH)2 compounds and for both considered
positions of the Cu dopant. Such models are a prerequisite for studying the
many-body phases that LK-99 is conjectured to realize. We emphasize that the
precise chemical composition, purity, and structure of the supposed compound
are far from being settled, and initial proposals may need to be re-examined.
Nevertheless, we study the _proposed_ structures here for completeness. We
find, in all cases and in agreement with independent calculations, that Cu
forms a high-density flat at the Fermi level. Based on this finding we ask
whether a flat band superconductivity scenario is viable. We perform the
calculations of the quantum geometry in these bands and find that they are
atomically localized and will likely not favor superconductivity.

In all cases, we find that the two bands straddling the Fermi level are atomic
and dominated by the Cu $d$-orbitals. In the Pb9Cu1(PO4)6O structure, a four-
band model is required due to the close proximity of the O bands, which
hybridize with the Cu orbitals particularly in the $k_{z}=\pi$ plane. For this
structure, the dispersion of the Cu bands is predominantly due to
hybridization with O. In contrast, Pb9Cu1(PO4)6(OH)2 shows well-isolated Cu
bands at the Fermi level, and a two-band model can be constructed. App. B
contains a complete discussion of the parameters, symmetry, and quantum
geometry of the bands.

A low-energy model of the flat bands can be constructed from the
$d_{yz},d_{xz}$ orbitals of Cu on the 1b/1c positions for Cu1/Cu2 doping, and
the $p_{x},p_{y}$ orbitals of O on the 1a position in space group $P3$ (see
App. B for conventions of the lattice). Our DFT calculations show degeneracies
at $\Gamma$ and $A$ (which form double Weyl points [9]) that cannot be
protected by $C_{3}$ alone since it forms an abelian group. We check that the
separate spin sectors possess a spin-less time-reversal symmetry
$\mathcal{T}^{2}=+1$ which protects a 2D complex irrep ${}^{1}E{}^{2}E$
arising from the $d$ and $p$ orbitals. The preservation of the spin-less
$\mathcal{T}$ in the FM phase comes from the fact that the magnetization is
taken into account as the local momentum term $M(\bm{r})s_{z}$ (with real
scalar $M(\bm{r})$). The spin-orbital coupling is negligible for the Cu and O
atoms.

###  IV.1  Pb9Cu(PO4)6O

Figure 11: Comparison of DFT and tight-binding model band structures for
Pb9Cu(PO4)6O. The relaxed DFT (red) and short-range tight-binding model (blue)
are shown for Cu1 doping $(a)$ and for Cu2 doping $(b)$. The 2D quantum metric
$g(\mathbf{k})$ is computed on the plane $k_{3}=\pi$ for the Cu1 (c) and Cu2
(d) models, showing peaked features where the Cu and O bands have the smallest
direct gap.

In the Pb9Cu(PO4)6O structure, it is necessary to construct a four-band model
due to the O-dominated bands that appear closely below (and hybridize with)
the Cu-dominated bands at the Fermi level. The model takes the form

$\displaystyle
h_{4}(\mathbf{k})=\begin{pmatrix}h_{C}(\mathbf{k})&h_{CO}(\mathbf{k})\\\
h^{\dagger}_{CO}(\mathbf{k})&h_{O}(\mathbf{k})\end{pmatrix}$ (2)

consistent with these symmetries. For Cu1 doping, it is sufficient to include
the three nearest-neighbor Cu-O hoppings ($h_{CO}(\mathbf{k})$), six O-O in-
plane hoppings for O and Cu, and the two O-O vertical hoppings
($h_{O}(\mathbf{k})$). In fact, the Cu hoppings are less than 4meV and can be
safely dropped, effectively setting $h_{C}(\mathbf{k})=0$. The Cu2 doping
structure exhibits a smaller gap between the O and Cu bands, and requires the
inclusion of next nearest neighbor hoppings, with full expressions given in
App. B. We find excellent agreement with the DFT spectrum and irreps in Fig.
11 within this approximation.

Our tight-binding model shows that the dispersion of the Cu bands at the Fermi
level arises essentially in its entirety from hybridization with O.
Nevertheless, symmetry eigenvalues indicate that the bands are topologically
trivial: O hybridization does not cause a topological change. To confirm this,
we compute the 2D Fubini-Study quantum metric (Fig. 11(c,d))

$\displaystyle g(\mathbf{k})$
$\displaystyle=\frac{1}{2}\sum_{i=x,y}\mathop{\mathrm{Tr}}\partial_{i}P(\mathbf{k})\partial_{i}P(\mathbf{k})$
(3)

where $P(\mathbf{k})=U(\mathbf{k})U^{\dagger}(\mathbf{k})$ is the projector
onto the eigenvector matrix $U(\mathbf{k})$ of the two Cu bands, and non-
abelian Wilson (Fig. 20) loop

$\displaystyle W(k_{2},k_{3})$
$\displaystyle=U^{\dagger}(2\pi,k_{2},k_{3})\prod_{k_{1}}^{2\pi\leftarrow
0}P(\mathbf{k})\,U(0,k_{2},k_{3})$ (4)

computed over the set of two bands at the Fermi level. Both show strongly
localized states, consistent with their elementary band representation.
Furthermore, the basis Wanniers functions are localized: square root of
Wannier spread is about $0.20a\approx 0.26c$ for Cu and about $0.33a\approx
0.43c$ for O). As we can see, the Wannier function is more localized for Cu
than that for O, expalining the fact why the hopping among Cu is smaller than
that among O.

###  IV.2  Pb9Cu(PO4)6(OH)2

We construct a 2-band model $h_{2}(\mathbf{k})$ for the two bands near the
Fermi level in Pb9Cu(PO4)6(OH)2 for both Cu1 doping and Cu2 doping and in both
PM and FM phases. The 2-band model is constructed with $d_{xz}$ and $d_{yz}$
on Cu (at 1b for Cu1 doping and at 1c for Cu2 doping), and it only contains NN
hopping along all three directions in addition to the onsite energy term. The
form of the model is the same for both doping and for both PM and FM phases,
since they all preserve the spinless TR and $C_{3}$ symmetries for the 2 bands
near the Fermi level according to the Wannierization of DFT data. With
parameter values determined from the Wannierization of the DFT data, we can
see that the model gives very similar bands as the DFT band structure (Fig.
12). (Detailed expression and parameter values of the 2-band model can be
found in Appendix. B.2.)

In particular, for one specific Cu doping, we choose the same parameter values
for the NN hoppings for PM and FM phases, since the DFT values in the two
phases are very close (difference smaller than $0.1meV$); the only non-
negligible difference between PM and FM is just an onsite energy shift,
indicating that the effect of FM on the states near Fermi level is just
uniformly shifting their energy. The small bandwidth of the bands near the
Fermi level comes from the small hopping magnitudes (maximum hopping about
12meV), which comes from the small spread of the Wannier function of the basis
(the square root of the Wannier spread of each orbital is about $0.24a\approx
0.32c$). (See Appendix. B.2 for details.) We can see the Wannier spread of the
2-band model is larger than that of the Cu in the 4-band model, which is
consistent with the larger hopping here compared to the hopping among Cu in
the 4-band model.

In the DFT bands structure, the symmetry-protected gapless points at $\Gamma$
and $A$ are double Weyl points with chirality $\pm 2$. (See Appendix. B.2 for
details.) However, the band splitting along $\Gamma$-A is very small (maximum
splitting about $2$meV), and thus in our simplified NN-hopping model, we
neglect the band splitting along $\Gamma$-A for simplicity. Although such
simplification in our model will make the two double Weyl points merge into an
accidental nodal line along $\Gamma$-A, it will be convenient for future study
on the correlated physics, since it makes the eigenvectors of the Hamiltonian
independent of $k_{z}$.

Both Ref. [16] and Ref. [9] contain 2-band PM models constructed from $d_{xz}$
and $d_{yz}$ on Cu. However, the 2-band PM model in Ref. [16] has mirror
symmetry along $y$, which we do not include in our model since the DFT
calculation indicates a considerable breaking of the mirror symmetry. (See
Appendix. B.2 for details.) The 2-band PM models in both Ref. [16] and Ref.
[9] have considerably different parameter values than ours. Ref. [9] shows
double Weyl points at $\Gamma$ and A in their 2 band model, for which they
choose to include the small band splitting along $\Gamma$-A in their model.

Figure 12:  Comparison of the relaxed DFT (red) and 2-band tight-binding model
band structures (blue) and irreps for Pb9Cu(PO4)6(OH)2 in the specified Cu
doping and magnetic properties. The relaxed DFT and the NN-hopping tight-
binding model (blue) are in close agreement. The expression of the model and
the values of the model parameters are specified in Appendix. B.2. Note that
this is a 2-band model, so the quantum geometric tensor of both bands, taken
together, vanishes.

##  V  Interacting Hamiltonian

We use the constraint random phase approximation (cRPA) method[73, 74, 75, 76]
to compute the screened Coulomb interaction for the $(d_{xz},d_{yz})$ orbitals
of Cu near $E_{f}$ for two Cu-doped structure of Pb9Cu1(PO4)6O and
Pb9Cu1(PO4)6(OH)2. In Table. 3, we list the values of Hubbard-Kanamori
parameters, i.e., the onsite intra-orbital Hubbard $\mathcal{U}$, inter-
orbital $\mathcal{U}^{\prime}$, and onsite exchange $\mathcal{J}$. The
interacting Hamiltonian can be constructed as

$\displaystyle\hat{H}_{\text{int}}$
$\displaystyle=\mathcal{U}\sum_{im}n_{im\uparrow}n_{im\downarrow}+\mathcal{U}^{\prime}\sum_{i,m\neq
m^{\prime}}n_{im\uparrow}n_{im^{\prime}\downarrow}$ (5)
$\displaystyle+(\mathcal{U}^{\prime}-\mathcal{J})\sum_{i,m<m^{\prime},\sigma}n_{im\sigma}n_{im^{\prime}\sigma}$
$\displaystyle+\mathcal{J}\sum_{i,m\neq
m^{\prime}}(-c_{im\uparrow}^{\dagger}c_{im\downarrow}c_{im^{\prime}\downarrow}^{\dagger}c_{im^{\prime}\uparrow}+c_{im\uparrow}^{\dagger}c_{im\downarrow}^{\dagger}c_{im^{\prime}\downarrow}c_{im^{\prime}\uparrow}),$

where $i$ is the site index and $m$ the orbital index. We remark that the ab-
initio Hubbard-Kanamori parameters may need to be renormalized as the two
quasi-flat bands near $E_{f}$ also have weights of other orbitals.

Phase | $\mathcal{U}$ | $\mathcal{U}^{\prime}$ | $\mathcal{J}$
---|---|---|---
Pb9Cu1(PO4)6O Cu1 | 2.75 | 1.71 | 0.52
Pb9Cu1(PO4)6O Cu2 | 3.53 | 2.38 | 0.58
Pb9Cu1(PO4)6(OH)2 Cu1 | 2.88 | 1.99 | 0.45
Pb9Cu1(PO4)6(OH)2 Cu2 | 1.85 | 0.96 | 0.14
Table 3: The ab-initio Hubbard-Kanamori parameters. In the table,
$\mathcal{U}$, $\mathcal{U}^{\prime}$, $\mathcal{J}$ denotes the onsite intra-
orbital Hubbard, inter-orbital Hubbard, and onsite exchange interaction. All
numbers are in eV.

###  V.1 Comments on Flat Bands and Interactions

Flat bands are thought to be beneficial for strongly correlated phases since
the interaction strength necessarily dominates over the single-particle
bandwidth. However, the quantum geometry of the flat bands is another
essential ingredient in determining the nature of the resulting strongly
correlated phases.

For instance in the single-band square lattice Hubbard model at half filling,
the flat band limit $t\to 0$ yields a fully decoupled lattice (the atomic
limit) which is a perfect paramagnet. It is $O(t^{2}/U)$ corrections that
stabilize an anti-ferromagnetic phase. If, however, the interaction strength
is much larger than a set of isolated bands but smaller than the gap between
these bands and their complement, a different set of phases can emerge. In the
repulsive case, ferromagnetism can be proven [43, 44], and superconductivity
(or phase separation) for attractive case [43, 34, 26]. Features of the many-
body phase, for instance the mass of the Cooper pair, stiffness of the spin
wave, and a lower bound on the electron-phonon coupling are determined by
quantum geometry [34, 26, 77].

To entertain the possibility of superconductivity in flat bands, we recall
that the mean-field critical temperature in flat bands will be proportional to
the interaction strength [22]. Although this scenario is exponentially
improved over one-band BCS theory, it still requires an attractive interaction
of roughly 100 meV. This order of magnitude is larger than the bandwidth of
the two-band model in the OH structure, and roughly equal to the bandwidth of
the upper Cu bands in the four-band model for the O structure. (We note in
this case that the band gap is small, and projecting the interaction to the
flat bands may not be justified.) The repulsive Coulomb interaction we
computed is eV scale, and is much larger than the total bandwidth of both
tight-binding models.

Two recent papers [9, 20] have proposed models where one of the Cu bands is
perfectly flat. Although the two Cu bands together form an indecomposable
elementary band representation with trivial topology, analyzing only one of
the two bands (per spin) away from their degeneracy points at $\Gamma$ and $A$
can show strong quantum geometry. It is only appropriate to consider strong
coupling groundstate a single flat (gapless) band in the limit where the
interaction strength is much smaller than the bandwidth. This limit, while
interesting and worthy of study, places an upper bound on the interaction
strength which limits its applicability to a high-temperature phase.

##  VI Further Verification and Tests

###  VI.1 Immediate Experiments

As we experimentally found that samples are multi-phase, short-term
experiments should focus on isolating the different phases and characterizing
their compositions and structures. Also several more synthesis should be
performed to confirm that the outcomes are repeatable. Once we know reliably
the outcome of the synthesis the sample needs to be thoroughly characterized
with a wide range of methods. SEM/EDX will be a good first attempt to analyze
how many different phases are in the materials and what their respective
compositions are. It would be ideal if those phases can be separated either
manually of in different synthesis attempts targeting the respective
compositions found via SEM/EDX. As we showed it is possible to pick single
crystals from at least one phase.

###  VI.2 Chemical Structure Verification

Once phase-pure materials are obtained, they can be thoroughly characterized.
If single crystals are obtained, SXRD is ideal in combination with chemical
analysis such as SEM/EDX or, if enough samples can be separated, ICP-OES.
Otherwise PXRD in combination with Rietveld refinement will be necessary. If
the phases are not crystalline, the structural characterization becomes more
complicated and chemical analysis is a more important first step, likely to be
followed by high-resolution electron microscopy.

Characterization of physical properties is most meaningful on single-phase
materials and can be performed after through characterization of all
components of the samples. Similarly, theoretical analysis of the electronic
structures is most meaningful if the final crystal structures are known.

###  VI.3 Theoretical Analysis

Once the chemical structure is firmly established and the set of bands at the
Fermi level is settled, the bands must be fully analyzed based on the
principles outlined in this paper, including orbital, quantum geometric,
topological, and localization content in all the multiple phases that result
from the reaction. The interacting Hamiltonian must then be derived and with
it the values of the Hubbard $U$. Then the ground state of the system must be
determined. Even if not superconducting, one must still explain the levitation
properties shown in initial experiments - either large diamagnetism or some
flavor of ferromagnetism. The phonons, and more importantly the electron-
phonon interaction need to be obtained; an understanding is needed of the
rather flat Pb phonon bands on the $k_{z}=\pi$ plane that seem to become
negative at low temperature. If indeed, however, the sample does turn out to
be superconducting, short of being a fundamental discovery, it would also
point out the limitation of our theoretical understanding of the mechanisms
that create it.

##  VII Conclusions

Assuming the reported structure and the location of the Cu dopants, we have
performed first principle calculations of the electronic structure of
Pb9Cu1(PO4)6O and Pb9Cu1(PO4)6(OH)2. We find that the bands are almost flat,
localized on the Cu atoms, and have weak quantum geometry and trivial
topology. Hence in a ”flat-band superconductivity” scenario, already alluded
to previously, our calculation of the quantum geometric properties of the
active bands shows that such a phase would have extremely small superfluid
stiffness even at zero temperature. Instead, their flatness and the strong
Coulomb repulsion is compatible with the (anti/)ferromagnetic ground states
found in calculations. The phonon spectra seem to contain imaginary phonons
for both high and low-temperature paramagnetic phases (with the reported
unrelaxed structure). The experimental findings also do not, as of yet,
strongly suggest the presence of superconductivity. Hence either the reported
structure is incorrect, or the ab-initio calculations and the structure are
correct but then furthemore give rise to a non-superconducting ground-state,
or the ab-initio calculations are incorrect due to large correlations or other
factors. In fact, our preliminary experiments and structural solutions point
to LK-99 being a multi-phase material, where the part that structurally agrees
with a doped version of Pb apatite is transparent and thus probably not
superconducting. Either way, this points to the difficulty of predicting
and/or explaining superconductors by ab initio methods, even when they exist.
It also suggests that extra cases be taken in the literature, both
experimental and theoretical. Experimentally, samples need to be much more
carefully analyzed with a wide range of diffraction and spectroscopic methods.
The individual phases should be isolated and their properties should be
studied separately. Theoretically, one must at the very minimum check the
correct symmetry of the states, their localization and topology, and try to
obtain superconductivity from the first principle Hamiltonian, rather than
introducing it by hand in a BdG formalism. These represent hard challenges
that only serious investigations can overcome.

##  VIII Acknowledgements

H.H. and Y.J. were supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (Grant Agreement
No. 101020833). D.C. acknowledges the hospitality of the Donostia
International Physics Center, at which this work was carried out. J.H-A. is
supported by a Hertz Fellowship. D.C. and B.A.B. were supported by the Simons
Investigator Grant No. 404513, the Gordon and Betty Moore Foundation through
Grant No. GBMF8685 towards the Princeton theory program, the Gordon and Betty
Moore Foundation’s EPiQS Initiative (Grant No. GBMF11070), Office of Naval
Research (ONR Grant No. N00014-20-1-2303), BSF Israel US foundation No.
2018226, NSF-MERSEC (Grant No. MERSEC DMR 2011750). J.Y. is supported by the
Gordon and Betty Moore Foundation through Grant No. GBMF8685 towards the
Princeton theory program. B.A.B. and C.F. are also part of the SuperC
collaboration. Y.J. and S.B-C. acknowledge financial support from the MINECO
of Spain through the project PID2021- 122609NB-C21 and by MCIN and by the
European Union Next Generation EU/PRTR-C17.I1, as well as by IKUR Strategy
under the collaboration agreement between Ikerbasque Foundation and DIPC on
behalf of the Department of Education of the Basque Government. BAB also
acknowledges support the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (grant agreement no.
101020833). J.H. and D.C. are supported by DOE Grant No. DE-SC0016239. J.H is
also supported by a Hertz Fellowship. C.H.H, P.S.B, and E.L.G are supported by
DMR-1956403, and the Research Corporation for Science Advancement (Cottrell
Scholar Program) for non-tenured faculty. We are grateful to T. T. Debela, A.
Walsh, D. Scanlon, A. Rosen, and C. Musgrave for helpful discussions related
to defect formation enthalpy calculations. The authors acknowledge the use of
Princeton’s Imaging and Analysis Center, which is partially supported by the
Princeton Center for Complex Materials, a National Science Foundation (NSF) –
MRSEC program (DMR-2011750). The authors would also like to acknowledge G.
Cheng for helping with carbon coating. LMS is supported by the Gordon and
Betty Moore Foundation’s EPIQS initiative through Grant No. GBMF9064, as well
as the David and Lucille Packard foundation. SBL is supported by the National
Science Foundation Graduate Research Fellowship Program under Grant No.
(DGE-2039656). Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily
reflect the views of the National Science Foundation.

###### Contents

  1. I Introduction
  2. II Chemical Structure
  3. III ab initio Results
    1. III.1 Electronic Structure
    2. III.2 Phonon Spectrum
  4. IV Tight-Binding Models
    1. IV.1  Pb9Cu(PO4)6O
    2. IV.2  Pb9Cu(PO4)6(OH)2
  5. V  Interacting Hamiltonian
    1. V.1 Comments on Flat Bands and Interactions
  6. VI Further Verification and Tests
    1. VI.1 Immediate Experiments
    2. VI.2 Chemical Structure Verification
    3. VI.3 Theoretical Analysis
  7. VII Conclusions
  8. VIII Acknowledgements
  9. A Additional DFT results
    1. A.1 Crystal structure
    2. A.2 Orbital projections
    3. A.3 Bands of unrelaxed structures
    4. A.4 Phonon Spectra
    5. A.5 DFT computational details
  10. B Single-Particle Hamiltonians
    1. B.1  Pb9Cu(PO4)6O
    2. B.2  Pb9Cu(PO4)6(OH)2
  11. C Experimental Details
    1. C.1 Synthesis
    2. C.2 Refinement
    3. C.3 Powder Diffraction Comparisons

##  Appendix A Additional DFT results

###  A.1 Crystal structure

In Table. 4, we list the experimental and DFT relaxed lattice constants in
undoped and Cu-doped phases of Pb10(PO4)6O and Pb10(PO4)6(OH)2. In Table. 5,
we give the relaxed atomic positions.

Compound | Phase | $a$ | $c$ | Volume
---|---|---|---|---
Pb-O | Experiment[45] | 9.865 | 7.431 | 626.251
undoped | 10.018 | 7.485 | 650.585
Cu1-doped | 9.795 | 7.339 | 609.820
Cu2-doped | 9.868 | 7.410 | 624.970
Pb-(OH)2 | Experiment[46] | 9.866 | 7.426 | 625.991
undoped | 9.866 | 7.426 | 625.991
Cu1-doped | 9.731 | 7.296 | 598.352
Cu2-doped | 9.718 | 7.301 | 597.151
Table 4: The experimental and relaxed lattice constants in undoped and Cu-doped phases Pb10(PO4)6O and Pb10(PO4)6(OH)2 (written as Pb-O and Pb-(OH)2 for short). The experimental structures are given by Ref.[45, 46].  Compound | Phase | Atom | Position
---|---|---|---
Pb-O | undoped | tri-O | $(0,0,0.726)$
Cu1-doped | tri-O | $(0,0,0.821)$
Cu1 | $(\frac{1}{3},\frac{2}{3},0.498)$
Cu2-doped | tri-O | $(0,0,0.746)$
Cu2 | $(\frac{2}{3},\frac{1}{3},0.502)$
Pb-(OH)2 | undoped | tri-O | $(0,0,0.040),(0,0,0.540)$
tri-H | $(0,0,0.900),(0,0,0.400)$
Cu1-doped | tri-O | $(0,0,0.154),(0,0,0.762)$
tri-H | $(0,0,0.019),(0,0,0.628)$
Cu1 | $(\frac{1}{3},\frac{2}{3},0.473)$
Cu2-doped | tri-O | $(0,0,0.154),(0,0,0.630)$
tri-H | $(0,0,0.021),(0,0,0.497)$
Cu2 | $(\frac{2}{3},\frac{1}{3},0.498)$
Table 5: The relaxed atomic positions in undoped and Cu-doped phases
Pb10(PO4)6O and Pb10(PO4)6(OH)2 (written as Pb-O and Pb-(OH)2 for short).

###  A.2 Orbital projections

In Fig. 13 and Fig. 14, we show the orbital projections of Pb10(PO4)6O and
Pb10(PO4)6(OH)2, respectively.

Figure 13: Orbital projections of Pb9Cu1(PO4)6O in Cu1\- and Cu2-doping, PM
and FM phases, where the first row is the projection of $(d_{xz},d_{yz})$
orbitals of Cu, and the second row is the $(p_{x},p_{y})$ orbitals of
trigonal-O.  Figure 14: Orbital projections of Pb9Cu1(PO4)6(OH)2 in Cu1\- and
Cu2-doping, PM and FM phases, where the first row is the projection of
$(d_{xz},d_{yz})$ orbitals of Cu, and the second row is the $(p_{x},p_{y})$
orbitals of trigonal-O.

###  A.3 Bands of unrelaxed structures

As mentioned in the main text, the original experimental structures have
fractional occupancy for O or (OH)2. Hence the unrelaxed structures used here
have fixed the position for O and (OH)2 and remove the fractional occupancy.
This means that these unrelaxed structures are not experimental and need
further relaxation in order to obtain stable structures for DFT. The Cu-doped
phases also do not have experimental structures and need further relaxation.
For completeness, we also show in Fig. 15 the bands of unrelaxed structure for
Pb10(PO4)6O and Pb10(PO4)6(OH)2, in undoped, Cu1-doped, and Cu2-doped phases.

Figure 15: The DFT bands of unrelaxed structures. (a)-(c) are undoped,
Cu1-doped, and Cu2-doped bands for Pb10(PO4)6O, and (d)-(f) are for
Pb10(PO4)6(OH)2.

(a) Pb-O at low-T.

(b) Pb-O at high-T.

Figure 16: Phonon spectrum for relaxed Pb10(PO4)6O without doping Cu. The
negative branches get harden at $k_{3}=0$ plane as the temperature goes
higher. Instead, the imaginary phonon at $k_{3}=\pi$ are still soft, which may
be caused by the short cutoff in $c$-direction.

(a) Cu1-doped Pb9Cu1(PO4)6O in PM phase at low-T.

(b) Cu1-doped Pb9Cu1(PO4)6O in PM phase at high-T.

(c) Cu1-doped Pb9Cu1(PO4)6O in FM phase at low-T.

(d) Cu1-doped Pb9Cu1(PO4)6O in FM phase at high-T.

(e) Cu2-doped Pb9Cu1(PO4)6O in PM phase at low-T.

(f) Cu2-doped Pb9Cu1(PO4)6O in PM phase at high-T.

(g) Cu2-doped Pb9Cu1(PO4)6O in FM phase at low-T.

(h) Cu2-doped Pb9Cu1(PO4)6O in FM phase at high-T.

Figure 17: Phonon spectrum for relaxed Pb9Cu1(PO4)6O in Cu1\- and Cu2-doping
state. Different from the undoped structure, the doped ones present imaginary
phonon contributed by O and Cu at low-T, while their contributions are
eliminated at high-T.

###  A.4 Phonon Spectra

Here, we present the phonon spectrum of relaxed Pb-O, the doped Pb-O, and
Pb-(OH)2 compounds. The temperature effect is included via the smearing method
in harmonic approximation level. Without specific indication, these
calculations are performed within 111 cells. As shown below, all these phonon
spectrums present imaginary modes in the 111-cell results. Similar to the
phonon of undoped Pb-(OH)2 in Sec. III.2, the phonon presents a hardening
trend as the temperature goes higher. For the imaginary modes at $k_{3}=\pi$
plane, we expect they can be eased in larger supercell calculation as distant
force constants are included. In the doped structure, the O1 and Cu sites also
contribute to the imaginary branches at low temperatures, which are absent in
the high-temperature results. This is much more evident in the FM phonon
spectrum. Since O and Cu are much lighter than Pb, the current imaginary
phonons may not be physical, and further relaxation and larger supercell
calculations may eliminate them, which we leave for future study. We also
emphasize that current instabilities don’t determine the stability issues of
these compounds since larger supercell takes longer time.

(a) Cu1-doped Pb-(OH)2 in PM phase at low-T.

(b) Cu1-doped Pb-(OH)2 in PM phase at high-T.

(c) Cu1-doped Pb-(OH)2 in FM phase at low-T.

(d) Cu1-doped Pb-(OH)2 in FM phase at high-T.

(e) Cu2-doped Pb-(OH)2 in PM phase at low-T.

(f) Cu2-doped Pb-(OH)2 in PM phase at high-T.

Figure 18: Phonon spectrum for relaxed Pb-(OH)2 in Cu1\- and Cu2-doping state.
As in Cu-doped Pb-O compounds, O also contributes to the imaginary phonon
modes at low-T, which are absent at high-T.

###  A.5 DFT computational details

We use the Vienna ab-initio Simulation Package (VASP)[78, 79, 80, 81, 82] to
perform the ab-initio computations. The generalized gradient approximation
(GGA) with Perbew-Burke-Ernzerhof (PBE) exchange-correlation potential[83] is
adopted. A cutoff energy of 500 eV is adopted and no extra Hubbard U
correction is applied in the calculation of band structures. The maximally
localized Wannier functions (MLWFs) are constructed using Wannier90[84, 85,
86, 87] to obtain onsite energies and hoppings parameters of orbitals. The
phonon spectrum is calculated through density functional perturbation theory
(DFPT) assisted by Phonopy package [88].

Defect formation energy calculations were performed using a 2x2x2 supercell,
computed with the same base function, a 500 eV cutoff, and a 4 eV Hubbard U
correction for Cu. Bulk Cu, Cu2S, and Cu3P were used for the competing
chemical phases.

##  Appendix B Single-Particle Hamiltonians

In this Appendix, we give a detailed, symmetry-based construction of the
single-particle bands for two possibles structures of Cu doped lead apatite.
Both models show partially occupied topologically trivial bands dominated by
Cu $d$-orbitals at the Fermi level. We then fit our parameters to ab initio
Wannier calculations.

###  B.1  Pb9Cu(PO4)6O

We first provide the details of the single-particle Hamiltonians in the
PB9Cu(PO4)6O structure. The location of the Cu dopant, which is proposed to
replace Pb (although experimental characterization of the system is crucial to
test this hypothesis), can occur either at the 1b (Cu1 doping) or 1c position
(Cu2 doping). Note that Cu doping and structural relaxation reduce the space
group of the original Pb9Cu(PO4)6O compound to $P31^{\prime}$ generated by
translations, $C_{3}$, and spin-less time-reversal $\mathcal{T}$.

In both cases, DFT orbital projections indicate that the dominant orbitals in
the bands near the Fermi surface are the $d$-orbitals of the Cu atom at the
1b/1c position and the $p$ orbitals of the O atom at the 1a position (shown in
Fig. 19). We use lattice vectors
$\mathbf{a}_{1}=a(1,0,0),\mathbf{a}_{2}=C_{3}\mathbf{a}_{1},\mathbf{a}_{3}=c(0,0,0)$
where $C_{3}$ is a three-fold rotation. The Wannier centers/locations of Cu
atoms are $\mathbf{r}_{1b}=(\mathbf{a}_{1}+2\mathbf{a}_{2})/3+z\mathbf{a}_{3}$
and
$\mathbf{r}_{1c}=(2\mathbf{a}_{1}+\mathbf{a}_{2})/3+z^{\prime}\mathbf{a}_{3}$.
The realistic structures show a small vertical displacement between O and Cu
atoms (which are not fixed by symmetry since $C_{3}$ is in-plane) given by
$z,z^{\prime}\sim.25$.

(a) (b)

Figure 19: Real Space lattice and hoppings. (a) Hopping structure of the
4-band model in 2D. The only relevant 3D term is an interlayer $t^{\prime}$
hopping between Cu orbitals. (b) Hopping structure of the 2-band model in 2D.

The Cu $d$ orbitals and O $p$ orbitals both transform in the 2D irrep
${}^{1}E{}^{2}E$ with angular momentum $\pm 1\,\mathrm{mod}\,3$. Topological
quantum chemistry [47] provides the following momentum space irreps when these
local orbitals are induced to the space group:

$\displaystyle{}^{1}E{}^{2}E_{1a}\uparrow P31^{\prime}$
$\displaystyle=\Gamma_{2}\Gamma_{3}+A_{2}A_{3}+K_{2}+K_{3}+H_{2}+H_{3},$ (6)
$\displaystyle{}^{1}E{}^{2}E_{1b}\uparrow P31^{\prime}$
$\displaystyle=\Gamma_{2}\Gamma_{3}+A_{2}A_{3}+K_{1}+K_{3}+H_{1}+H_{3}$
$\displaystyle{}^{1}E{}^{2}E_{1c}\uparrow P31^{\prime}$
$\displaystyle=\Gamma_{2}\Gamma_{3}+A_{2}A_{3}+K_{1}+K_{2}+H_{1}+H_{2}$

and the irrep notation is contained in the character tables below.

$\displaystyle\begin{array}[]{c|ccc}31^{\prime}&1&C_{3}\\\ \hline\cr A&1&1\\\
{}^{1}E{}^{2}E&2&-1\\\
\end{array}\qquad\begin{array}[]{c|ccc}31^{\prime}&1&C_{3}\\\
\hline\cr\Gamma_{1}&1&1\\\ \Gamma_{2}\Gamma_{3}&2&-1\\\
\end{array}\qquad\begin{array}[]{c|ccc}31^{\prime}&1&C_{3}\\\ \hline\cr
A_{1}&1&1\\\ A_{2}A_{3}&2&-1\\\
\end{array}\qquad\begin{array}[]{c|ccc}31^{\prime}&1&C_{3}\\\ \hline\cr
K_{1}&1&1\\\ K_{2}&1&e^{\frac{2\pi i}{3}}\\\ K_{3}&1&e^{-\frac{2\pi i}{3}}\\\
\end{array}\qquad\begin{array}[]{c|ccc}31^{\prime}&1&C_{3}\\\ \hline\cr
H_{1}&1&1\\\ H_{2}&1&e^{\frac{2\pi i}{3}}\\\ H_{3}&1&e^{-\frac{2\pi i}{3}}\\\
\end{array}$ (7)

Note that the $K^{\prime}$ and $H^{\prime}$ points can be obtained by time-
reversal, and the only irreps at $M$ and $L$ are the trivial irrep of
$\mathcal{T}$. More details can be found on the Bilbao Crystallographic Server
https://www.cryst.ehu.es/. We see that degeneracies occur at the $\Gamma$
point and $A$, and the rest of the bands are split.

We now build a short ranged tight-binding model from these orbitals. We will
primarily discuss the para-magnetic (spin-unpolarized) DFT results which
provide a single-particle band structure. Including Coulomb repulsion favors a
ferromagnetic groundstate which is expected due to the flat bands obtained in
the tight-binding models.

For both dopant positions, we write the Hamiltonian as

$\displaystyle
h_{4}(\mathbf{k})=\begin{pmatrix}h_{C}(\mathbf{k})&h_{CO}(\mathbf{k})\\\
h^{\dagger}_{CO}(\mathbf{k})&h_{O}(\mathbf{k})\end{pmatrix}$ (8)

describing the two Wannier orbitals dominated by Cu and O respectively, and
their coupling. The symmetries of the model have representations

$\displaystyle D[C_{3}]h(\mathbf{k})D^{\dagger}[C_{3}]$
$\displaystyle=h(C_{3}\mathbf{k}),\qquad
h^{*}(\mathbf{k})=h(-\mathbf{k}),\qquad
h(\mathbf{k}+\mathbf{G})=V[\mathbf{G}]h(\mathbf{k})V^{\dagger}[\mathbf{G}]$
(9)

where $D[C_{3}]=R_{2\pi/3}\oplus R_{2\pi/3}$,
$V[\mathbf{G}]=e^{-i\mathbf{G}\cdot\mathbf{r}_{1b}}\sigma_{0}\oplus\sigma_{0}$
and $R_{2\pi/3}=e^{-i\frac{2\pi}{3}\sigma_{2}}$ is the 2D rotation matrix.

The O block is dominated by a relative potential $\mu_{O}$, an in-plane
$\mathbf{a}_{1}$ hopping matrix $t_{O}$, and an out-of-plane $\mathbf{a}_{3}$
hopping matrix $t^{\prime}_{O}$. For higher accuracy, we also include
$\mathbf{a}_{1}\pm\mathbf{a}_{3}$ hoppings denoted $t^{\pm}_{O}$. The O-O
Hamiltonian is

$\displaystyle
h_{O}(\mathbf{k})=\mu_{O}\sigma_{0}+\left(t^{\prime}_{O}e^{-i\mathbf{k}\cdot\mathbf{a}_{3}}+\sum_{n=0}^{2}R_{2\pi/3}^{n}(t_{O}e^{-i\mathbf{k}\cdot
C_{3}^{n}\mathbf{a}_{1}}+t^{+}_{O}e^{-i\mathbf{k}\cdot
C_{3}^{n}(\mathbf{a}_{1}+\mathbf{a}_{3})}+t^{-}_{O}e^{-i\mathbf{k}\cdot
C_{3}^{n}(\mathbf{a}_{1}-\mathbf{a}_{3})})R^{n}_{-2\pi/3}+h.c.\right)\ .$ (10)

The Cu $d_{xz},d_{yz}$ orbitals have the same $C_{3}$ representations, and
their block has the analogous form (although we shall see that the hopping
strengths are much weaker):

$\displaystyle
h_{C}(\mathbf{k})=\mu_{Cu}\sigma_{0}+\left(t^{\prime}_{Cu}e^{-i\mathbf{k}\cdot\mathbf{a}_{3}}+\sum_{n=0}^{2}R_{2\pi/3}^{n}(t_{Cu}e^{-i\mathbf{k}\cdot
C_{3}^{n}\mathbf{a}_{1}}+t^{+}_{Cu}e^{-i\mathbf{k}\cdot
C_{3}^{n}(\mathbf{a}_{1}+\mathbf{a}_{3})}+t^{-}_{Cu}e^{-i\mathbf{k}\cdot
C_{3}^{n}(\mathbf{a}_{1}-\mathbf{a}_{3})})R^{n}_{-2\pi/3}+h.c.\right)\ .$ (11)

Lastly, the hybridization term contains the nearest-neighbor Cu-O hopping $T$,
as well as out-of-plane couplings $T_{\pm}$:

$\displaystyle
h_{CO}(\mathbf{k})=\sum_{n=0}^{2}R^{n}_{2\pi/3}(Te^{-i\mathbf{k}\cdot
C_{3}^{n}\boldsymbol{\delta}}+T_{+}e^{-i\mathbf{k}\cdot
C_{3}^{n}(\boldsymbol{\delta}+\mathbf{a}_{3})}+T_{-}e^{-i\mathbf{k}\cdot
C_{3}^{n}(\boldsymbol{\delta}-\mathbf{a}_{3})})R^{n}_{-2\pi/3}\ .$ (12)

Here $\boldsymbol{\delta}=-\mathbf{r}_{1}b/\mathbf{r}_{1}c$ for the Cu1 doping
model and Cu2 doping model respectively.

Incorporating the dominant terms from DFT, we use the following hopping
matrices for Cu1 doping:

$\displaystyle\mu_{O}$ $\displaystyle=3.803,\quad
t_{O}=-0.018(\sigma_{1}-i\sigma_{2})/2,\quad
t^{\prime}_{O}=-0.074\sigma_{0},\quad\mu_{C}=4.094,\quad
t_{C}=-0.0035(\sigma_{1}-i\sigma_{2})/2$ (13) $\displaystyle T$
$\displaystyle=\begin{pmatrix}0.056&0.0190\\\
-0.056&-0.045\end{pmatrix},\qquad\text{all other terms}=0$

and Cu2 doping:

$\displaystyle\mu_{O}$ $\displaystyle=3.624,\quad
t_{O}=-0.025(\sigma_{0}+\sigma_{3})/2,\quad
t^{\prime}_{O}=-0.021\sigma_{0},\quad
t_{O}^{+}=.007(\sigma_{1}-i\sigma_{2})/2,\quad
t_{O}^{-}=.003(\sigma_{1}-i\sigma_{2})/2$ (14) $\displaystyle\mu_{C}$
$\displaystyle=3.736,\quad t_{C}=R_{-2\pi/3}\begin{pmatrix}0&0\\\
.006&-.002\end{pmatrix}R_{2\pi/3},\quad t_{C}^{\prime}=-.004\sigma_{0},\quad
t_{C}^{+}=0,\quad t_{C}^{-}=-.003(\sigma_{0}+\sigma_{3})/2$ $\displaystyle T$
$\displaystyle=\left(\begin{array}[]{cc}-0.005&0.021\\\ -0.012&0.004\\\
\end{array}\right),\quad T_{+}=-.004(\sigma_{1}+i\sigma_{2})/2,\quad
T_{-}=-.006(\sigma_{1}+i\sigma_{2})/2\ .$

The single-particle second-quantized Hamiltonian is

$\displaystyle H_{4}$
$\displaystyle=\sum_{\mathbf{k},\alpha\beta,\sigma=\uparrow,\downarrow}c^{\dagger}_{\mathbf{k}\alpha,\sigma}[h_{4}(\mathbf{k})]_{\alpha\beta}c^{\dagger}_{\mathbf{k}\beta,\sigma}$
(15)

using $SU(2)$ spin symmetry due to the small spin-orbit coupling. The Fermi
level of the single-particle model occurs at $3/4$ filling of the upper two
bands in each spin sector, leading to a total of one hole per unit cell in the
spinful model.

In addition to an accurate modeling of the spectrum, it is crucial to study
the localization and topological properties of the bands, whose effect on the
many-body physics is magnified by a large density of states. To do so, we
compute the two-band non-abelian Wilson loop in Fig. 20 and the two-band
Fubini-Study metric at $k_{z}=0,\pi$ in Fig. 21. Our results are consistent
with tight localization, as expected from the atomic band representation
formed by the Cu $d$-orbitals.

We emphasize that, although the two Cu bands together are trivial and form an
EBR, the double degeneracy enforced by $C_{3}$ and $\mathcal{T}$ (the
$\Gamma_{2}\Gamma)3$ and $A_{2}A_{3}$ irreps) creates non-trivial quantum
geometry within a _single_ band. In particular, the Fubini-Study metric will
diverge at the touching, as required by the Berry curvature monopole (pointed
out in Ref. [9]) or double Weyl point.

Figure 20: Wilson loops. We compute the non-abelian Wilson loop
$W(k_{2},k_{3})=\exp i\oint dk_{1}A_{1}(\mathbf{k})$ for the Cu1 phase (a) and
Cu2 phase (c). The Wilson loop phases/Wannier centers $\vartheta(k_{2},k_{3})$
show no winding, but show greater variation in the $k_{z}=\pi$ plane. Direct
comparison of the Wilson loop spectra at $k_{3}=0$ (blue) and $k_{3}=\pi)$
(red) for Cu1 (b) and Cu2 (d). Stronger hybridization with the O bands at
$k_{3}=\pi$ is responsible for the dispersion of the Wilson loop.

Figure 21: The 2D Fubini-Study metric
$g(\mathbf{k})=\frac{1}{2}\mathop{\mathrm{Tr}}(\partial_{x}P\partial_{x}P+\partial_{y}P\partial_{y}P)$
(in unit of $a^{2}$ of the upper two bands for the DFT calculation and the
tight-binding model Eq. (2).

Fig. 11 shows that Cu and O DFT bands nearly touch at the $H$ point,
indicating close proximity to a band inversion. We check that tuning the
$t^{\prime}$ vertical hopping of the O orbitals realixes this band crossing,
resulting in a topological phase transition into a semi-metal (see Fig. 22).

Figure 22: Inverted Semimetal phase obtained for the Cu1 structure by
increasing the O out-of-plane hopping $t^{\prime}_{O}$ in $h_{4}(\mathbf{k})$.
We compare the dispersions of the $C_{3}$ irreps along the $KH$ line in the
normal (a) and inverted (b) phases. The inverted phase is a topological
semimetal due to the protected crossing of the different $C_{3}$ eigenvalues.
Note $C_{3}$ is preserved along all points in $KH$. (c) show the band
structure in the inverted phase.

DFT shows that ferromagnetism (and also a degenerate anti-ferromagnetism) is
favored over the paramagnetic (spin-degenerate) single-particle bands given
above. We find that the spin-polarized bands at the Fermi surface have
essentially identical parameters as Eq. (13), whereas small adjustments in the
parameters are required for Cu2 doping:

$\displaystyle\mu_{O}$ $\displaystyle=3.695,\quad
t_{O}=-0.025(\sigma_{0}+\sigma_{3})/2,\quad
t^{\prime}_{O}=-0.022\sigma_{0},\quad
t_{O}^{+}=.006(\sigma_{1}-i\sigma_{2})/2,\quad
t_{O}^{-}=.003(\sigma_{1}-i\sigma_{2})/2$ (16) $\displaystyle\mu_{C}$
$\displaystyle=3.751,\quad t_{C}=R_{-2\pi/3}\begin{pmatrix}0&0\\\
.006&-.002\end{pmatrix}R_{2\pi/3},\quad t_{C}^{\prime}=-.004\sigma_{0},\quad
t_{C}^{+}=0,\quad t_{C}^{-}=-.003(\sigma_{0}+\sigma_{3})/2$ $\displaystyle T$
$\displaystyle=\left(\begin{array}[]{cc}-0.005&0.021\\\ -0.013&0.005\\\
\end{array}\right),\quad T_{+}=0,\quad T_{-}=-.006(\sigma_{1}+i\sigma_{2})/2\
.$

Their main effect is to further shrink the gap between the O and Cu bands, as
can be seen in Fig. 23

Figure 23: The ferromagnetic bands are very similar to the paramagnetic bands
(see Main Text). The essential difference between the two cases is the Fermi
level, which fills $3/4$ of the spin-degenerate Cu bands in the paramagnetic
case, but fills $1/2$ of the spin-polarized Cu bands in the ferromagnetic
case.

###  B.2  Pb9Cu(PO4)6(OH)2

We provide the details of the single-particle Hamiltonians in the
Pb9Cu(PO4)6(OH)2 structure for both Cu1 and Cu2 doping and in both PM and FM
phases. For this structure, we have space group P3. In all the cases, we only
have two bands near the Fermi level Fig. 9, they are trivial atomic bands
given by $d_{xz}$ and $d_{yz}$ orbitals of Cu. (Cu atoms are at 1b for Cu1 and
at 1c for Cu2.) Wanneirization of the DFT bands would directly generate a DFT
2-band tight-binding model, which can generate band structure that perfectly
matches the DFT ones. The DFT 2-band tight-binding model shows that the
hoppings longer than NN are small (_i.e._ , less than 1meV), and thus we will
build a 2-band NN-hopping model for the two bands in each case. Briefly, the
form of the model is the same for all cases due to the same symmetries. The
hopping parameter values in the model are the same for the PM and FM phases
for one specific Cu doping, though the onsite energies are different. The
hopping values do change across different Cu doping, as indicated by the
different bands in Fig. 9 for different Cu doping. Details are discussed
below.

The 2-band NN hopping model that we build with $d_{xz}$ and $d_{yz}$ at the
same position in each unit cell reads

$H_{2band}=\sum_{\bm{k}}c^{\dagger}_{\bm{k},Cu}\left[E_{Cu}+2t_{z}\cos(k_{z}c)+\left(\sum_{n=0,1,2}R_{2\pi/3}^{n}t_{\bm{a}_{1}}R_{-2\pi/3}^{n}e^{-\mathrm{i}(C_{3}^{n}\bm{a}_{1})\cdot\bm{k}}+h.c.\right)\right]c_{\bm{k},Cu}\
,$ (17)

where e

$t_{\bm{a}_{1}}=\left(\begin{matrix}t_{xx}&t_{xy}\\\
t_{yx}&t_{yy}\end{matrix}\right)$ (18)

is the hopping matrix along $\bm{a}_{1}$,
$c^{\dagger}_{\bm{k},Cu}=(c^{\dagger}_{\bm{k},Cu,d_{xz}},c^{\dagger}_{\bm{k},Cu,d_{yz}})$.

For Cu1 doping, the parameter values are

$E_{Cu}=-0.0326\text{eV}\ ,\ t_{z}=-0.0062\text{eV}\ ,\
t_{xx}=-0.0094\text{eV}\ ,\ t_{xy}=0.0066\text{eV}\ ,\
t_{yx}=-0.0115\text{eV},t_{yy}=0\ .$ (19)

The in-plane hopping is shown in Fig. 19(b).Compared to the most general
symmetry-allowed NN hoppings, we directly neglect the off-diagonal part of the
$\bm{a}_{3}$ hopping (as it is about 0.5meV, equivalent to maximum band
splitting about $2$meV along $\Gamma$-A). We also set $t_{yy}=0$ since it is
smaller than 2meV. The band structure from the model Eq. (17) is shown in Fig.
12(a), which has a good agreement with the DFT band structure. The flatness of
the bands, which form an indecomposablt atomic representation of the Cu
orbitals, is merely due to the small hopping amplitudes/localized Wannier
functions. To be specific, the square root of the Wannier spread of each
orbital is $2.402\text{\AA}=0.2469a=0.3293c$, indicating that the Wannier
function is very localized. Owing to $t_{yx}/t_{xy}=-1.74\neq-1$, the mirror
symmetry along $y$ is broken.

Since this is a 2-band model, the geometric properties (_e.g._ , Fubini-Study
metric) of the two bands, taken together, vanishes. However, we may look at
the Wilson loop of the one of the two bands, for which we choose lower band.
The Wilson loop is ill-defiend for $k_{z}=0$ and $k_{z}=\pi$ planes, since the
two bands touch on the two planes at $\Gamma$ or A (enforced by symmetries),
and thus we look at $k_{z}=\pm\pi/2$ for the DFT 2-band tight-binding model,
which is shown in Fig. 24(a). Clearly, the lower band of the DFT 2-band model
has nonzero and opposite Chern numbers at $k_{z}=\pm\pi/2$, showing that the
gapless points at $\Gamma$ and A are double Weyl points with chirality $\pm
2$. However, in our simplified model (Eq. (17)), we neglect the band splitting
along $\Gamma$-A, since it is very small. Although such simplfiication merge
the two double Weyl points into an accidental nodal line along $\Gamma$-A, it
would be convenience for later study of correlated phases based on our model,
since it make the eigenvectors of the Hamiltonian independent of $k_{z}$.

Figure 24: The Wilson loop spetrum of the lower band of the DFT 2-band model
in each specified case on $k_{z}=\pi/2$. The Wilson loop $W$ is calcuated
along $\bm{b}_{2}$, and is ploted as a function of $k_{1}$ is along
$\bm{b}_{1}$.

Now we discuss the ferromagnetic case for Cu1 doping. With ferromagnetism,
only the spin-down bands remain near the Fermi level Fig. 9. The DFT spin-down
bands are very close to the DFT paramagnetic bands as shown in Fig. 12(a,b).
The effect of FM is approximately just a shift of the onsite energy for bands
near Fermi level, as the hopping values for the DFT FM spin-down two bands are
extremely similar to those for the DFT PM two bands (difference smaller
$0.1meV$) except of a considerable shift of the onsite energy. Therefore, we
can build a 2-band NN-hopping model that has the same form as Eq. (17), and
choose the same hopping parameter values as Eq. (19), while shifting onsite
energy to

$E_{Cu}=0.008907\text{eV}\ .$ (20)

Again, the band structure from the model has a good agreement with the DFT
band structure, as shown in Fig. 12. Again, the DFT 2-band tight-binding model
has double Weyl points with chirality $\pm 2$ at $\Gamma$ and A (Fig. 24(b)),
which we choose to merge into an accidental nodal line along $\Gamma$-A in our
simplified model.

We now discuss the Cu2 doping parameters, which are very similar to the Cu1
doping. A key difference is that the Cu is now at 1c, but this is simply a
choice of origin in the two-band model we build and does not affect the
parameterization of the Hamiltonian. In the PM phase, the hopping parameter
values read

$\displaystyle\begin{split}&E_{Cu}=-0.0207\text{eV},\ t_{z}=-0.0065\text{eV},\
t_{xx}=\frac{1}{4}\left(3t_{2}-\sqrt{3}t_{1}\right),\
t_{xy}=\frac{1}{4}\left(-3t_{1}-\sqrt{3}t_{2}\right)\\\
&t_{yx}=\frac{1}{4}\left(t_{1}-\sqrt{3}t_{2}\right),t_{yy}=\frac{1}{4}\left(\sqrt{3}t_{1}+t_{2}\right)\end{split}$
(21)

with

$t_{1}=0.0112\text{eV}\ ,\ t_{2}=-0.0027\text{eV}\ .$ (22)

The reason for us to parameterize the in-plane hoppings by only $t_{1}$ and
$t_{2}$, because along $\bm{a}_{2}$ direction the hopping matrix approximately
reads

$R_{2\pi/3}t_{\bm{a}_{1}}R_{-2\pi/3}=t_{\bm{a}_{2}}=\left(\begin{matrix}0&0\\\
t_{1}&t_{2}\end{matrix}\right)$ (23)

according to the Wannierization of the DFT data, where the neglected elements
are smaller than 1meV. In particular, $t_{xy}=-0.0072\text{eV}$ and
$t_{yx}=0.0040\text{eV}$ indicate the breaking of the mirror as $y$ since
$t_{xy}/t_{yx}=-1.81\neq-1$. The band structure from the model Eq. (17) again
has a good agreement with the DFT band structure, as shown in Fig. 12(b).
Small hopping again the small spread of the Wannier function of the basis: the
square root of the Wannier spread of each orbital is
$2.271\text{\AA}=0.2334a=0.3113c$. Again, the DFT 2-band tight-binding model
has double Weyl points with chirality $\pm 2$ at $\Gamma$ and A (Fig. 24(c)),
which we choose to merge into an accidental nodal line along $\Gamma$-A in our
simplified model.

In the ferromagnetic case, the model parameters again only defer from Eq.
(LABEL:eq:paramagnetic_2band_parameters_Cu_2) by an onsite energy shift:

$E_{Cu}=0.0023\text{eV}\ .$ (24)

Again, the band structure from the model has a good agreement with the DFT
band structure, as shown in Fig. 12(d). Again, the DFT 2-band tight-binding
model has double Weyl points with chirality $\pm 2$ at $\Gamma$ and A (Fig.
24(d)), which we choose to merge into an accidental nodal line along
$\Gamma$-A in our simplified model.

##  Appendix C Experimental Details

In this Appendix, we give detailed information on synthesis and
characterization of lead apatites.

###  C.1 Synthesis

Well homogenized powders of Pb2(SO4)O were synthesized by grinding
stoichiometric amounts of PbO (Sigma Aldrich, $>$99.0 percent) and PbSO4
(Sigma Aldrich, $>$99.0 percent) in an agate mortar and pestle. The mixture
was loaded in alumina crucibles in 3g amounts and then placed in a quartz tube
and sealed under dynamic vacuum. Experiments found that sealing with oxygen
present (i.e. not backfilling with Ar) resulted in the highest purity of
product. Sample purity was confirmed via powder X-ray diffraction using a STOE
Stadi P powder X-ray diffractometer equipped with a Mo K$\alpha$
($\lambda$$=0.71073$ Å) sealed-tube X-ray source and graphite monochromator at
room temperature in transmission geometry Fig. 25

Figure 25: Comparisons of experimental (black) and calculated powder X-ray
diffraction patters of Pb2(SO4)O (red).

Cu3P was grown phase pure by mixing elemental Cu (Alfa Aesar, -170+270 mesh,
99.9 percent) and elemental red Phosphorus (Alfa Aesar, 99.999+ percent) in a
3:1 stoichiometric ratio. The sample was placed in an alumina crucible and
subsequently sealed in a quartz tube. For this synthetic step the quartz tube
was evacuated three times and backfilled with argon to prevent any oxidation
before sealing. The tube was then loaded in a furnace heated first to 450 ∘C
over 2 hours and held there for 3 hours to prereact the phosphorus. The
furnace was then ramped to 1050 ∘C over 4 hours, kept at this temperature for
8 hours, and then shut off to cool quickly. This saved nearly 24 hours in
synthesis time compared to the original synthesis and results in a phase pure
product of the target compound, collected in transmission geometry on the same
instrument Fig. 26.

Figure 26: Comparisons of experimental (black) and calculated powder X-ray
diffraction patters of Cu3P (red).

Finally, a 1:1 molar ratio of Cu3P and Pb2(SO4)O were ground into a fine
powder in mortar and pestle. The powders were once again loaded into an
alumina crucible, placed in a quartz tube, and sealed under dynamic vacuum.
The sample was then heated to 950 ∘C, kept at this temperature for 6hrs, and
then shut off to cool quickly.

Figure 27: Comparisons of previously reported (black) and collected powder
X-ray data for proposed Pb10-xCux(PO4)6O (red). Figure 28: A representative
sample of the final product phase, ground into a powder. We observe at least 3
phases with seemingly different properties. Figure 29: EDS of the red and
white transparent phase showing pockets of Cu embedded in a phase consistent
with Pb10-xCux(PO4)6(OH)(SH). Figure 30: EDS of the metallic gray phase
showing Cu2S embedded in a phase consistent with Pb10-xCux(PO4)6(OH)(SH).

###  C.2 Refinement

Calculated precession images for data in refined herein are shown below.

Figure 31: Precession images for a thin transparent crystal. The only
systematic absences are [00l]=2n agreeing with the space group P63/$m$.

Collection and integration statistics for the full hemisphere data set are
shown below.

Table S 6: Collection and integration statistics for lead appetite structures. Crystal Dimension (mm) | 0.55 $\times$ 0.092 $\times$ 0.176 mm |  |  |
---|---|---|---|---
Temperature (K) | 298(2) |  |  |
Radiation Source, $\lambda$ (Å) | Mo K$\alpha$, $\lambda$=0.71073 Å |  |  |
Absorption Correction | analytical |  |  |
Space Group | P63/$m$ |  |  |
$a$ (Å) | 9.8508(1) |  |  |
$c$ (Å) | 7.4395(2) |  |  |
Cell Volume (Å3) | 625.198(19) |  |  |
Absorption Coefficient (mm-1) | 67.952 |  |  |
$\theta_{min}$ , $\theta_{max}$ | 2.39, 36.34 |  |  |
Number of Reflections | 32013 |  |  |
Unique Reflections (I >3$\sigma$, all) | 763, 1066 |  |  |
Rint(I>3$\sigma$, all) | 6.54, 7.40 |  |  |

Freely refining occupancies for both an O at [0 0 $1/4$] and S located at [0 0
0] results in site occupancy factors of 0.54 and 0.56, respectively. The
freely refined occupancies having a summation very close to 1 seems suggestive
that this could be a physically real interpretation of the structural
solution. After refining the thermal parameters anisotropically, a final
refinement converges with a goodness of fit (GoF) parameter of 1.22 and R
value of 3.94 compared to all reflection. Here, our site occupancies for O
(1.08) and S (1.12) sum to above 1. If we decide to restrict the occupancies
of the 2$b$ S and 2$a$ O to sum to a site occupancy factor of 1, we end
refinement with a final composition Pb10(PO4)6(OH)0.94(SH)1.06. This
refinement has nominally the same refinement statistics (GoF(all) $=1.22$,
R(all) $=3.95$).

Table S 7: Refinement Statistics for lead appetite structures without Cu doping. Refined Composition | Pb10(PO4)6(OH)0.94(SH)1.06 |  |  | Pb10(PO4)6(OH) 1.12 (SH) 1.08
---|---|---|---|---
Refinement Method | F2 |  |  | F2
Number of Parameters | 43 |  |  | 44
R(I>$3\sigma$), Rw(I>$3\sigma$) | 2.48,5.43 |  |  | 2.47, 5.43
R(all), Rw(all) | 3.95, 5.70 |  |  | 3.94, 5.69
S(I>$3\sigma$), S(all) | 1.39, 1.22 |  |  | 1.39, 1.22
$\Delta\rho_{max}$ , $\Delta\rho_{min}$ (e Å-3) | 3.38, -3.19 |  |  | 3.35, -3.21
Table S 8: Refined atomic coordinates for Pb10(PO4)6(OH)0.94(SH)1.06. Site |  | Wyckoff Position |  | x |  | y |  | z |  | Occupancy |  |
---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 6h |  | 0.24486(4) |  | 0.24667(4) |  | $1/4$ |  | 1 |  |
Pb2 |  | 4f |  | $2/3$ |  | $1/3$ |  | -0.00339(4) |  | 1 |  |
P1 |  | 6h |  | 0.3729(2) |  | 0.4020(2) |  | $3/4$ |  | 1 |  |
O1 |  | 6h |  | 0.4851(7) |  | 0.3366(8) |  | $3/4$ |  | 1 |  |
O2 |  | 12i |  | 0.2650(6) |  | 0.3480(6) |  | 0.5834(6) |  | 1 |  |
O3 |  | 6h |  | 0.5337(8) |  | 0.4166(7) |  | $1/4$ |  | 1 |  |
S1 |  | 2b |  | 0 |  | 0 |  | 0 |  | 0.53 |  |
O4 |  | 2a |  | 0 |  | 0 |  | $1/4$ |  | 0.47 |  |
Table S 9: Refined anisotropic displacement parameters for Pb10(PO4)6(OH)0.94(SH)1.06. Site |  | U11 |  | U22 |  | U33 |  | U12 |  | U13 |  | U23 |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 0.02052(15) |  | 0.01996(15) |  | 0.02682(15) |  | 0.01447(12) |  | 0.00000 |  | 0.00000 |  |
Pb2 |  | 0.01890(12) |  | 0.01890(12) |  | 0.01231(13) |  | 0.00945(6) |  | 0.00000 |  | 0.00000 |  |
P1 |  | 0.0088(7) |  | 0.0096(7) |  | 0.0125(7) |  | 0.0055(6) |  | 0.00000 |  | 0.00000 |  |
O1 |  | 0.022(3) |  | 0.028(3) |  | 0.024(3) |  | 0.020(3) |  | 0.00000 |  | 0.00000 |  |
O2 |  | 0.022(2) |  | 0.040(3) |  | 0.0212(19) |  | 0.018(2) |  | -0.0064(17) |  | -0.0101(19) |  |
O3 |  | 0.029(3) |  | 0.010(2) |  | 0.034(3) |  | 0.006(2) |  | 0.00000 |  | 0.00000 |  |
S1 |  | 0.015(2) |  | 0.015(2) |  | 0.044(5) |  | 0.0077(11) |  | 0.00000 |  | 0.00000 |  |
O4 |  | 0.014(6) |  | 0.014(6) |  | 0.054(16) |  | 0.007(3) |  | 0.00000 |  | 0.00000 |  |
Table S 10: Refined atomic coordinates for Pb10(PO4)6(OH)1.12(SH)${}_{1.08}.$ Site |  | Wyckoff Position |  | x |  | y |  | z |  | Occupancy |  |
---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 6h |  | 0.24486(4) |  | 0.24667(4) |  | $1/4$ |  | 1 |  |
Pb2 |  | 4f |  | $2/3$ |  | $1/3$ |  | -0.00339(4) |  | 1 |  |
P1 |  | 6h |  | 0.3729(2) |  | 0.4020(2) |  | $3/4$ |  | 1 |  |
O1 |  | 6h |  | 0.4851(7) |  | 0.3366(8) |  | $3/4$ |  | 1 |  |
O2 |  | 12i |  | 0.2650(6) |  | 0.3480(6) |  | 0.5834(6) |  | 1 |  |
O3 |  | 6h |  | 0.5337(8) |  | 0.4166(7) |  | $1/4$ |  | 1 |  |
S1 |  | 2b |  | 0 |  | 0 |  | 0 |  | 0.54 |  |
O4 |  | 2a |  | 0 |  | 0 |  | $1/4$ |  | 0.56 |  |
Table S 11: Refined anisotropic displacement parameters for Pb10(PO4)6(OH)1.12(SH)${}_{1.08}.$ Site |  | U11 |  | U22 |  | U33 |  | U12 |  | U13 |  | U23 |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 0.02052(15) |  | 0.01996(8) |  | 0.00796(9) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
Pb2 |  | 0.00687(11) |  | 0.00721(10) |  | 0.00817(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
P1 |  | 0.00854(10) |  | 0.00664(10) |  | 0.00948(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
O1 |  | 0.00854(10) |  | 0.00664(10) |  | 0.00948(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
O2 |  | 0.00854(10) |  | 0.00664(10) |  | 0.00948(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
O3 |  | 0.00854(10) |  | 0.00664(10) |  | 0.00948(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
S1 |  | 0.00854(10) |  | 0.00664(10) |  | 0.00948(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |
O4 |  | 0.00854(10) |  | 0.00664(10) |  | 0.00948(11) |  | 0.00000 |  | 0.00000 |  | 0.00000 |  |

We also investigated two ways of possible Cu doping within our refinement: one
in which the Cu atom substitutes on the Pb lattice site, and another which the
Cu atom is inserted along the chain. Attempts to dope both Pb sites with Cu
results in a refined composition of Pb9.55Cu0.45(PO4)6(OH)0.94(SH)1.06 with
similar statistics (GoF$=1.21$, R $=3.93$). We found the Cu has to be added in
as a split site after Pb position has been refined completely. The last two
steps of occupancy and anisotropic parameter refinement had to be done by
restraining the Pb/Cu split site position with automatic refining keys
switched off.

We can also obtain a reasonable refinement to the data if we place Cu into the
channel with composition Pb10(PO4)6(OH) 1.11 Cu0.49 (GoF $=1.22$, R $=3.93$).
Still, as the crystals are transparent, charge balance needs to be maintained
and thus this last solution is not chemically reasonable.

Table S 12: Refinement Statistics for lead appetite structures with Cu doping. Refined Composition | Pb9.55Cu0.45(PO4)6(OH)0.94(SH)1.06 |  |  | Pb10(PO4)6(OH) 1.11 Cu0.49
---|---|---|---|---
Refinement Method | F2 |  |  | F2
Number of Parameters | 44 |  |  | 42
R(I>$3\sigma$), Rw(I>$3\sigma$) | 2.45, 5.41 |  |  | 2.46, 5.42
R(all), Rw(all) | 3.93, 5.67 |  |  | 3.94, 5.68
S(I>$3\sigma$), S(all) | 1.38, 1.21 |  |  | 1.39, 1.22
$\Delta\rho_{max}$ , $\Delta\rho_{min}$ (e Å-3) | 3.34, -3.09 |  |  | 3.38, -3.17

Table S 13: Refined atomic coordinates for Pb9.55Cu0.45(PO4)6(OH)0.94(SH)1.06 . Site |  | Wyckoff Position |  | x |  | y |  | z |  | Occupancy |  |
---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 6h |  | 0.24487 |  | 0.24674 |  | $1/4$ |  | 0.956 |  |
Cu1’ |  | 6h |  | 0.24487 |  | 0.24674 |  | $1/4$ |  | 0.044 |  |
Pb2 |  | 4f |  | $2/3$ |  | $1/3$ |  | -0.00344 |  | 0.953 |  |
Cu2’ |  | 4f |  | $2/3$ |  | $1/3$ |  | -0.00344 |  | 0.047 |  |
P1 |  | 6h |  | 0.3728(2) |  | 0.4019(2) |  | $3/4$ |  | 1 |  |
O1 |  | 6h |  | 0.4853(7) |  | 0.3370(8) |  | $3/4$ |  | 1 |  |
O2 |  | 12i |  | 0.2652(6) |  | 0.3480(6) |  | 0.5834(6) |  | 1 |  |
O3 |  | 6h |  | 0.5341(8) |  | 0.4168(7) |  | $1/4$ |  | 1 |  |
S1 |  | 2b |  | 0 |  | 0 |  | 0 |  | 0.53 |  |
O4 |  | 2a |  | 0 |  | 0 |  | $1/4$ |  | 0.47 |  |
Table S 14: Refined anisotropic displacement parameters for Pb9.55Cu0.45(PO4)6(OH)0.94(SH)1.06. Site |  | U11 |  | U22 |  | U33 |  | U12 |  | U13 |  | U23 |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 0.02041(17) |  | 0.01986(17) |  | 0.02671(17) |  | 0.01442(13) |  | 0.00000 |  | 0.00000 |  |
Cu1’ |  | 0.02041(17) |  | 0.01986(17) |  | 0.02671(17) |  | 0.01442(13) |  | 0.00000 |  | 0.00000 |  |
Pb2 |  | 0.01871(14) |  | 0.01871(14) |  | 0.01213(15) |  | 0.00936(7) |  | 0.00000 |  | 0.00000 |  |
Cu2’ |  | 0.01871(14) |  | 0.01871(14) |  | 0.01213(15) |  | 0.00936(7) |  | 0.00000 |  | 0.00000 |  |
P1 |  | 0.0100(8) |  | 0.0107(8) |  | 0.0138(7) |  | 0.0060(6) |  | 0.00000 |  | 0.00000 |  |
O1 |  | 0.024(3) |  | 0.029(3) |  | 0.026(3) |  | 0.021(3) |  | 0.00000 |  | 0.00000 |  |
O2 |  | 0.024(2) |  | 0.042(3) |  | 0.023(2) |  | 0.019(2) |  | -0.0067(17) |  | -0.0102(19) |  |
O3 |  | 0.031(3) |  | 0.011(2) |  | 0.035(3) |  | 0.006(2) |  | 0.00000 |  | 0.00000 |  |
S1 |  | 0.017(2) |  | 0.017(2) |  | 0.045(5) |  | 0.0084(12) |  | 0.00000 |  | 0.00000 |  |
O4 |  | 0.013(6) |  | 0.013(6) |  | 0.061(17) |  | 0.007(3) |  | 0.00000 |  | 0.00000 |  |
Table S 15: Refined atomic coordinates for Pb10(PO4)6(OH) 1.11 Cu0.49. Site |  | Wyckoff Position |  | x |  | y |  | z |  | Occupancy |  |
---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 6h |  | 0.24486(4) |  | 0.24667(4) |  | $1/4$ |  | 1 |  |
Pb2 |  | 4f |  | $2/3$ |  | $1/3$ |  | -0.00339(4) |  | 1 |  |
P1 |  | 6h |  | 0.3729(2) |  | 0.4020(2) |  | $3/4$ |  | 1 |  |
O1 |  | 6h |  | 0.4851(7) |  | 0.3366(8) |  | $3/4$ |  | 1 |  |
O2 |  | 12i |  | 0.2650(6) |  | 0.3480(6) |  | 0.5834(6) |  | 1 |  |
O3 |  | 6h |  | 0.5336(8) |  | 0.4165(7) |  | $1/4$ |  | 1 |  |
Cu1 |  | 2b |  | 0 |  | 0 |  | 0 |  | 0.247 |  |
O4 |  | 2a |  | 0 |  | 0 |  | $1/4$ |  | 0.55 |  |
Table S 16: Refined anisotropic displacement parameters for Pb10(PO4)6(OH) 1.11 Cu0.49. Site |  | U11 |  | U22 |  | U33 |  | U12 |  | U13 |  | U23 |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Pb1 |  | 0.02052(15) |  | 0.01995(15) |  | 0.02683(15) |  | 0.01447(12) |  | 0.00000 |  | 0.00000 |  |
Pb2 |  | 0.01891(12) |  | 0.01891(12) |  | 0.01231(13) |  | 0.00946(6) |  | 0.00000 |  | 0.00000 |  |
P1 |  | 0.0088(7) |  | 0.0095(7) |  | 0.0125(7) |  | 0.0054(6) |  | 0.00000 |  | 0.00000 |  |
O1 |  | 0.022(3) |  | 0.028(3) |  | 0.024(3) |  | 0.020(3) |  | 0.00000 |  | 0.00000 |  |
O2 |  | 0.022(2) |  | 0.040(3) |  | 0.0211(19) |  | 0.019(2) |  | -0.0065(17) |  | -0.0102(19) |  |
O3 |  | 0.029(3) |  | 0.010(2) |  | 0.034(3) |  | 0.006(2) |  | 0.00000 |  | 0.00000 |  |
Cu1 |  | 0.012(2) |  | 0.012(2) |  | 0.041(5) |  | 0.0062(12) |  | 0.00000 |  | 0.00000 |  |
O4 |  | 0.020(8) |  | 0.020(8) |  | 0.07(2) |  | 0.010(4) |  | 0.00000 |  | 0.00000 |  |

###  C.3 Powder Diffraction Comparisons

In this section, we compare the experimental data from [1, 2] with simulated
patterns of experimentally determined Fig. 32 and computational structures
Fig. 33 used in this body of work. The first two comparisons show patterns
from crystallographic information files not determined in this paper. Upon
first glance, it is understandable why the originals authors came to the
conclusion that their structure resembles Pb10(PO4)6O. The direct overlap of
peaks resembles a good match of the experimental pattern, however, as
previously stated this structure is unstable at room temperature.

In contrast, the simulated pattern for Pb10(PO4)6(OH)2 seems chemically
reasonable and utilizes neutron diffraction to refine the positions of the H
atoms, yet falls short in matching well with the data. The slow migration of
simulated peaks to higher 2$\theta$ values is indicative of the simulated
crystal structure needs to expand slightly to overlap with the experimental
data.

Comparing to the simulated patterns of 3 of our structure solutions reported
in this paper, we see negligible changes in the way our data fits the
experimental work. This is despite the fact that one of our structures,
Pb10(PO4)6(OH)Cu0.49 is chemically infeasible because it cannot be charged
balanced. This leads us to the conclusion that even though the patterns match
well, the match of experimental data is still not enough to accurately verify
composition or structure.

Figure 32: Comparisons of experimental (black) and calculated powder X-ray
diffraction patters of structures in this work (red). The calculated patterns
are all based on refinements of diffraction experiments.

In Fig. 33, we compare the simulated patterns of computationally relaxed
structures with the experimental data. Intriguingly, as we relax the
structures of Pb10(PO4)6O and Pb10(PO4)6(OH)2, both compounds expand, leading
to Pb10(PO4)6(OH)2 looking like a proper fit. This structure is then likely
another chemically reasonable interpretation of the data.

As we probed the Cu doped structures, we noticed that all of these
computationally modified structures grow another peak around 15 degrees which
varies in intensity based on the anion identity and Cu doping position shown
zoomed in with Fig. 34. The highest expected intensity at this angle is found
in the the Cu2-substituted Pb9Cu(PO4)6O pattern, which, except for ignoring
this peak is the best-matched pattern to the experimental data of the four.
Once again, we see crystallographic uncertainties with ordered Cu dopants.
From these simulations, we are skeptical of any ordering of Cu dopants inside
the lead apatite structure.

Another interesting possibility, however, is the prospect of the Cu dopants
not having a preferred site selectivity. We note that our model for
Pb9.55Cu0.45(PO4)6(OH)0.94 and Pb10(PO4)6(OH) 1.11 Cu0.49 have somewhat
significant Cu doping, yet they do not have any predicted peaks below 18
degrees.

In conclusion, laboratory X-ray diffraction techniques such as SCXD and PXRD
are not yet enough to come to a conclusion on the exact structure of lead
apatite synthesized via the Lanarkite plus Cu3P route, but hope that symbiotic
techniques such as neutron scattering and X-ray photoelectron spectroscopy
will be used in the near future.

Figure 33: Comparisons of previous [1, 2] experimental (black) and simulated
powder X-ray diffraction patters of relaxed structures in this work (red).
Figure 34: Zoomed in low angle experimental (black) and simulated powder X-ray
diffraction patters of selected structures in this work (red).

## References

  * Lee  _et al._ [2023a] S. Lee, J.-H. Kim, and Y.-W. Kwon, The First Room-Temperature Ambient-Pressure Superconductor (2023a), arxiv:2307.12008 [cond-mat] .
  * Lee  _et al._ [2023b] S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, Superconductor Pb10-xCux(PO4)6O showing levitation at room temperature and atmospheric pressure and mechanism (2023b), arxiv:2307.12037 [cond-mat] .
  * Abramian  _et al._ [2023] P. Abramian, A. Kuzanyan, V. Nikoghosyan, S. Teknowijoyo, and A. Gulian, Some remarks on possible superconductivity of composition Pb9CuP6O25 (2023), arxiv:2308.01723 [cond-mat] .
  * Baskaran [2023] G. Baskaran, Broad Band Mott Localization is all you need for Hot Superconductivity: Atom Mott Insulator Theory for Cu-Pb Apatite (2023), arxiv:2308.01307 [cond-mat] .
  * Brückner  _et al._ [1995a] S. Brückner, G. Lusvardi, L. Menabue, and M. Saladini, Inorganica Chimica Acta 236, 209 (1995a).
  * Cabezas-Escares  _et al._ [2023] J. Cabezas-Escares, N. F. Barrera, C. Cardenas, and F. Munoz, Theoretical insight on the LK-99 material (2023), arxiv:2308.01135 [cond-mat] .
  * Griffin [2023] S. M. Griffin, Origin of correlated isolated flat bands in copper-substituted lead phosphate apatite (2023), arxiv:2307.16892 [cond-mat] .
  * Guo  _et al._ [2023] K. Guo, Y. Li, and S. Jia, Ferromagnetic half levitation of LK-99-like synthetic samples (2023), arxiv:2308.03110 [cond-mat] .
  * Hirschmann and Mitscherling [2023] M. M. Hirschmann and J. Mitscherling, Tight-binding models for SG 143 (P3) and application to recent DFT results on copper-doped lead apatite (2023), arxiv:2308.03751 [cond-mat] .
  * Hou  _et al._ [2023] Q. Hou, W. Wei, X. Zhou, Y. Sun, and Z. Shi, Observation of zero resistance above 100∘ K in Pb10-xCux(PO4)6O (2023), arxiv:2308.01192 [cond-mat] .
  * Kumar  _et al._ [2023a] K. Kumar, N. K. Karn, and V. P. S. Awana, Synthesis of possible room temperature superconductor LK-99:Pb9Cu(PO4)6O (2023a), arxiv:2307.16402 [cond-mat] .
  * Kumar  _et al._ [2023b] K. Kumar, N. K. Karn, Y. Kumar, and V. P. S. Awana, Absence of superconductivity in LK-99 at ambient conditions (2023b), arxiv:2308.03544 [cond-mat] .
  * Kurleto  _et al._ [2023] R. Kurleto, S. Lany, D. Pashov, S. Acharya, M. van Schilfgaarde, and D. S. Dessau, Pb-apatite framework as a generator of novel flat-band CuO based physics, including possible room temperature superconductivity (2023), arxiv:2308.00698 [cond-mat] .
  * Lai  _et al._ [2023] J. Lai, J. Li, P. Liu, Y. Sun, and X.-Q. Chen, First-principles study on the electronic structure of Pb10-xCux(PO4)6O ($x$=0, 1) (2023), arxiv:2307.16040 [cond-mat] .
  * Liu  _et al._ [2023] L. Liu, Z. Meng, X. Wang, H. Chen, Z. Duan, X. Zhou, H. Yan, P. Qin, and Z. Liu, Semiconducting transport in Pb10-xCux(PO4)6O sintered from Pb2SO5 and Cu3P (2023), arxiv:2307.16802 [cond-mat, physics:physics] .
  * Oh and Zhang [2023] H. Oh and Y.-H. Zhang, S-wave pairing in a two-orbital t-J model on triangular lattice: Possible application to Pb10-xCux(PO4)6O (2023), arxiv:2308.02469 [cond-mat] .
  * Si and Held [2023] L. Si and K. Held, Electronic structure of the putative room-temperature superconductor Pb9Cu(PO4)6O (2023), arxiv:2308.00676 [cond-mat] .
  * Sun  _et al._ [2023] Y. Sun, K.-M. Ho, and V. Antropov, Metallization and Spin Fluctuations in Cu-doped Lead Apatite (2023), arxiv:2308.03454 [cond-mat] .
  * Tao  _et al._ [2023] K. Tao, R. Chen, L. Yang, J. Gao, D. Xue, and C. Jia, The Cu induced ultraflat band in the room-temperature superconductor Pb10-xCux(PO4)6O4 ($x=0,0.5$) (2023), arxiv:2308.03218 [cond-mat] .
  * Tavakol and Scaffidi [2023] O. Tavakol and T. Scaffidi, Minimal model for the flat bands in copper-substituted lead phosphate apatite (2023), arxiv:2308.01315 [cond-mat] .
  * Wu  _et al._ [2023] H. Wu, L. Yang, B. Xiao, and H. Chang, Successful growth and room temperature ambient-pressure magnetic levitation of LK-99 (2023), arxiv:2308.01516 [cond-mat] .
  * Peotta and Törmä [2015] S. Peotta and P. Törmä, Nature Communications 6, 8944 (2015), arXiv:1506.02815 [cond-mat.supr-con] .
  * Törmä  _et al._ [2022] P. Törmä, S. Peotta, and B. A. Bernevig, Nature Reviews Physics 4, 528 (2022), arXiv:2111.00807 [cond-mat.supr-con] .
  * Zhou  _et al._ [2021] X. Zhou, W.-S. Lee, M. Imada, N. Trivedi, P. Phillips, H.-Y. Kee, P. Törmä, and M. Eremets, Nature Reviews Physics 3, 462 (2021).
  * Tovmasyan  _et al._ [2018] M. Tovmasyan, S. Peotta, L. Liang, P. Törmä, and S. D. Huber, Phys. Rev. B 98, 134513 (2018), arXiv:1805.04529 [cond-mat.str-el] .
  * Tovmasyan  _et al._ [2016] M. Tovmasyan, S. Peotta, P. Törmä, and S. D. Huber, Phys. Rev. B 94, 245149 (2016), arXiv:1608.00976 [cond-mat.str-el] .
  * Liang  _et al._ [2017] L. Liang, T. I. Vanhala, S. Peotta, T. Siro, A. Harju, and P. Törmä, Phys. Rev. B 95, 024515 (2017), arXiv:1610.01803 [cond-mat.supr-con] .
  * Huhtinen  _et al._ [2022] K.-E. Huhtinen, J. Herzog-Arbeitman, A. Chew, B. A. Bernevig, and P. Törmä, Phys. Rev. B 106, 014518 (2022).
  * Hu  _et al._ [2019] X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Phys. Rev. Lett. 123, 237002 (2019).
  * Jiang and Barlas [2023] G. Jiang and Y. Barlas, Phys. Rev. Lett. 131, 016002 (2023).
  * Mao and Chowdhury [2023] D. Mao and D. Chowdhury, Proceedings of the National Academy of Science 120, e2217816120 (2023), arXiv:2209.06817 [cond-mat.str-el] .
  * Hofmann  _et al._ [2023] J. S. Hofmann, E. Berg, and D. Chowdhury, Phys. Rev. Lett. 130, 226001 (2023), arXiv:2204.02994 [cond-mat.str-el] .
  * Törmä  _et al._ [2018] P. Törmä, L. Liang, and S. Peotta, Phys. Rev. B 98, 220511 (2018), arXiv:1810.09870 [cond-mat.supr-con] .
  * Herzog-Arbeitman  _et al._ [2022a] J. Herzog-Arbeitman, A. Chew, K.-E. Huhtinen, P. Törmä, and B. A. Bernevig, arXiv e-prints , arXiv:2209.00007 (2022a), arXiv:2209.00007 [cond-mat.str-el] .
  * Yu  _et al._ [2022] J. Yu, Y.-A. Chen, and S. Das Sarma, Phys. Rev. B 105, 104515 (2022).
  * Yu  _et al._ [2023] J. Yu, M. Xie, F. Wu, and S. Das Sarma, Phys. Rev. B 107, L201106 (2023).
  * Xie  _et al._ [2020] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Phys. Rev. Lett. 124, 167002 (2020), arXiv:1906.02213 [cond-mat.supr-con] .
  * Herzog-Arbeitman  _et al._ [2022b] J. Herzog-Arbeitman, V. Peri, F. Schindler, S. D. Huber, and B. A. Bernevig, Phys. Rev. Lett. 128, 087002 (2022b), arXiv:2110.14663 [cond-mat.mes-hall] .
  * Cǎlugǎru  _et al._ [2022] D. Cǎlugǎru, A. Chew, L. Elcoro, Y. Xu, N. Regnault, Z.-D. Song, and B. A. Bernevig, Nature Physics 18, 185 (2022), arXiv:2106.05272 [cond-mat.mes-hall] .
  * Hwang  _et al._ [2021] Y. Hwang, J.-W. Rhim, and B.-J. Yang, Phys. Rev. B 104, L081104 (2021).
  * Hwang  _et al._ [2021] Y. Hwang, J.-W. Rhim, and B.-J. Yang, Phys. Rev. B 104, 085144 (2021), arXiv:2106.13057 [cond-mat.mes-hall] .
  * Resta [2011] R. Resta, European Physical Journal B 79, 121 (2011), arXiv:1012.5776 [cond-mat.mtrl-sci] .
  * Lieb [1989] E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * Mielke and Tasaki [1993] A. Mielke and H. Tasaki, Communications in Mathematical Physics 158, 341 (1993), arXiv:cond-mat/9305026 [cond-mat] .
  * Krivovichev and Burns [2003] S. V. Krivovichev and P. C. Burns, Zeitschrift für Kristallographie-Crystalline Materials 218, 357 (2003).
  * Brückner  _et al._ [1995b] S. Brückner, G. Lusvardi, L. Menabue, and M. Saladini, Inorganica Chimica Acta 236, 209 (1995b).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Elcoro  _et al._ [2021] L. Elcoro, B. J. Wieder, Z. Song, Y. Xu, B. Bradlyn, and B. A. Bernevig, Nature communications 12, 5965 (2021).
  * Zhu  _et al._ [2023] S. Zhu, W. Wu, Z. Li, and J. Luo, arXiv e-prints , arXiv:2308.04353 (2023), arXiv:2308.04353 [cond-mat.supr-con] .
  * Tinkham [2004] M. Tinkham,  _Introduction to Superconductivity_, 2nd ed. (Dover Publications, 2004).
  * Yu  _et al._ [2016] F. Yu, M. Hirschberger, T. Loew, G. Li, B. J. Lawson, T. Asaba, J. B. Kemper, T. Liang, J. Porras, G. S. Boebinger, J. Singleton, B. Keimer, L. Li, and N. P. Ong, Proceedings of the National Academy of Sciences 113, 12667 (2016), https://www.pnas.org/doi/pdf/10.1073/pnas.1612591113 .
  * Liao  _et al._ [1999] C.-J. Liao, F.-H. Lin, K.-S. Chen, and J.-S. Sun, Biomaterials 20, 1807 (1999).
  * Negas and Roth [1968] T. Negas and R. Roth, Journal of Research of the National Bureau of Standards. Section A, Physics and Chemistry 72, 783 (1968).
  * Verbeeck  _et al._ [1980] R. M. H. Verbeeck, H. J. M. Heiligers, F. C. M. Driessens, and H. G. Schaeken, Zeitschrift für anorganische und allgemeine Chemie 466, 76 (1980), https://onlinelibrary.wiley.com/doi/pdf/10.1002/zaac.19804660109 .
  * Freysoldt  _et al._ [2014] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, Rev. Mod. Phys. 86, 253 (2014).
  * Walsh [2021] A. Walsh, npj Computational Materials 7, 72 (2021).
  * Freysoldt  _et al._ [2009] C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Rev. Lett. 102, 016402 (2009).
  * Bian  _et al._ [2022] H. Bian, T. Liu, D. Li, Z. Xu, J. Lian, M. Chen, J. Yan, and S. F. Liu, Chemical Engineering Journal 435, 135071 (2022).
  * Cerny  _et al._ [1994] R. Cerny, K. Yvon, M. Wakihara, and P. Fischer, Journal of Alloys and Compounds 209, L29 (1994).
  * Zhu  _et al._ [2023] S. Zhu, W. Wu, Z. Li, and J. Luo, First order transition in pb10-xcux(po4)6o ($0.9<x<1.1$) containing cu2s (2023), arXiv:2308.04353 [cond-mat.supr-con] .
  * Jain [2023] P. K. Jain, Phase transition of copper (i) sulfide and its implication for purported superconductivity of lk-99 (2023), arXiv:2308.05222 [cond-mat.supr-con] .
  * Inc. [2012] B. A. Inc., Apex2 (2012), accessed August 7, 2023.
  * Petricek  _et al._ [2023] V. Petricek, L. Palatinus, J. Plasil, and M. Dusek, Jana2020 – a new version of the crystallographic computing system jana (2023).
  * Palatinus L. [2007] C. G. Palatinus L., Superflip – a computer program for the solution of crystal structures by charge flipping in arbitrary dimensions. (2007).
  * Momma and Izumi [2008] K. Momma and F. Izumi, Vesta: a three-dimensional visualization system for electronic and structural analysis (2008).
  * Vil [a] Pb10(po4)6o (pb10[po4]6o) crystal structure: Datasheet from “pauling file multinaries edition – 2022” in springermaterials (https://materials.springer.com/isp/crystallographic/docs/sd_1616772) (a), copyright 2023 Springer-Verlag Berlin Heidelberg & Material Phases Data System (MPDS), Switzerland & National Institute for Materials Science (NIMS), Japan.
  * Vil [b] Pb10(po4)6(oh)2 (pb5[po4]3[oh]) crystal structure: Datasheet from “pauling file multinaries edition – 2022” in springermaterials (https://materials.springer.com/isp/crystallographic/docs/sd_1003884) (b), copyright 2023 Springer-Verlag Berlin Heidelberg & Material Phases Data System (MPDS), Switzerland & National Institute for Materials Science (NIMS), Japan.
  * Vil [c] Pb10(po4)6(oh)2 (pb5[po4]3[oh]) crystal structure: Datasheet from “pauling file multinaries edition – 2022” in springermaterials (https://materials.springer.com/isp/crystallographic/docs/sd_1714473) (c), copyright 2023 Springer-Verlag Berlin Heidelberg & Material Phases Data System (MPDS), Switzerland & National Institute for Materials Science (NIMS), Japan.
  * Zhu  _et al._ [2007] K. Zhu, K. Yanagisawa, R. Shimanouchi, A. Onda, K. Kajiyoshi, and J. Qiu, Journal of the Ceramic Society of Japan 115, 873 (2007).
  * Vil [d] Pb10(po4)6(oh)2 (pb5[po4]3[oh]) crystal structure: Datasheet from “pauling file multinaries edition – 2022” in springermaterials (https://materials.springer.com/isp/crystallographic/docs/sd_1631620) (d), copyright 2023 Springer-Verlag Berlin Heidelberg & Material Phases Data System (MPDS), Switzerland & National Institute for Materials Science (NIMS), Japan.
  * Zhu  _et al._ [2010] K. Zhu, J. Qiu, H. Ji, K. Yanagisawa, R. Shimanouchi, A. Onda, and K. Kajiyoshi, Inorganica Chimica Acta 363, 1785 (2010).
  * Kim  _et al._ [1997] J. Y. Kim, B. A. Hunter, R. R. Fenton, and B. J. Kennedy, Australian journal of chemistry 50, 1061 (1997).
  * Aryasetiawan  _et al._ [2004] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. Lichtenstein, Physical Review B 70, 195104 (2004).
  * Solovyev and Imada [2005] I. Solovyev and M. Imada, Physical Review B 71, 045103 (2005).
  * Aryasetiawan  _et al._ [2006] F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schönberger, Physical Review B 74, 125106 (2006).
  * Miyake  _et al._ [2009] T. Miyake, F. Aryasetiawan, and M. Imada, Physical Review B 80, 155134 (2009).
  * Yu  _et al._ [2023] J. Yu, C. J. Ciccarino, R. Bianco, I. Errea, P. Narang, and B. A. Bernevig, arXiv e-prints , arXiv:2305.02340 (2023), arXiv:2305.02340 [cond-mat.supr-con] .
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Computational materials science 6, 15 (1996a).
  * Kresse and Hafner [1993a] G. Kresse and J. Hafner, Physical Review B 48, 13115 (1993a).
  * Kresse and Hafner [1993b] G. Kresse and J. Hafner, Physical review B 47, 558 (1993b).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Physical Review B 49, 14251 (1994).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Physical review B 54, 11169 (1996b).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Physical review B 56, 12847 (1997).
  * Souza  _et al._ [2001] I. Souza, N. Marzari, and D. Vanderbilt, Physical Review B 65, 035109 (2001).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Reviews of Modern Physics 84, 1419 (2012).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, _et al._ , Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Togo  _et al._ [2008] A. Togo, F. Oba, and I. Tanaka, Phys. Rev. B 78, 134106 (2008).