On the dynamical stability of copper-doped lead apatite

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


INTRODUCTION
Copper-doped lead apatite Pb 10−x Cu x (PO 4 ) 6 O with 0.9 < x < 1.1, known as LK-99, has been recently claimed to exhibit superconductivity above room temperature and at ambient pressure [1,2].This remarkable claim is backed by magnetic (half-)levitation on a permanent magnet and by a sudden drop in resistivity at the claimed superconducting transition temperature.However, subsequent extensive experimental efforts by other groups have failed to confirm the superconductivity [3][4][5][6][7][8][9][10][11][12][13].The magnetic half-levitation is reproduced in some insulating samples, where it is attributed to soft ferromagnetism [6][7][8].A plausible explanation for the sudden resistivity drop is provided by a first order phase transition of Cu 2 S impurities [9,14], which is further supported by the highly insulating nature of a single crystalline sample without Cu 2 S impurities [13].
On the theoretical front, initial density functional theory calculations reported an electronic structure exhibiting relatively flat bands near the Fermi level for a simple model of copper-doped lead apatite [15][16][17][18].However, subsequent calculations showed that the inclusion of spin-orbit coupling or non-local correlations lead to an insulating electronic structure [19][20][21], a conclusion that is also reached with the inclusion of local correlations using dynamical mean-field theory [22][23][24].Different estimates of critical superconducting temperatures have so far delivered values significantly lower than room temperature [25][26][27].
These state-of-the-art electronic structure calculations all assume a specific structural model as a starting point, often suggested by experiment.However, other theoretical works have questioned the suitability of these structural models both in terms of the thermodynamic feasibility of copper doping [28] or the dynamical stability of the experimentally proposed structures [28][29][30][31][32][33].Indeed, one of the most basic quantities used to characterize a material is its dynamical stability.A dynamically stable structure corresponds to a local minimum of the potential (free) energy surface, and its phonon frequencies are real.A dynamically unstable structure corresponds to a saddle point of the potential (free) energy surface, and some of its phonon frequencies are imaginary with associated eigenvectors that encode atomic displacement patterns that lower the energy of the system.Only dynamically stable structures can represent real materials.Puzzlingly, recent computational works have claimed that the experimentally reported structures of the parent lead apatite [28][29][30][31] and of the copper-doped lead apatite compounds [28-30, 32, 33] are dynamically unstable at the harmonic level, which would imply that they cannot be the true structures of the materials underpinning LK-99, and would question the validity of most electronic structure calculations to date.
In this work, we demonstrate that both parent lead apatite and copper-doped lead apatite compounds are dynamically stable at room temperature.The parent compounds are largely stable at the harmonic level, with some exhibiting very slight instabilities which are suppressed by quartic anharmonic terms.For the copper-doped compounds, dynamical stability at the harmonic level depends on the doping site, but even those that are dynamically unstable at the harmonic level are overall stable at room temperature with the inclusion of anharmonic phonon-phonon interactions.

Lead apatite
Lead apatite is a compound that was first experimentally reported over 70 years ago [34].Figure 1 depicts an example of lead apatite, with a hexagonal lattice (space group P 6 3 /m) and general formula Pb 10 (PO 4 ) 6 X 2 , where X is either a halide atom or an OH group.The variant Pb 10 (PO 4 ) 6 O, which is claimed to be the parent structure of LK-99 [1,2], has also been known experimentally for decades [34][35][36][37].The X site corresponds to Wyckoff position 4e, giving a multiplicity of four in the unit cell, but these sites are only partially filled with an occupation of 1  2 for halide atoms and the OH group, and an occupation of 1  4 for O.We consider two representative cases, Pb 10 (PO 4 ) 6 O and Pb 10 (PO 4 ) 6 (OH) 2 , where the specific distribution of species on the X site results in space groups P 3 and P 6 3 , respectively.
The phonon dispersion of the parent Pb 10 (PO 4 ) 6 O and Pb 10 (PO 4 ) 6 (OH) 2 compounds is shown in Fig. 1.At the harmonic level, Pb 10 (PO 4 ) 6 O exhibits imaginary phonon frequencies at the zone boundary points M and K of the k z = 0 plane, and at the zone boundary point H of the k z = π c plane.However, the absolute values of these imaginary frequencies are less than 2.3 meV and the resulting anharmonic potentials have a dominant quartic term that strongly suppresses the instability to about 0.2 meV per formula unit (see Supplementary Figure 4).As a result, the calculation of self-consistent phonons including anharmonic interactions fully stabilizes the structure at the relatively low temperature of 50 K, and potentially lower.Earlier works reported that Pb 10 (PO 4 ) 6 O is dynamically unstable at the harmonic level [28][29][30], a result we confirm, but our work further demonstrates that Pb 10 (PO 4 ) 6 O is overall dynamically stable at room temperature driven by higher-order anharmonic terms.
Pb 10 (PO 4 ) 6 (OH) 2 is dynamically stable at the harmonic level.Earlier works reported that Pb 10 (PO 4 ) 6 (OH) 2 is dynamically unstable [28,29], an opposite conclusion that we attribute to unconverged harmonic calculations (see Supplementary Figure 3).We note that to fully converge the harmonic calculations of Pb 10 (PO 4 ) 6 (OH) 2 , the required coarse q-point grid includes all of Γ, M, K, and A points, which can only be accomplished with a regular grid of minimum size 6 × 6 × 2 or alternatively a nonuniform Farey grid [38] of minimum size (2 In our calculations, we use the latter as it is computationally more efficient.We highlight that our strategy, combining a nonuniform Farey grid [38] with nondiagonal supercells [39,40], offers a significant computational advantage in phonon calculations using the finite displacement method compared to the conventional diagonal supercell method with a regular grid used by most earlier works.Our approach drastically reduces the supercell sizes required to reach convergence, with the number of atoms decreasing from 3168 to only 132. Overall, we find that both lead apatite compounds Pb 10 (PO 4 ) 6 O and Pb 10 (PO 4 ) 6 (OH) 2 are dynamically stable, a conclusion that is in full agreement with multiple experimental reports of the structure of lead apatite over the past 70 years [34][35][36][37]41].
Copper doped lead apatite The claim of room temperature superconductivity in LK-99 is based on copper doping of lead apatite, with copper replacing about 1 in 10 lead atoms leading to a Pb 9 Cu(PO 4 ) 6 O stoichiometry.
There are two symmetrically distinct lead sites, labelled Pb(1) and Pb(2) in the literature (see Fig. 1a), and doping at these sites results in structures with the space groups P 3 and P 1, respectively (Fig. 2a,b).The original LK-99 work suggested that the doping site is Pb(1) [1,2], but subsequent experimental works have suggested that both Pb(1) and Pb(2) sites can be doped [13,29].
Computational works find that the relative energy between the two doping sites depends on the exchange-correlation functional and the magnitude of the Hubbard U parameter used on the copper atom, with most choices favouring doping at the Pb(2) site, a prediction we confirm with our own calculations.For completeness, in this work we explore doping at both sites.
The electronic structures of Pb 9 Cu(PO 4 ) 6 O with copper on the Pb(1) and Pb(2) sites are shown in Fig. 2c,d.For doping at the Pb(1) site, a non-magnetic calculation leads to a metallic state in which the Fermi energy crosses four relatively flat bands (a pair of doubly-degenerate bands).Inclusion of spin-orbit coupling while maintaining the non-magnetic configuration leads to a splitting of the pair of doubly-degenerate bands, and the Fermi energy crosses a pair of singly-degenerate relatively flat bands.A calculation including spin-orbit coupling and allowing a non-zero magnetic moment leads to a ferromagnetic configuration in which the system is gapped.We note that the spin-orbit coupling is not essential for a gap opening in the presence of ferromagnetic ordering.
The actual role of spin-orbit coupling, which lifts the orbital degeneracy, can be replicated by selecting a suitable level of theory, as evidenced by hybrid [20] or GW [21] calculations showcasing a gapped state in the presence of ferromagnetic ordering.The ferromagnetic configuration is the most energetically favourable, but may not be directly relevant for room temperature experiments as single crystal measurements suggest the material is a non-magnetic insulator exhibiting a diamagnetic response with potentially a small ferromagnetic component [13].Additionally, the ferromagnetic ordering may be an artifact of the DFT calculations, as dynamical mean-field theory calculations [22][23][24] suggest a gap opens due to a Mott-like band splitting without the need of ferromagnetic ordering.For doping at the Pb(2) site, we also find a metallic state with a single band crossing the Fermi level in non-magnetic calculations, and a gapped state in ferromagnetic calculations.For both doping sites, we find that the phonon dispersion is only weakly affected by the level of electronic structure theory used (see Supplementary Notes 3.1 and 3.2), so the discussion below should be largely independent of the precise electronic structure of the system.
We hereafter present the phonon dispersions of the non-magnetic state calculated using PBEsol+U without spin-orbit coupling.
The phonon dispersions of Pb 9 Cu(PO 4 ) 6 O with copper on the Pb(1) and Pb(2) sites are shown in Fig. 2e,f.Doping at the Pb(2) site leads to a dynamically stable structure at the harmonic level of theory.By contrast, doping at the Pb(1) site leads to a dynamically unstable structure at the harmonic level that exhibits two imaginary phonon branches of frequencies about 15i meV across the entire Brillouin zone.This harmonic instability is present irrespective of the level of theory used, including a Hubbard U parameter on the copper d orbitals, spin-orbit coupling, and ferromagnetic ordering (see Supplementary Figure 5).Importantly, anharmonic phonon-phonon interactions strongly suppress the instability and the structure becomes dynamically stable at 300 K. We reach similar conclusions for copper doping of Pb 10 (PO 4 ) 6 (OH) 2 (see Supplementary Figure 6).Overall, copper-doped lead apatite is dynamically stable at room temperature for doping at either site.The original paper claiming superconductivity in LK-99 suggested that copper doping of lead apatite occurs on the Pb(1) site [1,2].As the associated structure exhibits a dynamical instability at the harmonic level, we further explore its properties by considering the potential energy surface along the imaginary phonon modes at high symmetry points in the Brillouin zone (Fig. 3a).The dominant instability is driven by a Γ point phonon mode, and fully relaxing the structure along this instability leads to a distinct structure of P 1 symmetry (Fig. 3b; see also the link in Data availability statement for the optimized structure file).We find that the P 1 structure is dynamically stable at the harmonic level of theory (Fig. 3c), consistent with an earlier report [33].We ascribe the harmonic stability of the P 1 structure to a downward shift of the occupied part of the density of states compared to the P 3 structure (Fig. 3e), a trend similarly observed in the density of states of the harmonically stable Pb(2) doping case in Fig. 2f (see also Supplementary Figure 10 for the density of states of the Pb(2) doping case).This indicates that the harmonic stability is dominated by copper-derived orbitals.The four bands (a pair of doubly-degenerate bands) that cross the Fermi level in the P 3 structure (Fig. 2c) split under the distortion, such that the resultant P 1 structure has a metallic state with a single doubly-degenerate band crossing the Fermi level in the non-magnetic configuration (Fig. 3d).The distorted P 1 structure becomes an insulator in the presence of ferromagnetic ordering, which is consistent with previous results [33,42], similar to lead apatite with copper doping at the Pb(2) site (Fig. 2d).
Interestingly, the relative energy of the P 3 structure compared to the Γ-distorted P 1 structure is strongly dependent on both the volume and the electronic correlation strength as measured by the Hubbard U parameter.Specifically, we find that harmonic instabilities favouring the P 1 phase occur for large values of U and large volumes, while the harmonic instability of the P 3 phase completely disappears for small values of U and small volumes.This is evident from the phonon dispersion of the P 3 phase for different U values and volumes (Fig. 4a) and from the enthalpy difference between P 3 and P 1 phases indicated by the colour bar in the phase diagram in Fig. 4b.
These observations suggest that controlling volume, for example through hydrostatic pressure or strain, and controlling the degree of electronic correlation, for example by applying a gate voltage or doping, can be used to navigate the structural phase diagram of compounds based on lead apatite.Specifically, it may be possible to observe a temperature-driven structural phase transition between a low temperature P 1 phase and high temperature P 3 phase in a regime with a large harmonic dynamical instability.Finally, we note that electronic correlation beyond the static description provided by a Hubbard U correction may play an important role on this phase diagram [22][23][24], so further work is required to fully characterise it.

DISCUSSION
Since the original claim of room temperature superconductivity in LK-99, seven phonon dispersion calculations have been reported in the literature.Of these, parent [28][29][30][31] and copper-doped lead apatite with a P 3 space group [28-30, 32, 33] are claimed to be dynamically unstable at the harmonic level, while another work claims that the copper-doped lead apatite is dynamically stable at the harmonic level [27].
We attribute these puzzling and contradictory conclusions about the dynamical stability of lead apatite to the complexity of harmonic phonon calculations in this system, with a unit cell containing at least 41 atoms, and to the subtle interplay between volume, electronic correlation strength, and phonons.First, we find that fully converged phonon calculations for the parent compound Pb 10 (PO 4 ) 6 (OH) 2 require relatively large coarse q-point grids, specifically incorporating the Γ, M, K, and A points.However, none of the previously reported calculations include all these points, and as a result they incorrectly conclude that this compound is dynamically unstable at the harmonic level.Second, we find that dynamical stability for the copper-doped compounds at the harmonic level depends on both the value of the Hubbard U parameter and the volume of the system (Fig. 4; see also Supplementary Note 3.3).In this context, we rationalise the seemingly contradictory conclusions about dynamical stability of the copper-doped compounds by suggesting that different works use different volumes and different choices for the Hubbard U parameter.
Beyond clarifying the dynamical stability of lead apatite at the harmonic level, we have shown that anharmonic phonon-phonon interactions play a key role in stabilising multiple lead apatite compounds.Overall, our calculations indicate that both parent and copper-doped lead apatite compounds are dynamically stable at room temperature.
We believe that lead apatite is a nice example to illustrate the ability of state-of-the-art first principles methods to fully characterise a complex and experimentally relevant system.However, our work also demonstrates that reliable results and conclusions can only be reached with a careful consideration of convergence parameters, such as the size of the q-point grid, and physical models, such as the inclusion of anharmonic phonon-phonon interactions.
We show that the experimentally suggested structures of lead apatite and copper-doped lead apatite are dynamically stable at room temperature.Most structures are dynamically stable at the harmonic level, but some key structures, including the structure claimed to be responsible for superconductivity at ambient conditions, only becomes dynamically stable with the inclusion of anharmonic phonon-phonon interactions.Our results resolve a puzzling suggestion by multiple earlier computational works that claimed that the experimentally reported structures of both parent and copper-doped lead apatite compounds were dynamically unstable, and fully reconcile the current experimental and theoretical description of the structure of lead apatite.
We find that experimental lattice parameters agree well with those predicted by PBEsol, and the data presented in the main text has been obtained using PBEsol (see comparison between PBE and PBEsol in Supplementary Figure 5a).An on-site Hubbard interaction U is applied to the copper 3d orbitals based on the simplified rotationally invariant DFT+U method by Dudarev and co-workers [48].We have checked that DFT+U gives almost identical lattice parameters to DFT.
Converged results are obtained with a kinetic energy cutoff for the plane wave basis of 600 eV and a k-point grid of size 4 × 4 × 5 and 6 × 6 × 8 for the primitive cell of the parent and copper-doped lead apatite, respectively (see convergence test results in Supplementary Note 1).The geometry of the structures is optimised until all forces are below 0.01 eV per Å and the pressure is below 1 kbar.
Harmonic phonon calculations.-We perform harmonic phonon calculations using the finite displacement method in conjunction with nondiagonal supercells [39,40].We find that a coarse q-point grid of size 2 × 2 × 2 leads to converged phonon dispersions for the parent compound Pb 10 (PO 4 ) 6 O and the Cu-doped compounds.However, for the parent compound Pb 10 (PO 4 ) 6 (OH) 2 , a converged calculation requires a minimum coarse q-point grid including the high symmetry points Γ, M, K, and A, which we accomplish by means of a Farey nonuniform grid [38] of size (2×2×2)∪(3×3×1).
To evaluate the force derivatives, we use a three-point central formula with a finite displacement of 0.02 bohr.The underlying electronic structure calculations are performed using the same parameters as those described above.We have also cross-checked the phonon dispersions by performing additional calculations with castep [49] and Quantum Espresso [50] within the finite displacement method and density functional perturbation theory, respectively (see Supplementary Figure 2).
Anharmonic phonon calculations.-We perform anharmonic phonon calculations using the stochastic self-consistent harmonic approximation (SSCHA) [51][52][53], which accounts for anharmonic effects at both zero and finite temperature.The self-consistent harmonic approximation [54] is a quantum variational method on the free energy, and the variational minimization is performed with respect to a trial harmonic system.In its stochastic implementation, the forces on atoms are calculated in an ensemble of configurations drawn from the trial harmonic system.We use vasp to perform electronic structure calculations using the same parameters as those described above, and consider configurations commensurate with a 2 × 2 × 2 supercell.The number of configurations needed to converge the free energy Hessian is of the order of 4,000 for the parent lead apatite structure and of the order of 8,000 configurations for the copper-doped structure.

CULATIONS
We have tested various convergence parameters for the electronic structure calculations underpinning the calculation of the phonon dispersions.We find that an energy cutoff of 600 eV and a k-point grid size of 4 × 4 × 5 are converged for the parent lead apatite We have also performed a cross-check of the phonon dispersion using three different codes: vasp [1,2], castep [3] and Quantum Espresso [4].We have confirmed that these codes yield very similar results, as depicted in Fig. 2. The vasp and castep calculations are done using the finite displacement method, while the quantum espresso calculations are done using density functional perturbation theory.The PBEsol exchange-correlation functional and a 2 × 2 × 2 q-point grid are used.

CASTEP
Figure 3 shows the phonon dispersions of Pb 10 (PO 4 ) 6 O (top) and Pb 10 (PO 4 ) 6 (OH) 2 (bottom) for different choices of coarse q-point grid size.For Pb 10 (PO 4 ) 6 O, a qualitatively correct phonon dispersion is obtained with a coarse grid of size 2 × 2 × 2, consistent with a previous study [5].By employing larger grid sizes, we confirm that the imaginary frequencies at the K, M, and H points are physical rather than an artifact of Fourier interpolation.Supplementary Figure 3: Coarse q-point grid dependence of the harmonic phonon dispersions.a,b.Harmonic phonon dispersions of a Pb 10 (PO 4 ) 6 O and b Pb 10 (PO 4 ) 6 (OH) 2 using different coarse q-point grids.The maximum number of atoms used in each q-point grid within the finite displacement method in conjunction with nondiagonal supercells [6,7] is specified.The calculations are performed using the PBEsol exchange-correlation functional.
For Pb 10 (PO 4 ) 6 (OH) 2 , a qualitatively correct phonon dispersion is only obtained when all of Γ, M, K, and A points are included in the coarse q-point grid, which can only be accomplished with a regular grid of minimum size 6 × 6 × 2 or alternatively a non-uniform Farey grid [8] of minimum size (2×2×2)∪(3×3×1).We use the latter as it is computationally more efficient: within the nondiagonal supercell formalism [6,7], it requires calculations using supercells with up to 132 atoms, compared to up to the 264 atoms required for the regular 6 × 6 × 2 grid.We also highlight that our strategy, employing the nonuniform Farey grid with nondiagonal supercells, offers significant advantages in phonon calculations compared to the conventional diagonal supercell method with a regular grid.This approach drastically reduces the maximum number of atoms required for a supercell, from 3168 to only 132, enabling us to obtain the converged phonon dispersion for Pb 10 (PO 4 ) 6 (OH) 2 .
We note that all earlier phonon calculations for the Pb 10 (PO 4 ) 6 (OH) 2 compound in the literature use coarse q-point grids of sizes 1 × 1 × 1 [9] or 1 × 1 × 2 [10], and the imaginary phonon modes observed in these calculations are not physical but instead an artifact of Fourier interpolation caused by unconverged calculations.Indeed, the phonon dispersions reported in these works coincide with the corresponding unconverged calculations depicted in Fig. 3b.The harmonic phonon dispersion of Pb 10 (PO 4 ) 6 O exhibits imaginary frequencies at the M, K, and H points of the Brillouin zone, as depicted in the left panel of Fig. 4. The absolute values of these imaginary frequencies are below 2.3 meV, suggesting that the system is only marginally unstable at the harmonic level.To confirm this, we calculate the potential energy surface by displacing the atoms along the eigenvectors associated with the three imaginary modes, resulting in the double well potentials shown in the right panel of Fig. 4.
The anharmonic potentials have a dominant quartic term that suppresses the instability to a maximum of 0.2 meV per formula unit for the K point, with smaller instabilities for the M and H points.By comparison, the thermal energy associated with room temperature is about 26 meV, suggesting that these shallow double well potentials can be overcome by thermally-induced anharmonic vibrations.We confirm this in the main text by performing self-consistent harmonic calculations at 50 K that deliver a dynamically stable phonon dispersion.The structure is likely dynamically stable at even lower temperatures, possibly at 0 K where it would be stabilized by quantum fluctuations.similar dependence on the Hubbard U parameter, in which the magnitude of the imaginary frequency increases with increasing U , and a similarly robust imaginary frequency with the inclusion of spin-orbit coupling and ferromagnetic order.

Comparison of harmonic phonon dispersions between NM and FM states
In the main text, we present the phonon dispersions of the non-magnetic state.This is because (i) the ferromagnetic state is not consistent with experiment, and (ii) the dynamical stability properties of lead apatite are largely independent of the precise electronic structure of the system, as discussed in Section 3.1 above.In addition, in this section, we provide both electronic structures and phonon dispersions of both non-magnetic and ferromagnetic states with and without spin-orbit coupling, aiming for a full picture of the system.We find harmonic instabilities in both non-magnetic and ferromagnetic states (see Fig. 7), showing marginal spin-phonon coupling.Supplementary Figure 7: Comparison of electron band structure and phonon dispersion between non-magnetic and ferromagnetic states at various levels of theory for Pb 9 Cu(PO 4 ) 6 O. a Electron band structure and b phonon dispersion of non-magnetic (NM) and ferromagnetic (FM) states obtained using the PBE functional.c Electron band structure and d phonon dispersion of NM and FM states obtained using the PBEsol+U functional with spin-orbit coupling (SOC).In d, the phonon spectra is obtained using a 1 × 1 × 1 q-point grid, which implies that only the phonons at the Γ point are directly calculated.This is sufficient to investigate the Γ point instability which is the primary instability in this system as discussed in Fig. 3 of the main text.

FIG. 2 .
FIG. 2. Electron band structure and phonon dispersion of Pb 9 Cu(PO 4 ) 6 O. a,b Optimized crystal structures of Pb 9 Cu(PO 4 ) 6 O for copper doping at the a Pb(1) and b Pb(2) sites.c,d Electronic band structures of Pb 9 Cu(PO 4 ) 6 O for copper doping at the c Pb(1) and d Pb(2) sites.NM, FM, SOC refer to non-magnetic, ferromagnetic, and spin-orbit coupling, respectively.For doping at the Pb(2) site, the initial non-magnetic configuration converges to a ferromagnetic configuration in the presence of spin-orbit coupling.e Harmonic and anharmonic (300 K) phonon dispersions of the Pb 9 Cu(PO 4 ) 6 O structure with P 3 symmetry for doping at the Pb(1) site.f Harmonic phonon dispersion of the Pb 9 Cu(PO 4 ) 6 O structure with P 1 symmetry for doping at the Pb(2) site.The phonon dispersions for the copper doped cases are obtained using the NM state without SOC (see SI for phonon dispersions of the FM state).The data are obtained with U = 3 eV on the copper 3d orbital.
FIG.3.Γ-distorted P 1 structure and its electron band structure and phonon dispersion.a.Potential energy surface along the imaginary phonon modes at high symmetry points in the Brillouin zone of the P 3 structure for doping at the Pb(1) site (see its harmonic phonon dispersion in Fig.2e).b Crystal structure distorted along the Γ mode of the P 3 structure.c,d.c Harmonic phonon dispersion and d electron band structure of the Γ-distorted P 1 structure.In a-c, the data are obtained using the NM state without SOC.e. Non-magnetic density of states (DOS) and Cu partial DOS of both the P 3 and Γ-distorted P 1 structures.See also partial DOS of other atoms in Supplementary Figure10.The data are obtained with U = 3 eV on the copper 3d orbital.

FIG. 4 .
FIG. 4. Volume-Hubbard U phase diagram.a. Representative harmonic phonon dispersions of the non-magnetic P 3 structure doped at the Pb(1) site for different U values and volumes.b.Volume-Hubbard U phase diagram of Pb 9 Cu(PO 4 ) 6 O for doping at the Pb(1) site.The volume of the P 3 structure is presented on the horizontal axis as a reference, corresponding to a pressure range of 0 − 25 GPa.We find slightly larger volume changes in the Γ-distorted P 1 structure.The color bar indicates the enthalpy difference between the P 3 and P 1 structures (in meV per atom).SOC is not considered in a and b.

Pb 10 (Supplementary Figure 1 :
PO 4 ) 6 O, as illustrated in Fig. 1.For copper-doped lead apatite Pb 9 Cu(PO 4 ) 6 O we find that a k-point grid size of 6 × 6 × 8 leads to converged results.Convergence tests on the electronic structure calculations underpinning the phonon dispersions.a,b.Convergence tests of phonon dispersions with respect to a energy cutoff and b k-point grid for Pb 10 (PO 4 ) 6 O. c The k-point grid convergence test on the phonon dispersion for Pb 9 Cu(PO 4 ) 6 O.The phonon dispersions for the copper doped case are obtained using a NM state without SOC.The calculations are performed using the PBE and PBEsol+U (U = 3 eV) methods for Pb 10 (PO 4 ) 6 O and Pb 9 Cu(PO 4 ) 6 O, respectively.

Supplementary Figure 4 :
Potential energy surface of Pb 10 (PO 4 ) 6 O. a. Harmonic phonon dispersion of Pb 10 (PO 4 ) 6 O calculated using a Farey nonuniform grid of size (2 × 2 × 2) ∪ (3 × 3 × 1) which explicitly includes the K, M, and H points. b.Potential energy surface along the imaginary phonon modes at high symmetry points in the Brillouin zone.The calculations are performed using the PBEsol exchange-correlation functional.

Supplementary Figure 6 :
Harmonic phonon dispersion of Pb 9 Cu(PO 4 ) 6 (OH) 2 at various levels of theory.a. Hubbard U dependence of harmonic phonon band dispersion with U applied to the copper 3d orbital.b.The lowest imaginary Γ-mode energy in the harmonic phonon dispersion as a function of Hubbard U .In a-b, the phonon dispersion and the Γ-mode energy are obtained using a NM state without SOC.c. Spin-orbit coupling and ferromagnetism effects on the lowest imaginary Γ-mode energy.