# Electrostatic Performance of InSb, GaSb, Si and Ge p-channel Nanowires

Celso Martinez-Blanque  Dpto. de Electrónica y Tecnología de Computadores,
Facultad de Ciencias, Universidad de Granada. Av.Fuentenueva S/N, 18071
Granada, Spain  Enrique G. Marin  Dipartimento di Ingegneria
dell´Informazione, Università di Pisa, Via G. Caruso 16, 56122 Pisa, Italy
Alejandro Toral  José M. Gonzalez-Medina  Francisco G. Ruiz<EMAIL_ADDRESS>Andrés Godoy<EMAIL_ADDRESS>Francisco Gámiz  Dpto. de Electrónica y Tecnología
de Computadores, Facultad de Ciencias, Universidad de Granada. Av.Fuentenueva
S/N, 18071 Granada, Spain

###### Abstract

The electrostatic performance of p-type nanowires (NWs) made of InSb and GaSb,
with special focus on their gate capacitance behavior, is analyzed and
compared to that achieved by traditional semiconductors usually employed for
p-MOS such as Si and Ge. To do so, a self-consistent
$\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ simulator has been implemented to
achieve an accurate description of the Valence Band and evaluate the charge
behavior as a function of the applied gate bias. The contribution and role of
the constituent capacitances, namely the insulator, centroid and quantum ones
are assessed. It is demonstrated that the centroid and quantum capacitances
are strongly dependent on the semiconductor material. We find a good inherent
electrostatic performance of GaSb and InSb NWs, comparable to their Ge and Si
counterparts making these III-Sb compounds good candidates for future
technological nodes.

##  I Introduction

In the downscaling roadmap of Metal-Oxide-Semiconductor Field-Effect-
Transistors (MOSFETs), industry has already reached the milestone of
transition from planar to three-dimensional devices enforced by the necessity
of controlling the Short Channel Effects (SCEs). In this context, nanowires
(NWs) have been postulated as the most efficient structure to reduce SCEs as
they provide a higher gate electrostatic control of the channel compared to
other multigate devices. In addition, the requirements of both high
performance and low power devices demand the use of materials different to Si,
able to reduce the supply voltage while maintaining or increasing the on
current ($I_{\mathrm{ON}}$).

In this scenario, and due to their high bulk mobility, III-V materials have
attracted extensive research interest in recent years DelAlamo2011 , being
recognized as promising building blocks for the next generation of electronics
and photonics Riel2014 . Most of current works in the literature focus on
n-type devices, where materials such as GaAs, InAs and InGaAs have already
demonstrated impressive performance. Unfortunately, the enhancement in the
electron mobility of these materials is not accompanied by a similar
enhancement in the hole mobility, resulting in a poorer performance for pMOS
transistors. To achieve fully operational complementary-MOS (CMOS) circuits,
the tolerable asymmetry between both transistors cannot be very marked, and
thus, the search for a high performance nMOS is bound to find its equivalent
pMOS. Several materials are currently being investigated as technologically
feasible p-type channels, being Ge one of the preferred choices due to its
high bulk mobility Pillarisetty2011 . Nonetheless, increasingly more attention
has been focused on antimonide compounds, such as InSb and GaSb, owing to
their excellent bulk hole mobilities among the rest of III-V compounds
Wernersson2013 . Particularly, GaSb has demonstrated the ability to provide
both, high bulk hole and electron mobility,
$\sim$$1400\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
and
$7000\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$,
respectively Shimoida2013 .

Figure 1: Geometry of the cylindrical NW, where Tins is the gate insulator
thickness and R the semiconductor radius.

A high mobility, however, does not guarantee a good performance. A potential
substitute of silicon must also guarantee a large charge modulation in the
channel, namely a good electrostatic control of it. This characteristic of a
device is mainly determined by the gate capacitance and the inversion charge.
In this sense, III-V materials have been largely accused to be affected by a
Density of States (DOS) bottleneck. However, there are few studies dealing
with this issue and its impact on the overall performance of pMOS NWs.

To this purpose, in this work we address the electrostatic properties and gate
capacitance of GaSb and InSb NWs, and compare them to that achieved by Si and
Ge devices, assessing whether the expected high mobility of these compound
semiconductors may be degraded by a poorer electrostatic performance.

The outline of the paper is as follows. First, Section II describes the
numerical simulator developed for this study and the main parameters employed
to model each material. Next, in Section III, the main results regarding the
bandstructure, charge distribution and gate capacitance for the different NWs
are shown and analyzed, focusing on the main differences between the four
considered materials. That Section also provides insights on the influence of
the insulator capacitance, which is modified by changing the dielectric
material. Finally, the main conclusions of this work are drawn in Sec. IV.

##  II Numerical Simulator

| Si | Ge | GaSb | InSb
---|---|---|---|---
$E_{\mathrm{g}}\;$($\mathrm{eV}$) | $1.17$ | $0.742$ | $0.812$ | $0.235$
$\Delta_{\mathrm{SO}}\;$($\mathrm{eV}$) | $0.044$ | $0.286$ | $0.76$ | $0.81$
$E_{\mathrm{P}}\;$($\mathrm{eV}$) | $-$ | $-$ | $18.0$ | $23.3$
$m_{\mathrm{c}}\;$ | $-$ | $-$ | $0.039$ | $0.0135$
$\gamma_{1}$ | $3.55$ | $13.27$ | $13.4$ | $34.8$
$\gamma_{2}$ | $0.65$ | $4.32$ | $4.7$ | $15.5$
$\gamma_{3}$ | $1.26$ | $5.61$ | $6.0$ | $16.5$
$\varepsilon_{\mathrm{r}}$ | $11.7$ | $16.2$ | $15.7$ | $16.8$
$\Xi_{0}\;$ ($\mathrm{eV}$) | $4.05$ | $4.00$ | $4.06$ | $4.59$
Table 1: $\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ parameters for Si Shin2009
, Ge Fischetti1996 , GaSb and InSb Vurgaftman2001 .

In this study we are interested in the electrostatic behavior of cylindrical
gate-all-around (GAA) p-type NWs made of four different materials: two III-V
antimonides, GaSb and InSb, and two group IV semiconductors, Si and Ge.
Particularly, we have studied NWs with 5nm diameter oriented along the [111]
direction (Fig. 1), since this orientation exhibits the best transport
properties for holes in Si and Ge Neophytou2013c ; Niquet2012d ; Martinez-
Blanque2014 . To allow a fair comparison of the electrostatics, the same
high-$\kappa$ insulator, $\mathrm{Al}_{2}\mathrm{O}_{3}$, with thickness
$T_{\mathrm{ins}}=1.5$nm and relative dielectric constant
$\varepsilon_{\mathrm{r}}=9$, is considered in all cases. This way, the
contribution of the insulator capacitance ($C_{\mathrm{ins}}$) is identical
for all the NWs and the comparison can be focused on the channel material.

In order to describe the electrical performance of the semiconductor NWs under
study, we have used a self-consistent solver for the two-dimensional (2D)
Schrödinger and Poisson equations in the transversal cross-section of the
cylindrical GAA NWs. The effective mass approximation with non-parabolic
corrections has demonstrated to provide a simple but accurate description of
the band structure of III-V compounds close to the Conduction Band (CB)
minimum and has, therefore, been used in similar studies of n-type NWs
Marin2015 . However, this method fails in the description of the more complex
Valence Band (VB) of group IV (Si, Ge) and III-V semiconductors due to the
coupling of the different subbands encompassing the VB: Heavy Holes (HH),
Light Holes (LH) and Split-Off (SO). Thus, to achieve an accurate description
of the VB of such materials, it is mandatory to use a model which accounts for
the coupling between these subbands.

In this work, the $\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ method, which
explicitly accounts for the coupling between VB subbands has been implemented
to calculate the VB structure and the associated envelope wave functions
Lassen2006 ; Shin2010 ; RaviKishore2012 . The
$\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ simulation scheme employs a reduced
set of parameters which can be obtained semi-empirically as they are tuned
with other atomistic descriptions. The results presented in this work have
been benchmarked against atomistic methods, such as Tight-Binding (TB),
providing a good agreement around the band edges Shin2009 ; Shin2010 . For
indirect gap materials, such as Si and Ge, a six-band model is accurate enough
since the different subbands forming the VB are considered explicitly, whereas
the rest of bands, including the CB, are regarded as remote. However, for
direct materials with a small gap, such as the III-Sb compounds here analyzed,
the proximity of the CB to the VB has a non negligible influence on the
resulting VB structure. The CB can also be explicitly accounted for in an
eight-band $\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ model Bahder1990 . The
implemented model uses a set of semi-empirical parameters directly related to
the different anisotropic effective masses of the involved subbands, the so-
called Luttinger parameters, that can be found in the literature
Vurgaftman2001 . These parameters are valid for both, six-band and eight-band
$\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ models. However, for the latter,
they have to be reduced to exclude the effect of the CB, which is explicitly
taken into account in the model Bahder1990 . The
$\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ parameters used in this study, and
the references from where they have been extracted, are presented in Table 1.
The Kane coupling parameter $E_{\mathrm{P}}$ has been reduced for GaSb, as
suggested in Veprek2009 , to avoid the appearance of spurious solutions.

##  III Results and Discussion

###  III.1 Valence Band bandstructure

Figure 2: Bandstructure for 5nm diameter cylindrical NWs with channel
materials: Si (upper left), Ge (upper right), GaSb (bottom left) and InSb
(bottom right), at a gate overdrive voltage of 0V (left) and
$-0.4\text{\,}\mathrm{V}$ (right). The energy is referred to the VB edge, and
the Fermi level is depicted as a dashed line. Figure 3: DOS evaluated for 5nm
diameter cylindrical NWs with four different channel materials: Si (upper
left), Ge (upper right), GaSb (bottom left) and InSb (bottom right), at a gate
overdrive voltage of 0V (left) and $-0.4\text{\,}\mathrm{V}$ (right). The
energy is referred to the VB edge, and the Fermi level is depicted as a dashed
line.

Figure 2 depicts the VB of the four NWs with the different materials: from top
to bottom and left to right, Si, Ge, GaSb and InSb, respectively. Two
different gate overdrive voltages, $V_{\mathrm{G}}-V_{\mathrm{T}}$ (with
$V_{\mathrm{T}}$ the threshold voltage), are depicted for each of these
devices, 0V (left) and $-0.4\text{\,}\mathrm{V}$ (right). Here,
$V_{\mathrm{T}}$ has been evaluated as the voltage that maximizes
$\,\text{d}^{2}Q_{\mathrm{i}}/\,\text{d}V_{\mathrm{G}}^{2}$ Marin2013 , where
$Q_{\mathrm{i}}$is the semiconductor charge. The depicted bandstructures
energies are referred to the maximum of the VB, and the Fermi energy is
represented as an horizontal dashed line. Additionally, Fig. 3 plots the DOS
as a function of the energy for each of the devices and biases considered in
Fig. 2.

The VB structure calculated with the $\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$
method has been re-ordered using an algorithm that minimizes the changes of
the subband group velocity Zlatan2012 . As a result, although the VB structure
in the NW is the result of a mix between the bulk hole subbands (HH, LH and
SO), Fig. 2 shows that we can still sort the resulting them into HH-dominated
(HHD) subbands and LH-dominated (LHD) subbands, exhibiting two trends
corresponding to different curvatures. HHD subbands have lower curvature and
group velocity, and higher DOS than the LHD ones.

At a first glance, the Si device presents the highest DOS due to: i) the
larger number of subbands per unit energy, and ii) their lower curvature, both
for HHD and LHD subbands. The Ge, InSb and GaSb NWs exhibit similar
bandstructures, where the HH band is much less influential, resulting in just
one HHD subband and a lower DOS. When increasing the gate overdrive voltage
from $0\text{\,}\mathrm{V}$ to $-0.4\text{\,}\mathrm{V}$, a reduction in the
energy distance between the HHD and the LHD subbands, and also between the
first and second LHD subbands is observed for all materials (Fig. 2) resulting
in an increased DOS (Fig. 3).

###  III.2 Hole distribution and centroid

Employing the bandstructure and wave functions obtained from the
$\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ simulations, the charge distribution
can be evaluated. Our results show an almost isotropic behavior of the charge
for the four devices under study, due to the [111] orientation considered.
Fig. 4 depicts the radial hole distribution $p(r)$ in the NWs made of
different materials (dotted line for Si, solid line for Ge, dash-dotted for
GaSb and dashed for InSb) when two different gate overdrive voltages are
considered: 0V (a) and $-0.4\text{\,}\mathrm{V}$ (b).

In the subthreshold and near-threshold regimes, Fig. 4(a), the charge for all
the devices presents a quite similar trend, with the maximum located at the
device center. Nevertheless, when the applied gate voltage increases, Fig.
4(b), the charge distribution is shifted towards the semiconductor-insulator
interface, presenting a clear peak far from the device center for Si and Ge,
whereas for GaSb and InSb it exhibits a broader flat hole distribution at the
center of the NW.

Figure 4: Hole density along the radial axis for 5nm cylindrical NWs made of
different materials (Si dotted blue, Ge solid red, GaSb dash-dotted orange and
InSb dashed green) at a gate overdrive voltage of 0V (a) and
$-0.4\text{\,}\mathrm{V}$ (b). Vertical lines and symbols denote the position
of the centroid $\Delta$ for each material: circle for Si, triangle for Ge,
star for GaSb and square for InSb. $r=0$ corresponds to the center of the NW.

The displacement of the charge towards the semiconductor-insulator interface
can be quantitatively evaluated using the charge centroid Lopez-Villanueva2000
; Lee2012 . We will consider the logarithmic weighting of the hole
distribution proposed in Lee2012 , which is closely related to the gate
capacitance analyzed in the next Section. So that, the centroid ($\Delta$) of
the hole distribution, referred to the semiconductor-insulator interface, is
calculated as:

$\Delta=R\left(1-\exp\left(-2\pi x_{\mathrm{i}}\right)\right)$ (1)

with $x_{\mathrm{i}}$ defined as:

$x_{\mathrm{i}}=\frac{1}{2\pi}\frac{\int_{0}^{R}r\,\log(\frac{R}{r})\,p(r)\,\text{d}r}{\int_{0}^{R}r\,p(r)\,\text{d}r}\,\text{.}$
(2)

Figure 4 indicates the position of the centroid with a vertical line and a
symbol laying on their respective curves. At low charge densities
($V_{\mathrm{G}}-V_{\mathrm{T}}=$0\text{\,}\mathrm{V}$$), carriers are mainly
located close to the center of the device in all the cases and very similar
values are found for the centroid. For large overdrive voltages
($V_{\mathrm{G}}-V_{\mathrm{T}}=$-0.4\text{\,}\mathrm{V}$$), the charge
centroid moves closer to the interface, more markedly for Si than for the
other materials.

Figure 5: Variation of the centroid position with respect to the
semiconductor-insulator interface for the same devices as in Fig. 4.

Figure 5 provides deeper insight in the dependence of the centroid on
$V_{\mathrm{G}}-V_{\mathrm{T}}$ for the four devices under consideration. As
can be seen, the Si device has the lower value for the whole range of gate
overdrive voltages analyzed. In the subthreshold regime, the centroid remains
constant showing similar results regardless the material considered
consistently with the charge distribution presented in Fig. 4. When the device
gets into inversion, the charge centroid decreases, (as the charge is shifted
towards the interface) with different trends for each material. Si shows the
largest modulation of the centroid with $V_{\mathrm{G}}$, whilst InSb shows
the smallest. With the exception of Si, a noteworthy reduction on the centroid
control by the gate is observed when
$V_{\mathrm{G}}-V_{\mathrm{T}}<$-0.2\text{\,}\mathrm{V}$$. As it will be shown
next, this slope reduction has a direct effect on the gate capacitance of
those devices.

###  III.3 Gate Capacitance

The best figure to assess the gate control over the channel charge is the gate
capacitance, $C_{\mathrm{G}}$. With the reduction of the channel length, it
has become challenging to keep high values for $C_{\mathrm{G}}$, being this
one of the main reasons to employ GAA geometries. Thus, a comparison of this
magnitude for the four devices under consideration will provide a good measure
of their electrostatic performance.

For the cylindrical GAA under study, and assuming an isotropic distribution of
the charge density in the cross-section, the Gauss’s Law is straightforwardly
applicable. Thus, $C_{\mathrm{G}}$ can be regarded as the contribution of
three series capacitances:

$\frac{1}{C_{\mathrm{G}}}=\frac{1}{C_{\mathrm{ins}}}+\frac{1}{C_{\mathrm{c}}}+\frac{1}{C_{\mathrm{q}}}$
(3)

As aforementioned, the same dielectric (Al2O3) has been considered throughout
this work, ($T_{\mathrm{ins}}=$1.5\text{\,}\mathrm{nm}$$ and
$\varepsilon_{\text{ins}}=9\varepsilon_{0}$). So that, $C_{\mathrm{ins}}$ can
be evaluated as:

$C_{\mathrm{ins}}=\frac{2\pi\varepsilon_{\text{ins}}}{\ln{\left(1+\frac{T_{\mathrm{ins}}}{R}\right)}}$
(4)

where $\varepsilon_{\text{ins}}$ is the insulator dielectric constant, and the
rest of the parameters have been previously defined. $C_{\mathrm{c}}$ and
$C_{\mathrm{q}}$ are the centroid and quantum capacitance terms, respectively,
that can be evaluated as:

$C_{\mathrm{c}}=\frac{\partial
Q_{\mathrm{i}}}{\left(\partial\psi_{\mathrm{s}}-\psi_{\mathrm{c}}\right)}$ (5)

and

$C_{\mathrm{q}}=\frac{\partial Q_{\mathrm{i}}}{\partial\psi_{\mathrm{c}}}$ (6)

where $\psi_{\mathrm{c}}$ and $\psi_{\mathrm{s}}$ correspond to the potential
evaluated at the device center and the semiconductor-insulator interface,
respectively. $C_{\mathrm{c}}$ and $C_{\mathrm{q}}$ can be combined into the
so-called inversion capacitance $C_{\mathrm{inv}}=\partial
Q_{\mathrm{i}}/\partial\psi_{\mathrm{s}}$. $C_{\mathrm{q}}$, has been related
to the finite DOS resulting from quantization Luryi1988 ; Marin2013 . As for
$C_{\mathrm{c}}$, it can be calculated as Ruiz2010 ; Lee2012 :

$\frac{1}{C_{\mathrm{c}}}=\frac{x_{\mathrm{i}}}{\varepsilon_{\mathrm{s}}}+\frac{Q_{\mathrm{i}}}{\varepsilon_{\mathrm{s}}}\frac{\,\text{d}x_{\mathrm{i}}}{\,\text{d}Q_{\mathrm{i}}}$
(7)

where $\varepsilon_{\mathrm{s}}$ is the semiconductor dielectric constant and
$x_{\mathrm{i}}$ was defined in (2), remarking the close relation between
$C_{\mathrm{c}}$ and the charge centroid Lee2012 .

Figure 6: Gate capacitance $C_{\mathrm{G}}$ for 5nm cylindrical NWs with
different materials (Si in blue dotted line, Ge in red solid line, GaSb in
orange dash-dotted line and InSb in green dashed line) as a function of the
gate overdrive voltage $V_{\mathrm{G}}-V_{\mathrm{T}}$. The constituents
capacitances $C_{\mathrm{ins}}$, $C_{\mathrm{c}}$ and $C_{\mathrm{q}}$ are
also depicted.

In Fig. 6, the gate capacitances and their components are depicted for the
four materials considered here. We observe a similar trend on the total gate
capacitance for all of them. The Si device provides the highest value of
$C_{\mathrm{G}}$ in the whole range of gate voltages but the rest of materials
show a comparable performance in terms of $C_{\mathrm{G}}$with a maximum
variation of 15% in strong inversion. For a better understanding of the
behavior of each device, the different constituent capacitances must be
analyzed. As aforementioned, $C_{\mathrm{ins}}$ is the same for the four
devices, so we will restrict our analysis to $C_{\mathrm{q}}$ and
$C_{\mathrm{c}}$. Due to the extremely low values of $C_{\mathrm{q}}$ in the
subthreshold regime, the influence of $C_{\mathrm{ins}}$ and $C_{\mathrm{c}}$
in the total capacitance can be neglected in this range. For gate voltages
close to $V_{\mathrm{T}}$ and slightly higher, $C_{\mathrm{G}}$ is influenced
by both $C_{\mathrm{c}}$ and $C_{\mathrm{q}}$. In Si, the faster enhancement
of $C_{\mathrm{q}}$, due to the larger DOS, is spoiled by a low
$C_{\mathrm{c}}$. For Ge, InSb and GaSb, $C_{\mathrm{q}}$ increases more
softly, remaining below $20\text{\,}\mathrm{pF}\text{\,}{\mathrm{cm}}^{-1}$ up
to $V_{\mathrm{G}}-V_{\mathrm{T}}$ around $-0.2\text{\,}\mathrm{V}$. This DOS
bottleneck is, however, compensated by the higher $C_{\mathrm{c}}$ with
respect to Si, in weak inversion. The behavior of $C_{\mathrm{c}}$ is
controlled by the dielectric constant of the semiconductor and the centroid.
For Si, the lower semiconductor dielectric constant, makes the
$C_{\mathrm{c}}$ term lower than for Ge, InSb and GaSb for small gate
overdrives. Nevertheless, as the gate voltage becomes more negative, this
effect is compensated by the smaller value of the centroid (see Fig. 5). The
centroid warping in Fig. 5 also explains the degradation of $C_{\mathrm{c}}$
in Ge, InSb and GaSb. In spite of this degradation, GaSb provides similar
$C_{\mathrm{c}}$ values as compared to Si for
$V_{\mathrm{G}}-V_{\mathrm{T}}\leq-0.4$V.

From Fig. 6 it can be observed that the total gate capacitance is limited by
the low value of $C_{\mathrm{ins}}$, which veils the effect of the rest of
capacitances for large overdrive voltages. This capacitance does not depend on
the channel material, and can be raised by using higher $\kappa$ materials as
gate insulators. New simulations have thus been carried out with HfO2, a
high-$\kappa$ dielectric with $\varepsilon_{\text{ins}}=25\varepsilon_{0}$ and
$T_{\mathrm{ins}}=$1.5\text{\,}\mathrm{nm}$$. Fig. 7 compares $C_{\mathrm{G}}$
as a function of the overdrive voltage for the four semiconductors and the two
insulators, Al2O3 and HfO2. As can be seen, $C_{\mathrm{G}}$ values achieved
with HfO2 double those attained with Al2O3. Again, the Si device presents the
best performance in terms of gate capacitance, but GaSb holds the comparison
showing an excellent electrostatics behavior in the whole range of applied
bias. On the other hand, Ge and InSb NWs show similar trends with a slight
reduction as $V_{\mathrm{G}}-V_{\mathrm{T}}$ becomes more negative. Hence,
these results demonstrate the excellent electrostatic performance, in terms of
gate capacitance, of the p-channel antimonide NWs considered in this study. In
particular, GaSb is the most attractive alternative as it provides
$C_{\mathrm{G}}$ values similar to Si in addition to the already demonstrated
excellent transport properties.

Figure 7: Gate capacitance $C_{\mathrm{G}}$ for $5\text{\,}\mathrm{nm}$
cylindrical NWs with two different high-$\kappa$ insulators, Al2O3 and HfO2,
with a similar thickness of $1.5\text{\,}\mathrm{nm}$. Si in blue dotted line,
Ge in red solid line, GaSb in orange dash-dotted line and InSb in green dashed
line as a function of the gate overdrive voltage
$V_{\mathrm{G}}-V_{\mathrm{T}}$.

##  IV Conclusion

A comprehensive study of the electrostatic performance of 5nm diameter NWs has
been carried out focusing on p-type channels made of four different materials:
InSb, GaSb, Ge and Si. To do so, we have developed a numerical simulator that
self-consistently solves the Poisson equation and the eight-band
$\boldsymbol{k}\\!\cdot\\!\boldsymbol{p}$ method in the cross section of
cylindrical NWs. Our results indicate that, in spite of their lower effective
mass and smaller density of states, GaSb and InSb p-type NWs hold the
comparison with Si and, in the case of GaSb, outperforms Ge in terms of gate
capacitance and inversion charge. Thus, their good electrostatic performance
combined with the expectation of superior transport characteristics place GaSb
and InSb as attractive alternatives for p-type CMOS logic based on NWs.

## References

  * (1) J. A. del Alamo, “Nanometre-scale electronics with III–V compound semiconductors,” _Nature_ , vol. 479, no. 7373, pp. 317–323, 2011.
  * (2) H. Riel, L. E. Wernersson, M. Hong, and J. A. del Alamo, “III-V compound semiconductor transistors - from planar to nanowire structures,” _MRS Bulletin_ , vol. 39, no. 08, pp. 668–677, 2014.
  * (3) R. Pillarisetty, “Academic and industry research progress in germanium nanodevices,” _Nature_ , vol. 479, pp. 324–328, 2011.
  * (4) M. Borg and L. E. Wernersson, “Synthesis and properties of antimonide nanowires,” _Nanotechnology_ , vol. 24, no. 20, p. 202001, 2013.
  * (5) K. Shimoida, Y. Yamada, H. Tsuchiya, and M. Ogawa, “Orientational Dependence in Device Performances of InAs and Si Nanowire MOSFETs Under Ballistic Transport,” _IEEE Transactions on Electron Devices_ , vol. 60, no. 1, pp. 117–122, Jan. 2013.
  * (6) M. Shin, “Full-quantum simulation of hole transport and band-to-band tunneling in nanowires using the k$\cdot$p method,” _Journal of Applied Physics_ , vol. 106, no. 5, p. 054505, 2009.
  * (7) M. V. Fischetti and S. E. Laux, “Band structure, deformation potentials, and carrier mobility in strained Si, Ge, and SiGe alloys,” _Journal of Applied Physics_ , vol. 80, no. 4, p. 2234, 1996.
  * (8) I. Vurgaftman, J. R. Meyer, and L. R. Ram-Mohan, “Band parameters for III-V compound semiconductors and their alloys,” _Journal of Applied Physics_ , vol. 89, no. 11, p. 5815, 2001.
  * (9) N. Neophytou, O. Baumgartner, Z. Stanojević, and H. Kosina, “Bandstructure and mobility variations in p-type silicon nanowires under electrostatic gate field,” _Solid-State Electronics_ , vol. 90, pp. 44–50, Dec. 2013.
  * (10) Y. M. Niquet and C. Delerue, “Carrier mobility in strained Ge nanowires,” _Journal of Applied Physics_ , vol. 112, no. 8, p. 084301, 2012.
  * (11) C. Martinez-Blanque, F. G. Ruiz, A. Godoy, E. G. Marin, L. Donetti, and F. Gámiz, “Influence of alloy disorder scattering on the hole mobility of SiGe nanowires,” _Journal of Applied Physics_ , vol. 116, no. 24, p. 244502, Dec. 2014.
  * (12) E. G. Marin, F. G. Ruiz, A. Godoy, I. M. Tienda-Luna, and F. Gámiz, “Mobility and capacitance comparison in scaled ingaas versus si trigate mosfets,” _IEEE Electron Device Letters_ , vol. 36, no. 2, pp. 114–116, Feb 2015.
  * (13) B. Lassen, M. Willatzen, R. Melnik, and L. C. Lew Yan Voon, “Electronic structure of free-standing InP and InAs nanowires,” _Journal of Materials Research_ , vol. 21, no. 11, pp. 2927–2935, 2006.
  * (14) M. Shin, S. Lee, and G. Klimeck, “Computational Study on the Performance of Si Nanowire pMOSFETs Based on the k$\cdot$p Method,” _IEEE Transactions on Electron Devices_ , vol. 57, no. 9, pp. 2274–2283, Sep. 2010.
  * (15) V. V. Ravi Kishore, N. Čukarić, B. Partoens, M. Tadić, and F. M. Peeters, “Hole subbands in freestanding nanowires: six-band versus eight-band k$\cdot$p modelling.” _Journal of physics. Condensed matter : an Institute of Physics journal_ , vol. 24, no. 13, p. 135302, Apr. 2012.
  * (16) T. B. Bahder, “Eight-band k$\cdot$p model of strained zinc-blende crystals,” _Physical Review B_ , vol. 41, no. 17, pp. 11 992–12 001, 1990\.
  * (17) R. G. Veprek, S. Steiger, and B. Witzigmann, “Operator ordering, ellipticity and spurious solutions in k$\cdot$p calculations of III-nitride nanostructures,” _Optical and Quantum Electronics_ , vol. 40, no. 14-15, pp. 1169–1174, Jan. 2009.
  * (18) E. G. Marin, F. G. Ruiz, I. M. Tienda-Luna, A. Godoy, and F. Gámiz, “Analytical Gate Capacitance Modeling of III-V Nanowire Transistors,” _IEEE Transactions on Electron Devices_ , vol. 60, no. 5, pp. 1590–1599, May 2013.
  * (19) Z. Stanojević, V. Sverdlov, O. Baumgartner, and H. Kosina, “Subband engineering in n-type silicon nanowires using strain and confinement,” _Solid-State Electronics_ , vol. 70, pp. 73 – 80, 2012, selected Full-Length Papers from the {EUROSOI} 2011 Conference.
  * (20) J. A. Lopez-Villanueva, P. Cartujo-Cassinello, F. Gámiz, J. Banqueri, and A. J. Palma, “Effects of the inversion-layer centroid on the performance of double-gate MOSFET’s,” _IEEE Transactions on Electron Devices_ , vol. 47, no. 1, pp. 141–146, 2000.
  * (21) Y. Lee, K. Kakushima, K. Natori, and H. Iwai, “Gate capacitance modeling and diameter-dependent performance of nanowire MOSFETs,” _IEEE Transactions on Electron Devices_ , vol. 59, no. 4, pp. 1037–1045, 2012.
  * (22) S. Luryi, “Quantum capacitance devices,” _Applied Physics Letters_ , vol. 52, pp. 501–503, 1988.
  * (23) F. G. Ruiz, I. M. Tienda-Luna, A. Godoy, L. Donetti, and F. Gámiz, “A Model of the Gate Capacitance of Surrounding Gate Transistors: Comparison With Double-Gate MOSFETs,” _IEEE Transactions on Electron Devices_ , vol. 57, no. 10, pp. 2477–2483, Oct. 2010.


# Hypercontractivity of the semigroup
of the fractional laplacian on the $n$-sphere

Rupert L. Frank  Mathematics 253-37, Caltech, Pasadena, CA 91125, USA, and
Mathematisches Institut, Ludwig-Maximilans Universät München, Theresienstr.
39, 80333 München, Germany<EMAIL_ADDRESS>and  Paata Ivanisvili
Department of Mathematics, North Carolina State University, Raleigh, NC 27695
<EMAIL_ADDRESS>
###### Abstract.

For $1<p\leq q$ we show that the Poisson semigroup $e^{-t\sqrt{-\Delta}}$ on
the $n$-sphere is hypercontractive from $L^{p}$ to $L^{q}$ in dimensions
$n\leq 3$ if and only if $e^{-t\sqrt{n}}\leq\sqrt{\frac{p-1}{q-1}}$. We also
show that the equivalence fails in large dimensions.

###### Key words and phrases:

Hypercontractivity, Poisson Semigroup, n-sphere

###### 2010 Mathematics Subject Classification:

39B62, 42B35, 47A30

${}$${}$footnotetext: © 2021 by the author. This paper may be reproduced, in
its entirety, for non-commercial purposes.

##  1\. Introduction

###  1.1. Poisson semigroup on the sphere

Let

$\displaystyle\mathbb{S}^{n}=\\{x\in\mathbb{R}^{n+1}\,:\,\|x\|=1\\}$

be the unit sphere in $\mathbb{R}^{n+1}$, where
$\|x\|=\sqrt{x_{1}^{2}+\ldots+x_{n+1}^{2}}$ for
$x=(x_{1},\ldots,x_{n+1})\in\mathbb{R}^{n+1}$. Let $\Delta$ be the
Laplace–Beltrami operator on $\mathbb{S}^{n}$. We will be working with
spherical polynomials $f:\mathbb{S}^{n}\to\mathbb{C}$, i.e., finite sums

$f(\xi)=\sum_{d\geq 0}H_{d}(\xi),$

where $H_{d}$ satisfies

$\Delta H_{d}=-d(d+n-1)H_{d}.$

The heat semigroup $e^{t\Delta}$ is defined by $e^{t\Delta}f=\sum_{d\geq
0}e^{-d(d+n-1)t}H_{d}$. The hypercontractivity result for the heat semigroup
on $\mathbb{S}^{n}$ states that for any $1\leq p\leq q<\infty$, any integer
$n\geq 1$, and any $t\geq 0$ we have

(1) $\displaystyle\|e^{t\Delta}f\|_{q}\leq\|f\|_{p}\quad\text{for all}\
f\qquad\text{if and only if}\qquad e^{-tn}\leq\sqrt{\frac{p-1}{q-1}},$

where
$\|f\|_{p}^{p}=\|f\|_{L^{p}(\mathbb{S}^{n},d\sigma_{n})}^{p}=\int_{\mathbb{S}^{n}}|f|^{p}d\sigma_{n}$,
and $d\sigma_{n}$ is the normalized surface area measure of $\mathbb{S}^{n}$.
The case $n=1$ was solved independently in [9] and [10], and the general case
$n\geq 2$ was settled in [7]. We remark that the condition
$e^{-tn}\leq\sqrt{\frac{p-1}{q-1}}$ in (1) is different from the classical
hypercontractivity condition $e^{-t}\leq\sqrt{\frac{p-1}{q-1}}$ in Gauss space
due to Nelson [8], and on the hypercube due to Bonami [2]. The appearance of
the extra factor $n$ in (1) can be explained from the fact that the spectral
gap (the smallest nonzero eigenvalue) of $-\Delta$ equals $n$.

In [7] the authors ask what the corresponding hypercontractivity estimates are
for the Poisson semigroup on $\mathbb{S}^{n}$. As pointed out in [7], there
are two natural Poisson semigroups on $\mathbb{S}^{n}$ one can consider: 1)
$e^{-t\sqrt{-\Delta}}f$, and 2) $P_{r}f=\sum r^{d}H_{d}$, $r\in[0,1]$. Notice
that when $n=1$ both of these semigroups coincide (with $r=e^{-t}$). It was
conjectured by E. Stein that

$\|P_{r}f\|_{q}\leq\|f\|_{p}\quad\text{if and only if}\quad
r\leq\sqrt{\frac{p-1}{q-1}}$

holds on $\mathbb{S}^{n}$ for all $n\geq 1$. Besides the case $n=1$ mentioned
above, the case $n=2$ was confirmed in [4], and the general case $n\geq 2$ in
[1].

The question of hypercontractivity for the semigroup $e^{-t\sqrt{-\Delta}}$ on
$\mathbb{S}^{n}$ for $n\geq 2$, however, has remained open. Since the spectral
gap of $\sqrt{-\Delta}$ equals $\sqrt{n}$, it is easy to see that a necessary
condition for the estimate $\|e^{-t\sqrt{-\Delta}}f\|_{q}\leq\|f\|_{p}$ is
$e^{-t\sqrt{n}}\leq\sqrt{\frac{p-1}{q-1}}$; see Section 2.1. One might
conjecture that this necessary condition is also sufficient. Surprisingly, it
turns out the answer is positive in small dimensions and negative in large
dimensions.

######  Theorem 1.1.

Let $1<p<q$, $n\geq 1$, and $t\geq 0$. Then

(2)
$\displaystyle\textup{(i)}\;\;\|e^{-t\sqrt{-\Delta}}f\|_{q}\leq\|f\|_{p}\quad\text{for
all}\
f\qquad\text{implies}\qquad\textup{(ii)}\;\;e^{-t\sqrt{n}}\leq\sqrt{\frac{p-1}{q-1}}.$

Moreover, $(ii)$ implies $(i)$ in dimensions $n\leq 3$. Finally, for any
$q>\max\\{2,p\\}$, there exists $n_{0}=n_{0}(p,q)\geq 4$ such that (ii) does
not imply (i) in dimensions $n$ with $n\geq n_{0}$.

It remains an open problem to find a necessary and sufficient condition on
$t>0$ in dimensions $n\geq 4$ for which the semigroup $e^{-t\sqrt{-\Delta}}$
is hypercontractive from $L^{p}(\mathbb{S}^{n})$ to $L^{q}(\mathbb{S}^{n})$.

##  2\. Proof of Theorem 1.1

###  2.1. The necessity part $\textup{(i)}\Rightarrow\textup{(ii)}$

We recall this standard argument for the sake of completeness. Let
$f(\xi)=1+\varepsilon H_{1}(\xi)$ where $H_{1}$ is any (real) spherical
harmonic of degree $1$, i.e., $\Delta H_{1}=-nH_{1}$. Then
$e^{-t\sqrt{-\Delta}}f(\xi)=1+\varepsilon e^{-t\sqrt{n}}H_{1}(\xi)$. As
$\varepsilon\to 0$, we obtain

$\displaystyle\int_{\mathbb{S}^{n}}|1+\varepsilon
e^{-t\sqrt{n}}H_{1}(\xi)|^{q}d\sigma_{n}$
$\displaystyle=\int_{\mathbb{S}^{n}}\left(1+q\varepsilon
e^{-t\sqrt{n}}H_{1}(\xi)+\frac{q(q-1)}{2}\varepsilon^{2}e^{-2t\sqrt{n}}H^{2}_{1}(\xi)+O(\varepsilon^{3})\right)d\sigma_{n}$
$\displaystyle=1+\frac{q(q-1)}{2}\varepsilon^{2}e^{-2t\sqrt{n}}\|H_{1}\|_{2}^{2}+O(\varepsilon^{3}).$

Thus,

(3)
$\displaystyle\|e^{-t\sqrt{-\Delta}}f\|_{q}=1+\frac{q-1}{2}\varepsilon^{2}e^{-2t\sqrt{n}}\|H_{1}\|_{2}^{2}+O(\varepsilon^{3}).$

Similarly, we have

(4)
$\displaystyle\|f\|_{p}=1+\frac{p-1}{2}\varepsilon^{2}\|H_{1}\|_{2}^{2}+O(\varepsilon^{2}).$

Substituting (3) and (4) into the inequality
$\|e^{-t\sqrt{-\Delta}}f\|_{q}\leq\|f\|_{p}$, and taking $\varepsilon\to 0$ we
obtain the necessary condition $e^{-2t\sqrt{n}}\leq\frac{p-1}{q-1}$ which
coincides with (ii) in (2).

###  2.2. The sufficiency part $\textup{(ii)}\Rightarrow\textup{(i)}$ in
dimensions $n=1,2,3$.

Our goal is to show that if $1<p<q$ and if $t\geq 0$ is such that
$e^{-t2\sqrt{n}}\leq\frac{p-1}{q-1}$, then

(5) $\displaystyle\|e^{-t\sqrt{-\Delta}}f\|_{q}\leq\|f\|_{p}\quad\text{in
dimensions}\quad n=1,2,3.$

The case $n=1$ was confirmed in [10]. In what follows we assume
$n\in\\{2,3\\}$. First we need the fact that the heat semigroup $e^{t\Delta}$
has a nonnegative kernel. Indeed, for each $t>0$ there exists
$K_{t}:[-1,1]\to[0,\infty)$ such that

$e^{t\Delta}f(\xi)=\int_{\mathbb{S}^{n}}K_{t}(\xi\cdot\eta)f(\eta)d\sigma_{n}(\eta),$

where $\xi\cdot\eta=\sum_{j=1}^{n+1}\xi_{j}\eta_{j}$ for
$\xi=(\xi_{1},\ldots,\xi_{n+1})$ and $\eta=(\eta_{1},\ldots,\eta_{n+1})$, see,
for example, Proposition 4.1 in [7]. Next, we recall the subordination formula

(6) $\displaystyle
e^{-x}=\frac{1}{\sqrt{\pi}}\int_{0}^{\infty}e^{-y-x^{2}/(4y)}\frac{dy}{\sqrt{y}}\quad\text{valid
for all}\ x\geq 0,$

By the functional calculus, we deduce that the Poisson semigroup
$e^{-t\sqrt{-\Delta}}$ has a positive kernel with total mass $1$. The latter
fact together with the convexity of the map $x\mapsto|x|^{p}$ for $p\geq 1$
implies that $\|e^{-t\sqrt{-\Delta}}\|_{p}\leq\|f\|_{p}$ for all $t\geq 0$.
Thus, it suffices to verify (5) for those $t\geq 0$ for which
$e^{-2t\sqrt{n}}=\frac{p-1}{q-1}$.

Next we claim that it suffices to verify (5) only for the powers $p,q$ such
that $2\leq p\leq q$. Indeed, assume (5) holds for $2\leq p\leq q$. By duality
and the symmetry of the semigroup $e^{-t\sqrt{-\Delta}}$ we obtain
$\|e^{-t\sqrt{-\Delta}}f\|_{p^{\prime}}\leq\|f\|_{q^{\prime}}$ where
$p^{\prime}=\frac{p}{p-1}$, $q^{\prime}=\frac{q}{q-1}$, $1<q^{\prime}\leq
p^{\prime}\leq 2$. Notice that
$\frac{p-1}{q-1}=\frac{q^{\prime}-1}{p^{\prime}-1}$, thus we extend (5) to all
$p,q$ such that $1<p\leq q\leq 2$. It remains to extend (5) for those powers
$p,q$ when $p\leq 2\leq q$. To do so, let $p\leq 2\leq q$, and let $t\geq 0$
be such $e^{-2t\sqrt{n}}=\frac{p-1}{q-1}$. Choose $t_{1},t_{2}\geq 0$ so that
$t=t_{1}+t_{2}$ and $e^{-2t_{1}\sqrt{n}}=p-1$ and
$e^{-2t_{2}\sqrt{n}}=\frac{1}{q-1}$. Then we have

$\displaystyle\|e^{-t\sqrt{-\Delta}}f\|_{q}=\|e^{-t_{2}\sqrt{-\Delta}}(e^{-t_{1}\sqrt{-\Delta}}f)\|_{q}\leq\|e^{-t_{1}\sqrt{-\Delta}}f\|_{2}\leq\|f\|_{p}.$

In what follows we assume $2\leq p\leq q$. We will use a standard argument to
deduce the validity of the hypercontractivity estimate from a log Sobolev
inequality. Nonnegativity of the kernel for the Poisson semigroup combined
with the triangle inequality implies $|e^{-t\sqrt{-\Delta}}f|\leq
e^{-t\sqrt{-\Delta}}|f|$ for any $f$. Thus by continuity and standard density
arguments we can assume that $f\geq 0$, $f$ is not identically zero, and $f$
is smooth in $(\ref{hyp22})$.

The equality $e^{-2t\sqrt{n}}=\frac{p-1}{q-1}$ implies
$q=1+e^{2t\sqrt{n}}(p-1)$. Fix $p\geq 2$ and consider the map

$\displaystyle\varphi(t)=\|e^{-t\sqrt{-\Delta}}f\|_{q(t)}>0,\quad t\geq 0,$

where $q(t)=1+e^{2t\sqrt{n}}(p-1)$. If we show $\varphi^{\prime}(t)\leq 0$,
then we obtain $\varphi(t)\leq\varphi(0)=\|f\|_{p}$, and this proves the
sufficiency part. Let $\psi(t)=\ln\varphi(t)$. We have

$\displaystyle\frac{q^{2}}{q^{\prime}}\psi^{\prime}(t)=-\ln\left(\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q}d\sigma_{n}\right)+\frac{\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q}\left(\ln(e^{-t\sqrt{-\Delta}}f)^{q}+\frac{q^{2}}{q^{\prime}}\frac{\partial_{t}e^{-t\sqrt{-\Delta}}f}{e^{-t\sqrt{-\Delta}}f}\right)d\sigma_{n}}{\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q}d\sigma_{n}}.$

Clearly $\psi^{\prime}\leq 0$ if and only if

$\displaystyle\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q}\ln(e^{-t\sqrt{-\Delta}}f)^{q}d\sigma_{n}-\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q}d\sigma_{n}\ln\left(\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q}d\sigma_{n}\right)$
$\displaystyle\leq\frac{q^{2}}{q^{\prime}}\int_{\mathbb{S}^{n}}(e^{-t\sqrt{-\Delta}}f)^{q-1}\sqrt{-\Delta}(e^{-t\sqrt{-\Delta}}f)d\sigma_{n}.$

Let $g=e^{-t\sqrt{-\Delta}}f\geq 0$. Then we can rewrite the previous
inequality as

(7) $\displaystyle\int_{\mathbb{S}^{n}}g^{q}\ln
g^{q}d\sigma_{n}-\int_{\mathbb{S}^{n}}g^{q}d\sigma_{n}\ln\left(\int_{\mathbb{S}^{n}}g^{q}d\sigma_{n}\right)\leq\frac{q^{2}}{2(q-1)\sqrt{n}}\int_{\mathbb{S}^{n}}g^{q-1}\sqrt{-\Delta}gd\sigma_{n},$

where we used the fact that $q^{\prime}=2(q-1)\sqrt{n}$. Since
$e^{-t\sqrt{-\Delta}}$ is contractive in $L^{\infty}(\mathbb{S}^{n})$ with a
nonnegative, symmetric kernel, it follows that the validity of the estimate
(7) for $q=2$ implies (7) for all $q\in[2,\infty)$; see, e.g., Theorem 4.1 in
[3].

Let $g=\sum_{k\geq 0}H_{d}$ be the decomposition of $g$ into its spherical
harmonics. Then the estimate (7) for $q=2$ takes the form

$\displaystyle\int_{\mathbb{S}^{n}}g^{2}\ln
g^{2}d\sigma_{n}-\int_{\mathbb{S}^{n}}g^{2}d\sigma_{n}\ln\left(\int_{\mathbb{S}^{n}}g^{2}d\sigma_{n}\right)\leq\sum_{k\geq
0}2\sqrt{\frac{k(k+n-1)}{n}}\,\|H_{k}\|_{2}^{2}.$

It follows from Beckner’s conformal log Sobolev inequality [1] (which is a
consequence of Lieb’s sharp Hardy–Littlewood–Sobolev inequality [6]) that for
any smooth nonnegative $g=\sum_{k\geq 0}H_{k}$ we have

$\displaystyle\int_{\mathbb{S}^{n}}g^{2}\ln
g^{2}d\sigma_{n}-\int_{\mathbb{S}^{n}}g^{2}d\sigma_{n}\ln\left(\int_{\mathbb{S}^{n}}g^{2}d\sigma_{n}\right)\leq\sum_{k\geq
0}\Delta_{n}(k)\,\|H_{k}\|_{2}^{2}$

with $\Delta_{n}(k)=2n\sum_{m=0}^{k-1}\frac{1}{2m+n}$. Thus, the estimate (5)
is a consequence of the following lemma.

######  Lemma 2.1.

Let $n\in\\{2,3\\}$. Then for all integers $k\geq 1$ one has

$\displaystyle n\sum_{m=0}^{k-1}\frac{1}{2m+n}\leq\sqrt{\frac{k(k+n-1)}{n}}.$

###### Proof.

We first check the inequality for $k\leq 3$ by direct computation. Indeed, the
case $k=1$ is an equality. The case $k=2$ can be checked as follows,

$\displaystyle 1+\frac{n}{2+n}=\frac{2+2n}{2+n}\leq\sqrt{\frac{2+2n}{n}},$

which is true because $n(2+2n)\leq(2+n)^{2}$ holds for $n=2,3$. The case $k=3$
can be checked similarly:

$\displaystyle\frac{2+2n}{2+n}+\frac{n}{4+n}\leq\sqrt{\frac{6+3n}{n}}$

holds for $n=2,3$ (notice that this inequality fails for $n=4$).

Next, we assume $k\geq 4$. We have

$\displaystyle\sum_{m=0}^{k-1}\frac{1}{m+\frac{n}{2}}=\frac{2}{n}+\sum_{m=1}^{k-1}\frac{1}{m+\frac{n}{2}}\leq\frac{2}{n}+\int_{0}^{k-1}\frac{1}{x+\frac{n}{2}}dx=\frac{2}{n}+\ln\left(\frac{k+\frac{n}{2}-1}{\frac{n}{2}}\right).$

Thus it suffices to show

$\displaystyle\frac{2}{n}+\ln\left(\frac{k+\frac{n}{2}-1}{\frac{n}{2}}\right)-\frac{2}{n}\sqrt{\frac{k(k+n-1)}{n}}\leq
0.$

Notice that the left hand side, call it $h(k)$, is decreasing in $k$. Indeed,
we have

$\displaystyle
h^{\prime}(k)=\frac{1}{\frac{n}{2}+k-1}-\frac{2k+n-1}{n\sqrt{kn(k+n-1)}}\leq\frac{1}{\frac{n}{2}+k-1}-\frac{1}{\sqrt{kn}}\leq\frac{1}{2\sqrt{\frac{n}{2}(k-1)}}-\frac{1}{\sqrt{kn}}\leq
0.$

On the other hand, we have for $n=2,3$,

$\displaystyle
h(4)=\frac{2}{n}+\ln\left(\frac{6+n}{n}\right)-\frac{2}{n}\sqrt{\frac{12+4n}{n}}\leq
0.$

Indeed, if $n=2$, $h(4)=1+2\ln 2-\sqrt{10}<0$, and if $n=3$,
$h(4)=\frac{2+3\ln 3-4\sqrt{2}}{3}<0$. ∎

###  2.3. Counterexample to $\textup{(ii)}\Rightarrow\textup{(i)}$ in high
dimensions

Let $\lambda:=\frac{n-1}{2}$, and let $C_{d}^{(\lambda)}(x)$ be the Gegenbauer
polynomial

(8) $\displaystyle
C_{d}^{(\lambda)}(x)=\sum_{j=0}^{\lfloor\frac{d}{2}\rfloor}(-1)^{j}\frac{\Gamma(d-j+\lambda)}{\Gamma(\lambda)j!(d-2j)!}(2x)^{d-2j},$

where $\lfloor\frac{d}{2}\rfloor$ denotes the largest integer $m$ such that
$m\leq\frac{d}{2}$, and $\Gamma(x)$ is the Gamma function. Notice that if we
let $Y_{d}(\xi)=C_{d}^{(\lambda)}(\xi\cdot e_{1})$, where
$e_{1}=(1,0,\ldots,0)\in\mathbb{R}^{n+1}$, then $Y_{d}(\xi)$ is a spherical
harmonic of degree $d$ on $\mathbb{S}^{n}$. In particular, for $t\geq 0$ such
that $e^{-2t\sqrt{n}}=\frac{p-1}{q-1}$, the estimate
$\|e^{-t\sqrt{-\Delta}}f\|_{L^{q}(\mathbb{S}^{n})}\leq\|f\|_{L^{p}(\mathbb{S}^{n})}$
applied to $f=Y_{d}(\xi)$ is equivalent to the estimate

(9) $\displaystyle\frac{\|Y_{d}\|_{q}}{\|Y_{d}\|_{p}}\leq
e^{t\sqrt{d(d+n-1)}}=\left(\frac{q-1}{p-1}\right)^{\frac{1}{2}\sqrt{\frac{d(d+n-1)}{n}}}.$

Next, we need

######  Lemma 2.2.

For any $d\geq 0$ we have

(10)
$\displaystyle\lim_{n\to\infty}\,\frac{\|Y_{d}\|_{L^{q}(\mathbb{S}^{n},d\sigma_{n})}}{\|Y_{d}\|_{L^{p}(\mathbb{S}^{n},d\sigma_{n})}}=\frac{\|h_{d}\|_{L^{q}(\mathbb{R},d\gamma)}}{\|h_{d}\|_{L^{p}(\mathbb{R},d\gamma)}},$

where $d\gamma(y)=\frac{e^{-y^{2}/2}}{\sqrt{2\pi}}dy$ is the standard Gaussian
measure on the real line, and $h_{d}(x)$ is the probabilistic Hermite
polynomial

(11) $\displaystyle
h_{d}(x)=\sum_{j=0}^{\lfloor\frac{d}{2}\rfloor}\frac{(-1)^{j}d!}{j!(d-2j)!}\frac{x^{d-2j}}{2^{j}}.$

###### Proof.

Indeed, notice that

(12)
$\displaystyle\|Y_{d}\|_{p}^{p}=\int_{\mathbb{S}^{n}}|C_{d}^{(\lambda)}(\xi\cdot
e_{1})|^{p}d\sigma_{n}(\xi)=\int_{-1}^{1}|C_{d}^{(\lambda)}(t)|^{p}c_{\lambda}(1-t^{2})^{\lambda-\frac{1}{2}}dt,$

where
$c_{\lambda}=\frac{\Gamma(\lambda+1)}{\Gamma(\frac{1}{2})\Gamma(\lambda+\frac{1}{2})}$.
In particular, after the change of variables $t=\frac{s}{\sqrt{2\lambda}}$ in
(12), and multiplying both sides in (12) by $(d!/(2\lambda)^{d/2})^{p}$ we
obtain

$\displaystyle\left(\frac{d!}{(2\lambda)^{d/2}}\right)^{p}\|Y_{d}\|_{p}^{p}=\int_{\mathbb{R}}\left|\frac{d!}{(2\lambda)^{d/2}}C_{d}^{(\lambda)}\left(\frac{s}{\sqrt{2\lambda}}\right)\right|^{p}\frac{c_{\lambda}}{\sqrt{2\lambda}}\left(1-\frac{s^{2}}{2\lambda}\right)^{\lambda-\frac{1}{2}}\mathbbm{1}_{[-\sqrt{2\lambda},\sqrt{2\lambda}]}(s)ds,$

where $\mathbbm{1}_{[-\sqrt{2\lambda},\sqrt{2\lambda}]}(s)$ denotes the
indicator function of the set $[-\sqrt{2\lambda},\sqrt{2\lambda}]$. Notice
that by Stirling’s formula for any $j\geq 0$, and any $d\geq 0$ we have

(13)
$\displaystyle\lim_{\lambda\to\infty}\,\frac{1}{\lambda^{d-j}}\frac{\Gamma(d-j+\lambda)}{\Gamma(\lambda)}=1.$

Therefore, (11) and (8) together with (13) imply that for all $s\in\mathbb{R}$
we have

$\displaystyle\lim_{\lambda\to\infty}\frac{d!}{(2\lambda)^{d/2}}C_{d}^{(\lambda)}\left(\frac{s}{\sqrt{2\lambda}}\right)=h_{d}(s).$

Invoking Stirling’s formula again we have

$\displaystyle\lim_{\lambda\to\infty}\frac{c_{\lambda}}{\sqrt{2\lambda}}\left(1-\frac{s^{2}}{2\lambda}\right)^{\lambda-\frac{1}{2}}\mathbbm{1}_{[-\sqrt{2\lambda},\sqrt{2\lambda}]}(s)=\frac{e^{-s^{2}/2}}{\sqrt{2\pi}}\quad\text{for
all}\quad s\in\mathbb{R}.$

Finally, to apply Lebesgue’s dominated convergence theorem it suffices to
verify that for all $s\in\mathbb{R}$ and all $\lambda\geq\lambda_{0}$ we have
the following pointwise estimates

$\displaystyle
a)\quad\frac{c_{\lambda}}{\sqrt{2\lambda}}\left(1-\frac{s^{2}}{2\lambda}\right)^{\lambda-\frac{1}{2}}\mathbbm{1}_{[-\sqrt{2\lambda},\sqrt{2\lambda}]}(s)\leq
Ce^{-s^{2}/2}$ $\displaystyle
b)\quad\frac{d!}{(2\lambda)^{d/2}}C_{d}^{(\lambda)}\left(\frac{s}{\sqrt{2\lambda}}\right)\leq
c_{1}(d)(1+|s|)^{c_{2}(d)},$

where $\lambda_{0},C,c_{1}(d),c_{2}(d)$ are some positive constants
independent of $\lambda$ and $s$.

To verify a) it suffices to consider the case
$s\in[-\sqrt{2\lambda},\sqrt{2\lambda}]$. Since
$\lim_{\lambda\to\infty}\frac{c_{\lambda}}{\sqrt{2\lambda}}=\frac{1}{\sqrt{2\pi}}$
it follows that $\frac{c_{\lambda}}{\sqrt{2\lambda}}\leq C$ for all
$\lambda\geq\lambda_{0}$, where $\lambda_{0}$ is a sufficiently large number.
Next, the estimate $(1-\frac{s^{2}}{2\lambda})^{\lambda-1/2}\leq
C^{\prime}e^{-s^{2}/2}$ for $s\in[-\sqrt{2\lambda},\sqrt{2\lambda}]$ follows
if we show that $(1-\frac{1}{2\lambda})\ln(1-t)\leq
C^{\prime\prime}/\lambda-t$ for all $t:=\frac{s^{2}}{2\lambda}\in[0,1]$ where
$C^{\prime\prime}$ is a universal positive constant. The latter inequality
follows from $\ln(1-t)\leq-t$ for $t\in[0,1]$.

To verify b) it suffices to show that for all $\lambda\geq\lambda_{0}>0$ and
all integers $j$ such that $d\geq j\geq 0$ one has

$\displaystyle\frac{1}{\lambda^{d-j}}\frac{\Gamma(d-j+\lambda)}{\Gamma(\lambda)}\leq
C(d-j),$

where $C(d-j)$ depends only on $d-j$. The latter inequality follows from (13)
provided that $\lambda\geq\lambda_{0}$ where $\lambda_{0}$ is a sufficiently
large number.

Thus, it follows from the Lebesgue’s dominated convergence theorem that

$\lim_{n\to\infty}\frac{d!}{(n-1)^{d/2}}\|Y_{d}\|_{L^{p}(\mathbb{S}^{n},d\sigma_{n})}=\|h_{d}\|_{L^{p}(\mathbb{R},d\gamma)}.$

The lemma is proved. ∎

Now we fix $q>\max\\{p,2\\}$ and, in order to prove the failure of
$\textup{(ii)}\Rightarrow\textup{(i)}$ for all sufficiently large $n$, we
argue by contradiction and assume that there is a sequence of dimensions
$\\{n_{j}\\}_{j\geq 1}$ going to infinity such that
$\textup{(ii)}\Rightarrow\textup{(i)}$ in Theorem 1.1 does hold. Then, by
combining (9) and (10) we have

(14)
$\displaystyle\frac{\|h_{d}\|_{L^{q}(\mathbb{R},d\gamma)}}{\|h_{d}\|_{L^{p}(\mathbb{R},d\gamma)}}\leq\left(\frac{q-1}{p-1}\right)^{\frac{\sqrt{d}}{2}}.$

On the other hand, a consequence of the main result in [5] and the assumption
$q>\max\\{p,2\\}$ is that

$\displaystyle\lim_{d\to\infty}\left(\frac{\|h_{d}\|_{L^{q}(\mathbb{R},d\gamma)}}{\|h_{d}\|_{L^{p}(\mathbb{R},d\gamma)}}\right)^{1/d}=\left(\frac{q-1}{\max\\{p,2\\}-1}\right)^{\frac{1}{2}},$

which is in contradiction with (14).

######  Remark 2.1.

Let $B(x,y)$ be the Beta function. The estimate (9) for $p=2$ and $q=4$ takes
the form

(15)
$\displaystyle\int_{-1}^{1}|C_{d}^{(\frac{n-1}{2})}(t)|^{4}(1-t^{2})^{\frac{n-2}{2}}dt\leq
9^{\sqrt{\frac{d(d+n-1)}{n}}}\frac{(n-1)^{2}B(1/2,n/2)}{d^{2}(2d+n-1)^{2}B^{2}(n-1,d)},$

where we used the fact that
$\|Y_{d}\|^{2}_{L^{2}(\mathbb{S}^{n})}=\frac{n-1}{d(2d+n-1)B(n-1,d)}$. The
numerical computations show that the inequality (15) already fails for $d=7$
and $n=13$.

### Acknowledgements

Partial support through US National Science Foundation grants DMS-1363432 and
DMS-1954995 (R.L.F.) as well as DMS-2052645, DMS-1856486, and CAREER-
DMS-2052865, CAREER-DMS-1945102 (P.I.) is acknowledged.

## References

  * [1] W. Beckner, Sobolev inequalities, the Poisson semigroups, and analysis on the sphere $\mathbb{S}^{n}$. Proc. Natl. Acad. Sci. USA 89 (1992), 4816–4819.
  * [2] A. Bonami, Étude des coefficients de Fourier des fonctions de $L^{p}(G)$. Ann. Inst. Fourier 20 (1970), 335–420.
  * [3] L. Gross, Logarithmic Sobolev inequalities and contractivity properties of semigroups.In: Dell’Antonio G., Mosco U. (eds) Dirichlet Forms. Lecture Notes in Mathematics, vol 1563. Springer, Berlin, Heidelberg.
  * [4] S. Janson, On hypercontractivity for multipliers on orthogonal polynomials. Ark. Mat. 21 (1983), 97–110.
  * [5] L. Larsson-Cohn, $L_{p}$-norms of Hermite polynomials and an extremal problem on Wiener chaos. Ark. Mat. 40 (2002), 133–144.
  * [6] E. H. Lieb, Sharp constants in the Hardy–Littlewood–Sobolev and related inequalities. Ann. of Math. 118 (1983), no. 2, 349–374.
  * [7] C. Mueller, F. Weissler, Hypercontractivity for the Heat Semigroup for Ultraspherical Polynomials and on the $n$-Sphere. Journal of Functional Analysis 48 (1982), 252–282.
  * [8] E. Nelson, The free Markoff field. Journal of Functional Analysis 12 (1973), 211–227.
  * [9] O. Rothaus, Logarithmic Sobolev inequalities and the spectrum of Sturm–Liouville operators. Journal of Functional Analysis 39 (1980), 42–56.
  * [10] F. Weissler, Logarithmic Sobolev inequalities and hypercontractivity estimates on the circle. Journal of Functional Analysis 37 (1980), 218–234.


# Origin of correlated isolated flat bands in copper-substituted lead
phosphate apatite

Sinéad M. Griffin1,2 1Materials Sciences Division, Lawrence Berkeley National
Laboratory, Berkeley, California, 94720, USA  2Molecular Foundry Division,
Lawrence Berkeley National Laboratory, Berkeley, California, 94720, USA

(August 28, 2024)

###### Abstract

A recent report of room temperature superconductivity at ambient pressure in
Cu-substituted apatite (‘LK99’) has invigorated interest in the understanding
of what materials and mechanisms can allow for high-temperature
superconductivity. Here I perform density functional theory calculations on
Cu-substituted lead phosphate apatite, identifying correlated isolated flat
bands at the Fermi level, a common signature of high transition temperatures
in already-established families of superconductors. I elucidate the origins of
these isolated bands as arising from a structural distortion induced by the Cu
ions and a chiral charge density wave from the Pb lone pairs. These results
suggest that a minimal two-band model can encompass much of the low-energy
physics in this system. Finally, I discuss the implications of my results on
possible superconductivity in Cu-doped apatite.

††preprint: APS/123-QED

##  I Introduction

High-TC superconductors are arguably the holy grail of condensed matter
physics with huge potential applications for an energy-efficient future. The
first class of superconductors that were considered to be high-TC were the
cuprates which were discovered by Bednorz and Müller in 1987 [1]. The cuprates
have been subsequently followed by several new classes including the Fe-
pnictides in 2008 [2] and the nickelates [3].

While significant strides have been made in the discovery and understanding of
high-TC superconductors, and we continue to unearth novel examples within
established classes [4], a definitive roadmap to achieving room-temperature TC
under ambient pressures has remained elusive. Common to many of these high-TC
superconducting families are strongly correlated bands which can give rise to
unconventional mechanisms for Cooper pair formation [5, 6], and proximity to
multiple competing interactions such as antiferromagnetism, charge density
waves and spin-density waves. These phases can compete or coexist with
superconductivity where fluctuations between these states are believed to play
a significant role for achieving high-TC. Searching for these features in new
materials systems is therefore a promising route for finding new classes of
high-TC superconductors. For instance, the nickelate superconductors were
originally predicted in theory by Anisimov, Bukhvalov and Rice as an analogy
to the cuprate superconductors [7]. Similar approaches have also been proposed
to selectively design a material with the sought-after isolated d-manifold
that is associated with strong correlations [8], and have inspired high-
throughput searches for good candidate materials, further expanding the
horizons of high-Tc superconductivity [9].

The recent report of possible room temperature superconductivity at ambient
pressures in Cu-substituted apatite (also known as ‘LK99’)[10, 11] motivates
the need for a thorough understanding of the structure-property relationships
in these compounds to begin to unravel their potential correlated physics. In
this Letter, I use ab initio calculations to elucidate the key competing
interactions in Cu-doped apatite at the mean-field density functional level.

##  II Methods

I used the Vienna Ab initio Simulation Package (VASP) [12, 13, 14, 15] for all
density functional theory (DFT) calculations with full calculation details
given in the SI. I applied a Hubbard-U correction to account for the
underlocalization of the Cu-$d$ states. I tested values of U between 2 eV and
6 eV, finding my results were similar for all values calculated. The results
in the main text are for U = 4 eV which gives lattice parameters within 1% of
experiment [10, 16].

##  III Results

###  III.1 Structural Properties

Apatites are materials with the general formula A10(TO4)6X2±x, where A =
alkaline or rare earth metal; T = Ge, Si, or P; and X = halide, O, or OH. The
name ‘apatite’ derives from the Greek apatē meaning ‘deceit’ as a result of
the diverse range of forms it can take [17]. Here I consider the lead-
phosphate apatite Pb10(PO4)6(OH)2. Taking its structure reported from X-ray
diffraction in Ref. [16] as the starting point, its structure following a full
optimization is depicted in Fig. 1. It adopts the typical crystal structure of
various apatite chemistries, namely it forms a network comprising PbO6 prisms
that are corner shared with PO4 tetrahedra. I refer to these Pb as Pb(1), in
keeping with the convention in literature [18]. This framework is filled with
Pb6(OH)2 where the (OH)2 forms a chain in the center of a hexagonal structure
defined by these Pb(2) atoms. While here I consider OH-1 filling the hexagonal
column, I obtained similar results for O only in the column. Typically
apatites adopt the hexagonal $P6_{3}/m$ space group – here our resulting
structure has $P6_{3}$ owing to the breaking of reflection from the OH
molecule ordering. However, small structural deviations from the $P6_{3}/m$
commonly occur depending on its constituents.

Figure 1: (a) Lead-phosphate apatite structure with two inequivalent Pb sites
as described in the main text. Columns of O or OH sit in the center column
defined by Pb(2) hexagonal structure. (b) Calculated electronic localization
function for Pb10(PO4)6OH2. Oxygens surrounding Pb(2) are repelled by the lone
pair.

While the Pb(1) forms the overall framework with the PO4 tetrahedra, the Pb(2)
play a crucial role in Pb-O connectivity and polyhedra tilts throughout the
structure. In fact, while both Pb(1) and Pb(2) possess 6s2 lone pairs, only
the latter is stereochemically active. I verify this by calculating the
electronic localization function (ELF) as shown in Fig. 1(b), finding that the
Pb(2) lone pairs form a chiral arrangement that make an angle of $\approx
105^{\circ}$ with the a-axis (see SI for more details). The arrangement of the
lone pairs sets the resulting oxygen coordination of the Pb(2), namely its six
oxygens are arranged asymmetrically as a result of their repulsion from the
lone pair forming a chiral charge density wave. Since these oxygen are corner
shared with the PO4, the structural distortions associated with the Pb(2) lone
pair activity propagate throughout the structure. Such an effect has been
discussed previously in various apatite materials, noting that the lone pair
ordering differs depending on the specific system [18, 19].

I next consider the substitution of Cu on a Pb(1) site resulting in
CuPb9(PO4)6(OH)2, with the fully optimized structure shown in Fig. 2(b). I
find several changes to the structure with the inclusion of Cu. Firstly, all
lattice parameters decrease with $a$ going from 9.875 Å to 9.738 Å and $c$
going from 7.386 Å to 7.307 Å. While my calculated lattice parameters agree
well with those reported in Ref.[10], I find a greater structure collapse than
their report of $a$ going from 9.865 Å without Cu to 9.843 Å with Cu, and $c$
going from 7.431 Å without Cu to 7.428 Å with Cu. Interestingly, I find Cu
substitution results in a global structural distortion that results in a
change in coordination from nine to six (with a cutoff of $3$ Å) not only for
the Cu on the Pb(1) site, but also for all other Pb(1) sites. This results in
a modification of the polyhedral tilts throughout the structure, most notably
in the PO4 polyhedra, as seen in Fig. 2(b). I classify the distortion through
analysis of the symmetry-adapted phonon modes before and after Cu substitution
finding the structural distortion to be caused by $\Gamma_{1}$ and
$\Gamma_{2}$ modes of amplitudes 1.19 Å and 1.78 Å, respectively.
Visualization of these symmetry-adapted modes is given in the SI (Fig. S1);
this analysis confirms that the structural distortion caused by the Cu
substitution on Pb(1) is primarily driven by polyhedral tilts of the PO4 and
their corner-shared oxygen neighbors.

Figure 2: (a) Lead-phosphate apatite structure showing nine-coordinated Pb(1)
sites. (b) Cu-substituted structure showing six-coordinated Cu and Pb(1) sites
with distorted trigonal prism coordination with two different bondlengths are
a rigid twist of $\approx 24^{\circ}$ between the upper and lower triangles. A
cartoon of the resulting crystal-field diagram for Cu-$d^{9}$ is given on the
right.

Looking closely at the change in coordination of Cu on the Pb(1) site, I see
the Cu2+ is now six-coordinated with oxygen forming a distorted Jahn-Teller
trigonal prism. In an ideal trigonal prism the anions are arranged at the
vertices of a regular triangular prism with the transition metal at the center
with $D_{3h}$ symmetry. However, breaking in-plane mirror symmetry results in
Jahn-Teller distorted trigonal crystal field with $C_{3v}$ symmetry, such as
is found in the Janus compounds [20]. Here I find the mirror is already broken
by the next-nearest neighbor P ions – the Cu-O distance that has nearby P is
$2.06$ Å, whereas the Cu-O without surround P ions is $2.35$ Å. Such out-of-
plane asymmetry in the Cu environment results in a local dipole in $z$ which
can also influence the resulting electronic structure [20]. In addition to the
mirror symmetry breaking, the O triangles are also rotated at a small angle of
$\sim 24^{\circ}$ with respect to each other, referred to as a Bailar twist in
molecules[21], giving the final distorted trigonal prism. Interestingly, the
same structural distortion has been observed upon substitution of U on the
Ca(1) site in fluorapatite Ca10(PO4)6F2 as measured from X-ray absorption
spectroscopy [22]. There the authors found that the U substitution on the
Ca(1) sites causes an usual 6-coordinated structure – a trigonal prism with
six equal bondlengths of 2.06Å with an undetermined Bailar twist angle,
however.

###  III.2 Electronic Structure

Figure 3: Calculated spin-polarized electronic band structure (left) and
corresponding density of states (right). The spin-up bands are depicted in
solid orange, and the spin-down bands are dashed blue. The total density of
states is shaded grey with projections shown of the Cu-$d$ orbitals (pink) and
its neighboring O-$p$ orbitals (green). In both plots the Fermi level is set
to 0 eV and is marked by the dashed line.

I present the calculated spin-polarized electronic structure in Fig. 3.
Remarkably, I find an isolated set of flat bands crossing the Fermi level,
with a maximum bandwidth of $\sim$130 meV (see Fig.4) that is separated from
the rest of the valence manifold by 160 meV. Such a narrow bandwidth is
particularly indicative of strongly correlated bands. These especially flat
bands are consistent with the Cu-O coordination where I find Cu-O bondlengths
of 2.35 Å and 2.06 Å in the distorted trigonal prism. For comparison, the Cu-O
bondlengths in the cuprate superconductors are typically $\leq$ 2Å (in-plane)
and $\leq$ 2.3Å (apical) [3], giving further evidence of the unusual
coordination and resulting band localization in this isolated Cu-$d$ manifold.

Figure 4: Calculated spin-polarized electronic band structure in smaller
energy range around the Fermi level showing the isolated two-band Cu-$d$
manifold. The Fermi level is set to 0 eV and is marked by the dashed line.

The crystal-field splitting for the trigonal prism comprises a single
$d_{z^{2}}$, doubly degenerate $d_{xy}$ and $d_{x^{2}-y^{2}}$, and doubly
degenerate $d_{yz}$ and $d_{xz}$. In fact, the relative positioning of the
$d_{z^{2}}$ and the $d_{yz}/d_{xz}$ is set by the anisotropy induced by the
mirror-symmetry breaking – large values of asymmetry can cause the $d_{z^{2}}$
to be destabilized. In our case for Cu2+ with a $d^{9}$ configuration I expect
half filling of the doubly degenerate $d_{yz}/d_{xz}$ bands – this is
corroborated with our calculations in Fig. 3 where I find two bands of
$d_{yz}/d_{xz}$ character at the Fermi level that are are half filled. My
results suggest that low-energy physics of this system can be described by a
two-band $d_{yz}/d_{xz}$ model, similar to that originally suggested for Fe-
pnictide superconductors [23, 24]. However, unlike other correlated-$d$ band
superconductors, in this system the Cu-$d$ bands are particularly flat – there
is minimal band broadening from neighboring oxygen ions. If previous
assumptions about band flatness driving superconductivity are correct, then
this result would suggest a much more robust (higher temperature)
superconducting phase exists in this system, even compared to well-established
high-TC systems.

While the calculations presented here are performed for GGA+U with a U = 4 eV
applied to the Cu-$d$, I find the same qualitative features of the band
structures for a wide range of U values, as presented in the SI. This suggests
that the presence of such an isolated manifold of flat Cu-$d$ bands is a
feature of the structural reconstruction and new crystal field environment
provided by the apatite network. I further confirm this by calculating the
band structure of the Cu-substituted compound without performing a structural
relaxation, that is, taking the structure of the ideal Pb10(PO4)6(OH)2 and
replacing a Pb(1) with Cu without allowing the structure to relax. In this
case, I find slight spin-polarization in the band structure, however all
Cu-$d$ states are now in the bulk valence manifold and do not form an isolated
manifold (see SI Fig. S3). Finally, the top of the valence band is made up of
the Pb(2)-s states, confirming their contribution to the stereochemical
activity of the lone pairs (see SI).

I next investigate the exchange interactions between Cu ions in different unit
cells by constructing doubled unit cells in the in-plane and out-of-plane
directions. I find that there is a preference of 2 meV/Cu for ferromagnetic
coupling in the out-of-plane direction (with a Cu-Cu separation of $c=7.307$
Å), whereas a preference of 7 $\mu$eV/Cu for antiferromagnetic coupling in the
in-plane direction (with a Cu-Cu separation of $a=9.738$Å). However, while
this result is indicative of the potential exchange interactions in the
system, it relies on the unrealistic assumption that the Cu ions will sit on
the same substitution position in each unit cell. A full study of potential
exchange interactions for various Cu locations is out of the scope of this
work, and will crucially be informed by the Cu locations determined in
experiment.

So far all calculations have been for Cu on the Pb(1) site, which is the site
occupancy reported in Ref. [10]. I also calculated the structural and
electronic properties of Cu in the Pb(2) site with the resulting structure and
bands given in the SI. In this location, the Cu interrupts the hexagonal
network that is characteristic of the apatite structure, which causes a
structural rearrangement to the lower $P1$ symmetry, which would have a
profound effect on structural measures such as x-ray diffraction. The Cu
substituted in this position is now tetrahedrally coordinated by oxygen, and
has a significantly different electronic structure with no correlated $d$
bands crossing the Fermi level, as detailed in the SI. In fact, I find that Cu
on this Pb(2) is 1.08 eV more energetically favorable than Cu on the Pb(1)
site, suggesting possible difficulties in robustly obtaining Cu substituted on
the Pb(1) site.

##  IV Discussion

These theoretical results suggest that the apatite structure provides a unique
framework for stabilizing highly localized Cu-$d^{9}$ states that form a
strongly correlated flat band at the Fermi level. The central role of
stereochemically active 6s2 lone pairs of Pb(2) manifests in the formation of
a chiral charge density wave and the propagation of structural distortions
with connected polyhedra. When Cu is substituted on a Pb(1) site, the result
is a cascade of structural alterations, including reduced lattice parameters,
changes in coordination, and modified polyhedral tilts, leading to a local
Jahn-Teller distorted trigonal prism around Cu. This results in an unusually
flat set of isolated $d_{yz}/d_{xz}$ bands with half-filling.

I briefly note that achieving such a crystal field environment should also be
possible in intercalated twisted heterogeneous bilayers where selection of
different heterobilayers can provide the mirror symmetry breaking, while moiré
twist can provide an arbitrary rotation of the upper and lower triangles. In
fact such a platform would be ideal for probing the physics found here given
its broad range of tunability and the state-of-the-art characteriza- tion
probes for their interrogation [25].

I now discuss the potential implications of these features of Cu-substituted
apatite for possible high-TC superconductivity, and in particular the role of
the flat bands and the presence of competing magnetic interactions and
phonons. Isolated, flat bands have long been a target for achieving high-TC as
predicted from BCS theory[26, 27, 28]. The critical temperature TC given by by
T${}_{C}\propto\exp(-1/|U|\rho_{0}(E_{F}))$, where $|U|$ is the magnitude of
the effective attractive interaction and $\rho_{0}(E_{F})$ is the density of
states at the Fermi surface. In a flat band where the density of states
diverges, TC is proportional to $|U|$. This suggests that, at low interaction
strengths, the critical temperature in flat bands could be significantly
enhanced, as is currently being explored in various moiré systems[29]. In our
case, the correlated bands have a bandwidth of $\sim$ 130 meV which is less
than their separation from the rest of the valence manifold (160 meV) at the
level of density functional theory. While strong correlations will need to be
taken into account adequately to accurately describe these scenarios, these
are promising hints for these proposals that predict that the TC is
proportional to the interaction strength.

In contrast to conventional superconductors, where phonon-mediated
interactions govern Cooper pair formation, there is wide consensus that other
bosonic excitations are responsible for the attractive pair formation in
high-$T_{C}$ superconductors ranging from paramagnons to spin- and charge-
density waves [30]. In this system, I have identified several potential
sources of fluctuations that could contribute to pairing. Firstly, I
identified a charge density wave (CDW) driven by chiral lone pair ordering on
the Pb(2) sites – the presence of this CDW is strongly connected to the
structural rearrangement that occurs when Cu is incorporated into the Pb(1)
lattice sites. In addition to this, I identified two zone-center phonon modes
that trigger the global structural deformation that occurs as a result of the
Cu substitution, suggesting potentially strong electron-phonon coupling for
these modes. Finally, I calculated the relative exchange interactions between
Cu in neighboring unit cells. Interestingly for the out-of-plane coupling,
that is, along the Cu-Pb-Cu one-dimensional chains, I find ferromagnetic
coupling is favored by 2 meV/Cu over antiferromagnet coupling, even though the
Cu are over 7 Å apart, suggesting that spin fluctuations could also play a key
role.

Finally, the calculations presented here suggest that Cu substitution on the
appropriate (Pb(1)) site displays many key characteristics for high-TC
superconductivity, namely a particularly flat isolated $d$-manifold, and the
potential presence of fluctuating magnetism, charge and phonons. However,
substitution on the other Pb(2) does not appear to have such sought-after
properties, despite being the lower-energy substitution site. This result
hints to the synthesis challenge in obtaining Cu substituted on the
appropriate site for obtaining a bulk superconducting sample. Nevertheless, I
expect the identification of this new material class to spur on further
investigations of doped apatite minerals given these tantalizing theoretical
signatures and experimental reports of possible high-TC superconductivity.

#### Acknowledgements-.

I am grateful to David Prendergast, Adam Schwartzberg, Shuhada’ Sadaqat and
John Vinson for insightful discussions and encouragement. I also thank D.
Kwabena Bediako, Donny Evans, Katherine Inzani, Adam Schwartzberg and John
Vinson for feedback on this draft, and Jeff Neaton and Richard Martin for
further comments. This work was funded by the U.S. Department of Energy,
Office of Science, Office of Basic Energy Sciences, Materials Sciences and
Engineering Division under Contract No. DE-AC02-05-CH11231 within the Theory
of Materials program. Computational resources were provided by the National
Energy Research Scientific Computing Center and the Molecular Foundry, DOE
Office of Science User Facilities supported by the Office of Science, U.S.
Department of Energy under Contract No. DEAC02-05CH11231. The work performed
at the Molecular Foundry was supported by the Office of Science, Office of
Basic Energy Sciences, of the U.S. Department of Energy under the same
contract.

##  V Supplemental Material

###  V.1 Calculation Details

I used the Vienna Ab initio Simulation Package (VASP) [12, 13, 14, 15] for all
density functional theory (DFT) calculations with projector augmented wave
(PAW) pseudopotentials [31, 32] including Pb 5d106s26p2, Cu 3d104s1, P 3s23p5,
O 2s22p4, and H 1s1 as valence electrons. I used a plane wave cut-off energy
of 600 eV and a $6\times 6\times 8$ Gamma-centered k-point grid for structural
optimizations and a $10\times 10\times 12$ Gamma-centered k-point grid for the
density of states and electron localization functional calculations. All
calculations are done using the generalized gradient approximation (GGA) based
exchange-correlation functional PBEsol [33], with a Hubbard-U of 4 eV applied
to the Cu-d states as implemented with the Dudarev approach unless otherwise
stated [34]. This gives lattice parameters within 1% of experiment [10] as
summarized in Table LABEL:table-structure. The electronic convergence
criterion is set to $10^{-6}$ eV and the force convergence criterion is set to
0.01 eV / Å. We included spin-orbit coupling self-consistently in the
electronic structure calculations where specified. Electronic structure plots
were generated using the sumo software package [35].

###  V.2 Structural Distortion

We used the AMPLIMODES software from the Bilbao Crystallographic Server[36] to
analyze the symmetry-adapted phonon modes in going from the stoichiometric
Pb10(PO4)6(OH)2 with space group $P6_{3}$ to the Cu-substituted compound with
Cu on the Pb(1) site, resulting in space group $P3$. A summary of the
resulting distortion is given in Table 1 with a visualization of the two
modes, $\Gamma_{1}$ and $\Gamma_{2}$ depicted in Fig. 10.

Irrep | K-vector | Direction | Subgroup | Amplitude (Å)
---|---|---|---|---
$\Gamma_{1}$ | (0,0,0) | (a) | $P6_{3}$ | 1.192
$\Gamma_{2}$ | (0,0,0) | (a) | $P3$ | 1.778
Table 1: Symmetry-adapted distortion information in going from stoichiometric
Pb10(PO4)6(OH)2 to Cub9(PO4)6(OH)2 where Cu is on the Pb(1) site. Figure 5:
Summary of displacements of symmetry-adapted modes upon Cu substitution on the
Pb(1) site. (a) and (b) show the $\Gamma_{1}$ mode which retains the $P6_{3}$
space group. (c) and (d) show $\Gamma_{2}$ mode which reduces the $P6_{3}$
space group to $P3$. The arrows show the ionic motion under each phonon mode
with their length being proportional to the amplitude of the displacement.
Figure 6: (Top panel) Calculated spin-polarized band structure for different
values of U on Cu-d states for CuPb9(PO4)6(OH)2 with (a) U = 2 eV, (b) U = 4
eV, and (c) U = 6 eV. In all plots the majority spin channel is shown as a
solid orange line, and the minority as a dashed blue line. (Middle panel)
Zoomed-in version of the top plot of spin-polarized bands with a narrower
energy range. (Bottom panel) Calculated band structures without spin-
polarization. The Fermi level is set to 0 eV in all plots and is marked by the
dashed line.  Figure 7: Calculated spin-polarized band structures for (a)
stoichiometric Pb10(PO4)6(OH)2 and (b) CuPb9(PO4)6(OH)2 without structural
optimization, that is, with the structure of (a) with one Cu replaced for
Pb(1). In all plots the majority spin channel is shown as a solid orange line,
and the minority as a dashed blue line. The Fermi level is set to 0 eV in all
plots and is marked by the dashed line.  Figure 8: Calculated electronic
localization function for a 2$\times$2$\times$1 supercell of Pb10(PO4)6OH2
showing chiral charge density wave induced by Pb(2) lone pairs. Figure 9:
Results for CuPb9(PO4)6(OH)2 with Cu on the Pb(2) site. (a) The resulting
fully relaxed structure for Cu on the Pb(2) site. The inset shows the new
local environment of the Cu ion and resulting Cu-O bond lengths. (b) The
corresponding calculated spin-polarized electronic structure showing the bands
(left) and spin-polarized density of states (right). The total and orbital-
projected density of states are shown for the Cu-d and O-p states. The Fermi
level is set to 0 eV in all plots and is marked by the dashed line. . Figure
10: alculated non-spin-polarized electronic structures for (a) CuPb9(PO4)6O
and (b) CuPb9(PO4)6(OH)2. For both calculations the Cu is placed on the Pb(1)
site. The total and orbital-projected density of states are shown for replaced
for Pb(1). In all plots the majority spin channel is shown as a solid orange
line, and the minority as a dashed blue line. The Fermi level is set to 0 eV
in all plots and is marked by the dashed line. Note that these were both
calculated using the fully relaxed structures that included spin-polarization.

##  VI Structure Files

The resulting relaxed structures in POSCAR format are given below. The
specific geometry is given in the headers.

    Cu on Pb(1) for OH columns
    1.00000
       9.73753   0.00000  -0.00000
      -4.86876   8.43294   0.00000
       0.00000   0.00000   7.30653
    Cu Pb P O H
    1 9 6 26 2
    Direct(44) [A1B9C6D26E2]
       0.33333   0.66666   0.99594  Cu
       0.66666   0.33333   0.97556  Pb
       0.66666   0.33333   0.47828  Pb
       0.33333   0.66666   0.47831  Pb
       0.24304   0.99176   0.23619  Pb
       0.74033   0.00357   0.74019  Pb
       0.00823   0.25127   0.23619  Pb
       0.99642   0.73675   0.74019  Pb
       0.74872   0.75695   0.23619  Pb
       0.26324   0.25966   0.74019  Pb
       0.39800   0.37662   0.22932  P
       0.59606   0.63623   0.76062  P
       0.62337   0.02138   0.22932  P
       0.36376   0.95982   0.76062  P
       0.97861   0.60199   0.22932  P
       0.04017   0.40393   0.76062  P
       0.33463   0.49455   0.22138  O
       0.63904   0.50466   0.73061  O
       0.50544   0.84007   0.22138  O
       0.49533   0.13438   0.73061  O
       0.15992   0.66536   0.22138  O
       0.86561   0.36095   0.73061  O
       0.58275   0.46627   0.20829  O
       0.42113   0.55711   0.83037  O
       0.53372   0.11647   0.20829  O
       0.44288   0.86401   0.83037  O
       0.88352   0.41724   0.20829  O
       0.13598   0.57886   0.83037  O
       0.32305   0.24965   0.07561  O
       0.70965   0.76259   0.89995  O
       0.75034   0.07339   0.07561  O
       0.23740   0.94705   0.89995  O
       0.92660   0.67694   0.07561  O
       0.05294   0.29034   0.89995  O
       0.61424   0.72712   0.57832  O
       0.36008   0.28424   0.41289  O
       0.27287   0.88712   0.57832  O
       0.71575   0.07584   0.41289  O
       0.11287   0.38575   0.57832  O
       0.92415   0.63991   0.41289  O
       0.00000   0.00000   0.70585  O
       0.00000   0.00000   0.30956  O
       0.00000   0.00000   0.44500  H
       0.00000   0.00000   0.84006  H

    Cu on Pb(2) for OH columns
    1.00000
    9.72342   0.04200   0.00197
    -4.82301   8.48804   0.00345
    0.00025   0.00382   7.33474
    Cu   Pb   P    O    H
    1     9     6    26     2
    Direct
    0.29671   0.05116   0.33235
    0.32381   0.67469   0.98111
    0.65139   0.32730   0.98009
    0.66758   0.35114   0.48809
    0.32169   0.66632   0.47791
    0.75782   0.00049   0.75985
    0.99392   0.24366   0.24551
    0.99739   0.77380   0.76869
    0.74898   0.76332   0.23058
    0.22801   0.23596   0.76545
    0.40018   0.38468   0.22144
    0.59137   0.62723   0.73080
    0.61312   0.02576   0.23816
    0.37685   0.97309   0.72186
    0.96697   0.59735   0.22558
    0.01790   0.40361   0.72964
    0.31214   0.48018   0.22872
    0.66266   0.51575   0.73906
    0.51223   0.84427   0.23528
    0.47793   0.15426   0.72594
    0.14692   0.66732   0.22652
    0.83604   0.33147   0.73496
    0.58332   0.49499   0.22798
    0.40826   0.53051   0.72984
    0.50719   0.10481   0.25528
    0.48254   0.89402   0.72791
    0.88433   0.41328   0.23340
    0.10201   0.58649   0.73148
    0.36287   0.28748   0.04392
    0.64492   0.73807   0.90133
    0.70876   0.08669   0.05709
    0.26165   0.90861   0.88557
    0.90639   0.64180   0.05051
    0.07638   0.34936   0.89853
    0.65158   0.73564   0.56080
    0.35772   0.27368   0.38830
    0.27627   0.91556   0.54365
    0.73438   0.08266   0.39631
    0.07010   0.34970   0.55607
    0.90705   0.65323   0.39062
    0.00689   0.00336   0.92206
    0.08251   0.98342   0.30799
    0.03614   0.87556   0.35825
    0.04711   0.99828   0.04551

    Cu on Pb(1) for O columns
        1.00000
        9.62586   -0.00000   -0.00000
       -4.81293    8.33624    0.00000
        0.00000    0.00000    7.21588
       Cu   Pb   P    O
         1     9     6    25
    Direct
      0.66667  0.33333  0.98511
      0.00692  0.77593  0.26511
      0.99756  0.23907  0.75708
      0.22407  0.23099  0.26511
      0.76093  0.75849  0.75708
      0.76901  0.99308  0.26511
      0.24151  0.00244  0.75708
      0.33333  0.66667  0.01001
      0.33333  0.66667  0.50100
      0.66667  0.33333  0.49255
      0.62789  0.59460  0.23001
      0.38566  0.39352  0.76637
      0.40540  0.03329  0.23001
      0.60648  0.99214  0.76637
      0.96671  0.37211  0.23001
      0.00786  0.61434  0.76637
      0.49081  0.63441  0.25509
      0.48303  0.30416  0.77686
      0.36559  0.85640  0.25509
      0.69584  0.17887  0.77686
      0.14360  0.50919  0.25509
      0.82113  0.51697  0.77686
      0.75643  0.71787  0.09385
      0.25828  0.33902  0.92184
      0.28213  0.03856  0.09385
      0.66098  0.91926  0.92184
      0.96144  0.24357  0.09385
      0.08074  0.74172  0.92184
      0.29192  0.36023  0.58153
      0.71402  0.61660  0.42013
      0.63977  0.93170  0.58153
      0.38340  0.09742  0.42013
      0.06830  0.70808  0.58153
      0.90258  0.28598  0.42013
      0.56029  0.42167  0.15861
      0.49845  0.57912  0.78156
      0.57833  0.13861  0.15861
      0.42088  0.91933  0.78156
      0.86139  0.43971  0.15861
      0.08067  0.50155  0.78156
      0.00000  0.00000  0.35299

    Cu on Pb(2) for O columns
    1.00000
    9.78803 0.12357 0.00526
    -4.78808 8.26801 -0.00858
    0.00359 -0.00485 7.37499
    Cu Pb P O
    1 9 6 25
    Direct
    0.95321 0.73690 0.35934
    0.00223 0.25142 0.75098
    0.21445 0.19820 0.27381
    0.74490 0.74531 0.76387
    0.77953 0.98644 0.28067
    0.25110 0.00434 0.75792
    0.68083 0.34384 0.00125
    0.32306 0.66524 0.00484
    0.33825 0.67242 0.48957
    0.65744 0.32831 0.49365
    0.61175 0.58921 0.23731
    0.37015 0.39196 0.74887
    0.41117 0.02923 0.24495
    0.60851 0.98175 0.74867
    0.97935 0.39884 0.25463
    0.02783 0.63040 0.73552
    0.50362 0.66504 0.24183
    0.47249 0.30930 0.75483
    0.34797 0.84661 0.24748
    0.69473 0.16747 0.74982
    0.16028 0.50174 0.25577
    0.84707 0.53213 0.73907
    0.70683 0.63213 0.05889
    0.26126 0.35126 0.91635
    0.34673 0.07836 0.07899
    0.65481 0.91395 0.91340
    0.92330 0.31167 0.06924
    0.09741 0.73767 0.90429
    0.25995 0.33621 0.58013
    0.72714 0.65832 0.40169
    0.65107 0.91793 0.57478
    0.36022 0.08525 0.42001
    0.08931 0.74346 0.56613
    0.92479 0.26882 0.40749
    0.51861 0.40431 0.25717
    0.47997 0.57775 0.73775
    0.59365 0.12144 0.24141
    0.42729 0.91309 0.75013
    0.90084 0.50623 0.28580
    0.09909 0.51545 0.72278
    0.99782 0.95063 0.34472

## References

  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Kamihara et al. [2008] Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, Journal of the American Chemical Society 130, 3296 (2008).
  * Li et al. [2019] D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Nature 572, 624 (2019).
  * Pan et al. [2022] G. A. Pan, D. Ferenc Segedin, H. LaBollita, Q. Song, E. M. Nica, B. H. Goodge, A. T. Pierce, S. Doyle, S. Novakov, D. Córdova Carrizales, et al., Nature materials 21, 160 (2022).
  * Khodel and Shaginyan [1990] V. Khodel and V. Shaginyan, Jetp Lett 51, 553 (1990).
  * Heikkilä and Volovik [2016] T. T. Heikkilä and G. E. Volovik, Basic Physics of Functionalized Graphite pp. 123–143 (2016).
  * Anisimov et al. [1999] V. Anisimov, D. Bukhvalov, and T. Rice, Physical Review B 59, 7901 (1999).
  * Griffin et al. [2016] S. M. Griffin, P. Staar, T. C. Schulthess, M. Troyer, and N. Spaldin, Physical Review B 93, 075115 (2016).
  * Isaacs and Wolverton [2019] E. B. Isaacs and C. Wolverton, Physical Review X 9, 021042 (2019).
  * Lee et al. [2023a] S. Lee, J.-H. Kim, and Y.-W. Kwon, arXiv preprint arXiv:2307.12008 (2023a).
  * Lee et al. [2023b] S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, arXiv preprint arXiv:2307.12037 (2023b).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993), ISSN 01631829.
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994), ISSN 0163-1829, URL http://link.aps.org/doi/10.1103/PhysRevB.49.14251.
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Comp. Mater. Sci. 6, 15 (1996a), ISSN 09270256\.
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b), ISSN 0163-1829, URL http://link.aps.org/doi/10.1103/PhysRevB.54.11169.
  * Brückner et al. [1995] S. Brückner, G. Lusvardi, L. Menabue, and M. Saladini, Inorganica Chimica Acta 236, 209 (1995).
  * Roycroft and Cuypers [2015] P. D. Roycroft and M. Cuypers, Irish Journal of Earth Sciences 33, 71 (2015).
  * Peet et al. [2017] J. Peet, A. Piovano, M. Johnson, and I. R. Evans, Dalton Transactions 46, 15996 (2017).
  * Antao and Dhaliwal [2018] S. M. Antao and I. Dhaliwal, Journal of Synchrotron Radiation 25, 214 (2018).
  * Er et al. [2018] D. Er, H. Ye, N. C. Frey, H. Kumar, J. Lou, and V. B. Shenoy, Nano letters 18, 3943 (2018).
  * Bailar Jr [1964] J. C. Bailar Jr, Prep. Inorg. React 1, 1 (1964).
  * Rakovan et al. [2002] J. Rakovan, R. J. Reeder, E. J. Elzinga, D. J. Cherniak, C. D. Tait, and D. E. Morris, Environmental Science & Technology 36, 3114 (2002).
  * Han et al. [2008] Q. Han, Y. Chen, and Z. Wang, Europhysics Letters 82, 37007 (2008).
  * Raghu et al. [2008] S. Raghu, X.-L. Qi, C.-X. Liu, D. J. Scalapino, and S.-C. Zhang, Phys. Rev. B 77, 220503 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.220503.
  * Van Winkle et al. [2023] M. Van Winkle, I. M. Craig, S. Carr, M. Dandu, K. C. Bustillo, J. Ciston, C. Ophus, T. Taniguchi, K. Watanabe, A. Raja, et al., Nature Communications 14, 2989 (2023).
  * Heikkilä et al. [2011] T. T. Heikkilä, N. B. Kopnin, and G. E. Volovik, JETP letters 94, 233 (2011).
  * Kopnin et al. [2011] N. Kopnin, T. Heikkilä, and G. Volovik, Physical Review B 83, 220503 (2011).
  * Hofmann et al. [2020] J. S. Hofmann, E. Berg, and D. Chowdhury, Physical Review B 102, 201112 (2020).
  * Törmä et al. [2022] P. Törmä, S. Peotta, and B. A. Bernevig, Nature Reviews Physics 4, 528 (2022).
  * Bertel and Menzel [2016] E. Bertel and A. Menzel, Symmetry 8, 45 (2016).
  * Blöchl [1994] P. E. Blöchl, Physical Review B 50, 17953 (1994), ISSN 0163-1829, URL https://link.aps.org/doi/10.1103/PhysRevB.50.17953.
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999), ISSN 0163-1829, URL https://link.aps.org/doi/10.1103/PhysRevB.59.1758.
  * Perdew et al. [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical Review Letters 100, 136406 (2008), ISSN 0031-9007, URL https://link.aps.org/doi/10.1103/PhysRevLett.100.136406.
  * Dudarev et al. [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998), ISSN 1098-0121.
  * Ganose et al. [2018] A. M. Ganose, A. J. Jackson, and D. O. Scanlon, Journal of Open Source Software 3, 717 (2018).
  * Perez-Mato et al. [2010] J. Perez-Mato, D. Orobengoa, and M. Aroyo, Acta Crystallographica Section A: Foundations of Crystallography 66, 558 (2010).