2022

[1]Guo-yun Shao

1]School of Physics, Xi’an Jiaotong University, Xi’an, 710049, China 2]School
of Physics, Peking University, Beijing, 100871, China

# Speed of sound in QCD matter at finite temperature and density

<EMAIL_ADDRESS>Xin-ran Yang  Chong-long Xie  Wei-bo He  [  [

###### Abstract

The speed of sound in QCD matter at finite temperature and density is
investigated within the Polyakov loop improved Nambu–Jona-Lasinio (PNJL)
model. The spinodal structure associated with the chiral first-order chiral
phase transition is considered to describe the continuous variation of the
speed of sound. The behaviors of the squared sound speed in different phases,
including the stable, metastable and unstable phases, are derived. The
relation between speed of sound and QCD phase transitions is systematically
explored. In particular, the boundary of vanishing sound velocity is derived
in the temperature-density phase diagram, and the region where the sound wave
equation being broken is pointed out. Some interesting features of speed of
sound under different definitions are also discussed.

###### keywords:

Speed of sound, Quark matter, Chiral phase transition

##  1 Introduction

Quark-gluon plasma (QGP) can be created in heavy-ion collision (HIC)
experiments at relativistic energies. A crucial topic relevant is to explore
the equation of state (EOS) and phase transition from QGP to hadronic matter.
The hydrodynamic simulation provides a method to study the EOS of QGP Song11 ;
Song112 ; Deb16 . During the space-time evolution of QCD matter, the speed of
sound is one of the crucial physical quantities. Its dependence on environment
(temperature, density, chemical potential, etc.) carries important information
in describing the evolution of the fireball and final observables. Recently,
the studies in Refs. Gardim20 ; Sahu21 ; Biswas20  show that the speed of
sound as a function of charged particle multiplicity $\langle
dN_{ch}/d\eta\rangle$ can be extracted from heavy-ion collision data. In Ref.
Sorensen21  the authors try to build a connection between the sound speed and
baryon number cumulants to study the QCD phase structure.

The speed of sound in neutron star has also received a lot of attention (e.g.,
Refs. Reed20 ; Kanakis20 ; Han20 ). The density dependent behavior of sound
velocity influences the mass-radius relation, the tidal deformability, and
provides a sensitive probe of the EOS of neutron star matter. To obtain a two
solar mass neutron star, some studies find that it is essential for neutron
star matter to have a density range where the EOS is very stiff and the
corresponding squared speed of sound is significantly larger than $1/3$ Tews18
; Greif19 ; Forbes19 ; Drischler20 ; Essick20 ; Han19 ; Kojo . The study in
Ref. Jaikumar21  indicates that the speed of sound is crucial for the
gravitational wave frequencies induced by the $g$-mode oscillation of a
neutron star. It is also interesting to study the gravitational wave induced
by the cosmic QCD phase transition in which the speed of sound plays a
significant role.

As an important quantum in describing the evolution of strongly interacting
matter, the relation between the speed of sound and QCD phase transition is
worth exploring. The speed of sound has been calculated, e.g., in lattice QCD
Aoki06 ; Borsanyi14 ; Bazavov14 ; Philipsen13 ; Borsanyi20 , (P)NJL model
Motta18 ; Ghosh06 ; Marty13 ; Deb16 ; Saha18 ; zhao20 , quark-meson coupling
model Schaefer10 ; Abhishek18 , hadron resonance gas (HRG) model Venugopalan92
; Bluhm14 , field correlator method (FCM) Khaidukov18 ; Khaidukov19  and
quasiparticle model Mykhaylova21 . In previous studies, the main focus is put
on the region of high temperature and vanishing or small chemical potential.
In Ref. he2022 , we give an intensive study on the speed of sound in QCD
matter in the full temperature-chemical potential phase diagram. The numerical
results indicate that the dependence of sound speed on temperature and
chemical potential is indicative of QCD phase transition.

However, only the sound speed in the stable phase is considered in Ref. he2022
. There are still some crucial issues that need to be clarified. First of all,
the spinodal structure may be involved in heavy-ion collision experiments with
the decrease of collision energy Mishustin1999 ; Randrup2004 ; Koch2005 ;
Sasaki2007 ; Sasaki2008 ; Randrup2009 ; Steinheimer2012 ; Li2016 ;
Steinheimer2016 ; Steinheimer2017 ; Shao2020 . A complete evolution of sound
speed in the metastable and unstable phases needs to be explored to give a
distinct description of the fireball expansion. Secondly, it is found that the
sound speed takes small values at the CEP and on the boundaries of the first-
order phase transition near the CEP he2022 . A question aroused is that where
is the boundary of vanishing speed of sound. Furthermore, the behavior of
sound speed in the temperature and density phase diagram is still not
explored.

On the other hand, the values of speed of sound under different conditions are
involved in dealing with different problems in nuclear physics, such as the
gravitational signal from cosmic QCD phase transition Ahmadvand2018 ;
Tenkanen2022 , the bulk viscosity of strongly interacting matter Deb16 , the
equation of state of neutron star matter Reed20 ; Kanakis20 ; Han20  and the
evolution of QGP in HIC experiments. In this work, we will give a systematic
study on the relation between speed of sound and QCD phase transitions at
finite temperature and density under several different constraint conditions.
This work is helpful in dealing with the physics problems mentioned above.

The paper is organized as follows. In Sec. II, we derive the formulae of speed
of sound under different definitions in the temperature and density space, and
then briefly introduce the 2+1 flavor PNJL quark model. In Sec. III, we
present the numerical results of squared sound speed and discuss the relations
with the QCD phase structure. A summary is finally given in Sec. IV.

##  2 Speed of sound and the PNJL quark model

The general definition of speed of sound is

$c^{2}_{X}=\left(\frac{\partial p}{\partial\epsilon}\right)_{X}.$ (1)

A specifying constant quantity $X$ is required to describe the propagation of
the compression wave through a medium. To indicate the different profiles of
QCD matter, $X$ can be chosen as $s/\rho_{B},s,\rho_{B},T,\mu_{B}$. Different
definitions of speed of sound are taken in practice in dealing with different
physics issues.

For a fireball created in relativistic heavy-ion collisions, it evolves with a
constant entropy density per baryon $s/\rho_{B}$ if it is taken as an ideal
fluid. Therefore, it is meaningful to calculate the speed of sound along the
isentropic curve

$c_{s/{\rho_{B}}}^{2}=\bigg{(}\frac{\partial
p}{\partial\epsilon}\bigg{)}_{s/{\rho_{B}}}.$ (2)

The dependence of $c_{s/{\rho_{B}}}^{2}$ on parameters, e.g., temperature and
density, can indicate the variation of sound speed during the evolution and
provide important knowledge of interaction, phase transition and the EOS of
QGP.

The speed of sound with constant baryon number density or entropy density are
taken in describing the intermediate process of a hydrodynamic evolution Deb16
,

$c_{\rho_{B}}^{2}=\bigg{(}\frac{\partial
p}{\partial\epsilon}\bigg{)}_{\rho_{B}}\quad\textrm{and}\quad
c_{s}^{2}=\bigg{(}\frac{\partial p}{\partial\epsilon}\bigg{)}_{s}.$ (3)

For example, the temporal derivatives of temperature and chemical potential
are functions of $c_{\rho_{B}}^{2}$ and $c_{s}^{2}$, as

$\partial_{0}\mu_{B}=-c_{s}^{2}\mu_{B}\,\mathbb{\nabla}\cdot\mathbb{u},$ (4)

and

$\partial_{0}T=-c_{\rho_{B}}^{2}T\,\mathbb{\nabla}\cdot\mathbb{u},$ (5)

where $\mathbb{u}$ denotes the space component of four-velocity. The values of
$c_{\rho_{B}}^{2}$ and $c_{s}^{2}$ are directly connected to the bulk
viscosity coefficient.

It is also interesting to calculate the sound speed with a fixed temperature
or chemical potential

$c_{T}^{2}=\bigg{(}\frac{\partial
p}{\partial\epsilon}\bigg{)}_{T},\,\,\,\,\,\,\,\,c_{\mu_{B}}^{2}=\bigg{(}\frac{\partial
p}{\partial\epsilon}\bigg{)}_{\mu_{B}},$ (6)

In Ref. Sorensen21  the authors estimate $c^{2}_{T}$ as a function of the
logarithmic derivative with respect to the baryon density of QCD matter, and
try to build a connection with the baryon number cumulants to aid in detecting
the QCD critical endpoint. Besides, $c^{2}_{T}$ is also usually taken to study
the speed of sound in neutron star matter.

In this study, we will explore the speed of sound under different definitions
in the full temperature-density space. Since the general definitions can only
be used to calculate the sound speed on special trajectories, it is necessary
to derive the corresponding formulae in terms of $T$ and $\rho_{B}$. With the
fundamental thermodynamic relations, the sound speed formulae under different
constraint conditions can be derived as

$\\!c^{2}_{s/\rho_{B}}\\!=\\!\frac{\\!s^{2}\\!+\\!\rho_{B}^{2}\\!\left[\\!\left(\frac{\\!\partial\mu_{B}}{\\!\partial\rho_{B}}\\!\right)_{\\!T}\\!\left(\frac{\partial
s}{\\!\partial
T}\\!\right)_{\\!\rho_{\\!B}}\\!-\\!\left(\frac{\partial\mu_{\\!B}}{\partial
T}\\!\right)_{\\!\rho_{\\!B}}\\!\left(\frac{\partial
s}{\partial\rho_{\\!B}}\\!\right)_{\\!T}\\!\right]\\!+\\!s\rho_{\\!B}\\!\left[\\!\left(\frac{\partial\mu_{\\!B}}{\partial
T}\\!\right)_{\\!\rho_{\\!B}}\\!-\\!\left(\frac{\partial
s}{\partial\rho_{\\!B}}\\!\right)_{\\!T}\right]}{\left(Ts+\mu_{B}\rho_{B}\right)\left(\frac{\partial
s}{\partial T}\right)_{\rho_{B}}},$ (7)
$c_{s}^{2}=\frac{\rho_{B}\left[\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}-\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}\right]-s\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}}{\mu_{B}\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}},$ (8)
$c_{\rho_{B}}^{2}=\frac{s+\rho_{B}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}}{T\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}},\quad\quad\quad
c_{T}^{2}=\frac{\rho_{B}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}}{T\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}+\mu_{B}},$ (9)

and

$c_{\mu_{B}}^{2}=\frac{s\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}}{T\left[\left(\frac{\partial
s}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}-\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}\right]-\mu_{B}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}}.$ (10)

The details for deriving these formulae are affiliated in the appendix A. The
above formulae are only correct for isospin symmetric matter. The
corresponding formulae will be much more complicated for isospin asymmetric
matter.

To demonstrate the relation between the speed of sound under different
definition and QCD phase structure, we take the 2+1 flavor PNJL quark model in
the calculation. The Lagrangian density is given by

$\displaystyle\mathcal{L}$ $\displaystyle\\!=$
$\displaystyle\\!\bar{q}(i\gamma^{\mu}D_{\mu}\\!+\\!\gamma_{0}\hat{\mu}\\!-\\!\hat{m}_{0})q\\!+\\!G\sum_{k=0}^{8}\big{[}(\bar{q}\lambda_{k}q)^{2}\\!+\\!(\bar{q}i\gamma_{5}\lambda_{k}q)^{2}\big{]}$
(11)
$\displaystyle-K\big{[}\texttt{det}_{f}(\bar{q}(1+\gamma_{5})q)+\texttt{det}_{f}(\bar{q}(1-\gamma_{5})q)\big{]}$
$\displaystyle-U(\Phi[A],\bar{\Phi}[A],T),$

where $q$ denotes the quark fields with three flavors, $u,\ d$, and $s$;
$\hat{m}_{0}=\texttt{diag}(m_{u},\ m_{d},\ m_{s})$ in flavor space; $G$ and
$K$ are the four-point and six-point interacting constants, respectively. The
$\hat{\mu}=diag(\mu_{u},\mu_{d},\mu_{s})$ are the quark chemical potentials.

The covariant derivative in the Lagrangian is defined as
$D_{\mu}=\partial_{\mu}-iA_{\mu}$. The gluon background field
$A_{\mu}=\delta_{\mu}^{0}A_{0}$ is supposed to be homogeneous and static, with
$A_{0}=g\mathcal{A}_{0}^{\alpha}\frac{\lambda^{\alpha}}{2}$, where
$\frac{\lambda^{\alpha}}{2}$ is $SU(3)$ color generators. The effective
potential $U(\Phi[A],\bar{\Phi}[A],T)$ is expressed with the traced Polyakov
loop $\Phi=(\mathrm{Tr}_{c}L)/N_{C}$ and its conjugate
$\bar{\Phi}=(\mathrm{Tr}_{c}L^{\dagger})/N_{C}$. The Polyakov loop $L$ is a
matrix in color space

$L(\vec{x})=\mathcal{P}exp\bigg{[}i\int_{0}^{\beta}d\tau
A_{4}(\vec{x},\tau)\bigg{]},$ (12)

where $\beta=1/T$ is the inverse of temperature and $A_{4}=iA_{0}$.

The Polyakov-loop effective potential is

$\displaystyle\frac{U(\Phi,\bar{\Phi},T)}{T^{4}}$ $\displaystyle=$
$\displaystyle-\frac{a(T)}{2}\bar{\Phi}\Phi+b(T)\mathrm{ln}\big{[}1-6\bar{\Phi}\Phi+4(\bar{\Phi}^{3}+\Phi^{3})-3(\bar{\Phi}\Phi)^{2}\big{]},$
(13)

where

$\\!a(T)\\!=\\!a_{0}\\!+\\!a_{1}\big{(}\frac{T_{0}}{T}\big{)}\\!+\\!a_{2}\big{(}\frac{T_{0}}{T}\big{)}^{2}\,\,\,\texttt{and}\,\,\,\,\,b(T)\\!=\\!b_{3}\big{(}\frac{T_{0}}{T}\big{)}^{3}.$
(14)

The parameters $a_{i}$, $b_{i}$ listed in Table. 1 are fitted according to the
lattice simulation of QCD thermodynamics in pure gauge sector. The $T_{0}=210$
MeV is implemented in the calculation.

Table 1: Parameters in the Polyakov-loop potential Robner07  $a_{0}$ | $a_{1}$ | $a_{2}$ | $b_{3}$
---|---|---|---
$3.51$ | -2.47 | 15.2 | -1.75

The constituent quark mass in the mean field approximation can be derived as

$M_{i}=m_{i}-4G\phi_{i}+2K\phi_{j}\phi_{k}\ \ \ \ \ \ (i\neq j\neq k),$ (15)

where $\phi_{i}$ stands for quark condensate of the flavor $i$.

The thermodynamical potential of bulk quark matter is derived as

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle-2T\sum_{i=u,d,s}\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}(\mathcal{Q}_{1}+\mathcal{Q}_{2})-2\int_{\Lambda}\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}3(E_{u}+E_{d}+E_{s})$
(16)
$\displaystyle+2G\left({\phi_{u}}^{2}+{\phi_{d}}^{2}+{\phi_{s}}^{2}\right)-4K\phi_{u}\,\phi_{d}\,\phi_{s}+U(\bar{\Phi},\Phi,T)$

where $\mathcal{Q}_{1}=\mathrm{ln}(1+3\Phi
e^{-(E_{i}-\mu_{i})/T}+3\bar{\Phi}e^{-2(E_{i}-\mu_{i})/T}+e^{-3(E_{i}-\mu_{i})/T})$,
$\mathcal{Q}_{2}=\mathrm{ln}(1+3\bar{\Phi}e^{-(E_{i}+\mu_{i})/T}+3\Phi
e^{-2(E_{i}+\mu_{i})/T}+e^{-3(E_{i}+\mu_{i})/T})$, and
$E_{i}=\sqrt{\vec{p}^{\,2}+M_{i}^{2}}$ is the dispersion relation.
$\mu_{i}=\mu_{B}/3$ is taken for $u,d,s$ quark flavors. The pressure $p$ and
energy density $\epsilon$ can be derived using the thermodynamic relations in
the grand canonical ensemble as

$P=-\Omega,\,\,\,\,\,\,\,\,\epsilon=-P+Ts+\sum\mu_{i}\rho_{i},$ (17)

where $s$ is the entropy density and $\rho_{i}$ the quark number density of
flavor $i$.

For given $T$ and baryon density $\rho_{B}$, the values of
$\phi_{u},\phi_{d},\phi_{s},\Phi$, $\bar{\Phi}$ and $\mu_{B}$ are determined
by solving the equations by minimizing the thermodynamical potential

$\frac{\partial\Omega}{\partial\phi_{u}}=\frac{\partial\Omega}{\partial\phi_{d}}=\frac{\partial\Omega}{\partial\phi_{s}}=\frac{\partial\Omega}{\partial\Phi}=\frac{\partial\Omega}{\partial\bar{\Phi}}=0,$
(18)

and the relevant constraint condition. Other physical quantities can be then
derived using thermodynamic relations. The numerical results of speed of sound
under different conditions can be then derived according to Eqs. (7)-(10).

In the numerical calculation, a cut-off $\Lambda$ is implemented in 3-momentum
space for divergent integrations. We take the model parameters obtained in
Rehberg96 : $\Lambda=602.3$ MeV, $G\Lambda^{2}=1.835$, $K\Lambda^{5}=12.36$,
$m_{u,d}=5.5$ and $m_{s}=140.7$ MeV, determined by fitting $f_{\pi}=92.4$ MeV,
$M_{\pi}=135.0$ MeV, $m_{K}=497.7$ MeV and $m_{\eta}=957.8$ MeV.

##  3 Numerical results and discussions

In this section, we present the numerical results of the speed of sound under
different constraint conditions and discuss the relations with the QCD phase
transitions.

###  3.1 Sound velocity at constant $s/\rho_{B}$

Firstly, we plot the QCD phase diagram, including the first-order phase
transition (black solid line) and the spinodal structure (blue dashed line),
which separate the phase diagram into the stable, metastable and unstable
phases. The first-order phase transition line is obtained according to the
thermodynamic conditions for two-phase equilibrium, i.e., $T_{1}=T_{2}$,
$\mu_{1}=\mu_{2}$, $P_{1}=P_{2}$ for two stable phases. The spinodal line is
derived with the mechanical unstable condition. The corresponding inflection
points of pressure as a function of density can be determined for a given $T$.
For more details to derive the phase boundaries, one can refer to Refs.
Costa10 ; Shao2018 . The spinodal phase decomposition plays a dominant role in
the experimental exploration of the first-order nuclear liquid-gas
transitionChomaz2004 ; Shao20202 . It has inspired the anticipation to
identify the first-order chiral transition in high-energy heavy-ion collisions
through the spinodal phase separation Mishustin1999 ; Randrup2004 ; Koch2005 ;
Sasaki2007 ; Sasaki2008 ; Randrup2009 ; Steinheimer2012 ; Li2016 ;
Steinheimer2016 ; Steinheimer2017 . The recent simulation suggests that the
spinodal instability can be triggered within a certain energy range
Steinheimer2016 ; Steinheimer2017  . We also demonstrate in Fig. 1 the
isentropic curves with $s/\rho_{B}=0.1,1,3,5,5.9,10,50,100,300$ in the
$T-\rho_{B}$ plane to indicate the evolutionary trajectories of an ideal fluid
at different collision energies.

Figure 1:  QCD phase diagram and the isentropic curves for
$s/\rho_{B}=0.1,1,3,5,5.9,10,50,100,300$ in $T-\rho_{B}$ plane.

Figure 2:  Equations of state of QCD matter for $s/\rho_{B}=0.1,1,3,5,5.9$.
The squares on each curve are the inflection points where $(\frac{\partial
p}{\partial\epsilon})_{s/\rho_{B}}$ changes the sign. The red dashed line is
the profile of these inflection points.

We present in Fig. 2 the equations of state with chosen parameter
$s/\rho_{B}=0.1,1,3,5,5.9$ that pass through the first-order phase transition.
For each curve with $s/\rho_{B}<5.9$, there are two inflection points where
$\partial p/\partial\epsilon$ changes the sign. The red dashed curve in Fig. 2
is the connections of these inflection points, which lies in the spinodal
boundary associated with the first-order phase transition, i.e., in the
interior of the unstable phase. We will present the relation clearly in the
contour map of speed of sound in the $T-\rho_{B}$ panel soon. The inflection
points also correspond to the locations where the sound velocity vanishes in
the phase diagram.

The square of speed of sound can be directly derived with the definition given
in Eq. (2). In Fig. 3, we plot the curve of squared sound speed
$c^{2}_{s/\rho_{B}}$ as functions of energy density along the evolutionary
trajectories for $s/\rho_{B}=0.1,1,3,5,5.9$. Figure 3 shows that there exists
one peak and two valleys on each curve, which indicates that the speed of
sound are closely related to temperature and density. In particular, for the
case of $c^{2}_{s/\rho_{B}}<0$ in Fig. 3, it corresponds to $\partial
p/\partial\epsilon<0$ for a fixed $s/\rho_{B}$, as shown in Fig. 2.

Figure 3:  Values of $c^{2}_{s/\rho_{B}}$ as functions of energy density for
$s/\rho_{B}=0.1,1,3,5,5.9$.

To show more clearly the relation between the speed of sound and QCD phase
transitions, we present the contour map and 3D map of $c^{2}_{s/\rho_{B}}$ as
functions of $T$ and $\rho_{B}$ in Fig. 4 and Fig. 5, respectively. The phase
structure including the chiral crossover, chiral first-order, spinodal and
deconfinement phase transitions. The chiral crossover line and deconfinement
line are determined by requiring $\partial\phi_{i}/\partial T$ and
$\partial\Phi_{i}/\partial T$ taking extreme values for a given chemical
potential.

The two figures indicate that the region around the peaks of
$c^{2}_{s/\rho_{B}}$ in Fig. 3 is located in the region where the chrial
symmetry of $u,d$ quark is approximately restored already but still confined.
The valley at high energy density side in Fig. 3 lies in the region where the
chiral condensate of strange quark changes quickly. After the chiral
restoration of strange quark the sound speed increases again towards high
density. The valley in the low density side is located in the spinodal region
of the first-order phase transition, which is closely related to the chiral
condensate of $u,d$ quark. For each value of ${s/\rho_{B}}$ smaller than 5.9,
there exist a range of $(\frac{\partial p}{\partial\epsilon})_{s/\rho_{B}}$
taking negative values. The red dashed line is the boundary of vanishing sound
velocity with $c^{2}_{s/\rho_{B}}=0$ derived with Eq. (7) .

Figure 4:  Contour map of $c^{2}_{s/\rho_{B}}$ in the $T-\rho_{B}$ plane. The
red dashed line is the boundary of vanishing sound velocity .

Figure 5:  3D map of $c^{2}_{s/\rho_{B}}$ as functions of tempeature and
density. The black solid line is the first-order chiral phase transition line.
The blue dashed line is the spinodal line. The red dashed line is the boundary
of vanishing sound velocity.

Figure 4 and 5 also show that, at the high temperature or very high density,
with the restoration of chiral symmetry $c^{2}_{s/\rho_{B}}$ approaches to
$1/3$, the value of noninteracting gas. The value of $c^{2}_{s/\rho_{B}}$ at
lower density descends with the decrease of temperature. A rapid decrease
occurs in the chiral crossover region of $u,d$ quark, as shown in Fig. 5. It
indicates that the value of speed of sound is sensitive to the change of
dynamical quark mass. In the low-density region, a minimum of
$c^{2}_{s/\rho_{B}}$ appears near the deconfinement phase transition for a
given density (chemical potential). A similar behavior exists in lattice QCD
at zero chemical potential Aoki06 ; Borsanyi14 ; Bazavov14 ; Philipsen13 ;
Borsanyi20 . However, such a feature does not appear in the NJL model which
cannot describe the confinement-deconfinement phase transition Ghosh06 ;
Marty13 ; Deb16 ; Saha18 , which indicates that the color confinement also
plays an important role on the speed of sound near the crossover phase
transition line.

The value of $c^{2}_{s/\rho_{B}}$ is relatively smaller in the region of low
temperature and density. The red dashed line is the boundary of vanishing
sound velocity. Inside this boundary $c^{2}_{s/\rho_{B}}<0$, it means
physically that $(\frac{\partial p}{\partial\epsilon})_{s/\rho_{B}}<0$. In
this region, the mechanically stable condition is broken and the corresponding
sound wave equation becomes a decay function. A perturbance can not be
propagated like a sound wave in this situation. It can be seen that such a
region lies in the interior of the unstable phase of the spinodal structure.
Figure 4 and 5 also indicate that the speed of sound at the critical endpoint
is small but not zero in the mean field approximation.

The value of $c^{2}_{s/\rho_{B}}$ reflects the speed of sound in an ideal
fluid which can be approximately realized in heavy-ion collision experiments.
On the other hand, $s/\rho_{B}$ is connected with the collision energy. If the
value of $c^{2}_{s/\rho_{B}}$ at a fixed energy can be extracted from the
charged particle multiplicity $\langle dN_{ch}/d\eta\rangle$ Gardim20 ; Sahu21
; Biswas20 , we can access the information of phase transition using the
relation between $c^{2}_{s/\rho_{B}}$ and QCD phase diagram. Furthermore,
combining with the beam energy scan experiments, it provides a possible way to
diagnose the QCD phase structure. It is also inspiring for study on the
gravitational signal from the cosmic QCD phase transition in which a constant
speed of sound $\sqrt{1/3}$ is usually taken in literature.

###  3.2 Sound velocity at constant $\rho_{B}$ and $s$

We present in Fig. 6 the contour map of $c_{\rho_{B}}^{2}$ in the $T-\rho_{B}$
panel. Besides at the high-temperature side, this figure shows that
$c_{\rho_{B}}^{2}$ take relatively larger values in the region of low
temperature and density. The value is even larger than $1/3$, in particular,
in the metastable phase and unstable phase. There also exists a wide region
(inside the red line filled with the blue color) of $c_{\rho_{B}}^{2}<0$. More
physically, it means that the $(\frac{\partial
p}{\partial\epsilon})_{\rho_{B}}<0$ in this region, i.e., the pressure
decreases with the increase of energy density along the line of constant
density. The red line shows the boundary of vanishing sound speed at constant
baryon density.

Figure 6:  Contour map of $c^{2}_{\rho_{B}}$ in the $T-\rho_{B}$ plane. The
red dashed line is the boundary of vanishing sound velocity.

Figure 7:  Contour map of $c^{2}_{s}$ in the $T-\rho_{B}$ plane. The red
dashed line is the boundary of vanishing sound velocity. The green curves
shows the paths of $s=0.1,0.5,1.0,5$.

The contour map of $c_{s}^{2}$ is presented in Fig. 7. The behavior of
$c_{s}^{2}$ at low temperature and high density is similar with that of
$c^{2}_{s/\rho_{B}}$, because the curves at constant $s$ and $s/\rho_{B}$ in
the $T-\rho_{B}$ diagram are both roughly parallel with the density axis, as
indicated in Fig. 1. and Fig. 7. However, the curves at constant $s$ at low
density are almost perpendicular to those at constant $s/\rho_{B}$, the
resulting behaviors of $c_{s}^{2}$ and $c_{s/\rho_{B}}^{2}$ are quite
different in the corresponding region. A distinct characteristic is the
location of vanishing sound speed. The red dashed line in Fig. 7 is the
boundary of $c_{s}^{2}=0$. $(\frac{\partial p}{\partial\epsilon})_{s}$ takes
minus values inside the boundary, which includes a wide range of the first-
order phase transition and a region around the CEP.

Figure 8:  Regions of $(\frac{\partial\mu_{B}}{\partial T})_{s/\rho_{B}}<0$ in
the $T-\rho_{B}$ panel.

The region of $c_{s}^{2}<0$ in Fig. 7 and $c_{\rho_{B}}^{2}<0$ in Fig. 6 are
connected with the formula

$\bigg{(}\frac{\partial\mu_{B}}{\partial
T}\bigg{)}_{s/\rho_{B}}=\frac{\mu_{B}}{T}\frac{(\frac{\partial
p}{\partial\epsilon})_{s}}{(\frac{\partial
p}{\partial\epsilon})_{\rho_{B}}}=\frac{\mu_{B}}{T}\frac{c^{2}_{s}}{c^{2}_{\rho_{B}}}.$
(19)

When the condition $(\frac{\partial\mu_{B}}{\partial T})_{s/\rho_{B}}<0$ is
fulfilled, one of the two physical quantities $c_{s}^{2}$ and
$c_{\rho_{B}}^{2}$ takes a negative value. We show in Fig. 8 the regions of
$(\frac{\partial\mu_{B}}{\partial T})_{s/\rho_{B}}<0$. The numerical results
indicate that there indeed exists the regions where
$(\frac{\partial\mu_{B}}{\partial T})_{s/\rho_{B}}<0$ in the $T-\rho_{B}$
diagram. Comparing Fig. 8 with the negative value regions in Fig. 6 and Fig.
7, we can conclude that these numerical results confirm the formula in Eqs.
(19).

The behaviors of $c_{\rho_{B}}^{2}$ and $c_{s}^{2}$ in the phase diagram can
be used to study the fluid properties of quark gluon plasma. The values of
$c_{\rho_{B}}^{2}$ and $c_{s}^{2}$ are important parameters to indicate the
intermediated process in the evolution of a fluid. Eqs. (4) and (5) clearly
show that $c_{\rho_{B}}^{2}$ and $c_{s}^{2}$ are connected with the temporal
derivatives of temperature and chemical potential, respectively. Moreover,
$c_{\rho_{B}}^{2}$ and $c_{s}^{2}$ are related to the bulk viscosity of a
fluid. In particular they directly connect with the bulk viscosity
coefficient. Exploring the relation between the bulk viscosity and phase
transition is attractive to study the dissipation in the evolution of QGP. A
further research in this respect is undergoing.

###  3.3 Sound velocity at constant $T$ and $\mu_{B}$

The contour maps of $c^{2}_{T}$ at constant temperature and $c^{2}_{\mu_{B}}$
at constant chemical potential in the $T-\rho_{B}$ panel are demonstrated in
Fig. 9 and Fig. 10, respectively. The two figures show that $c^{2}_{T}$ and
$c^{2}_{\mu_{B}}$ are both close to $1/3$ at high temperature.

The contour of $c^{2}_{T}$ looks in general like that of $c^{2}_{s/\rho_{B}}$,
because the curves of constant $s/\rho_{B}$ in the $T-\rho_{B}$ panel are
almost parallel to the density axis in a wide range. The relative larger
deviation lies in the range with densities smaller than the boundary of the
first-order transition on the low-density side. The deviation produces
different behaviors between $(\frac{\partial p}{\partial\epsilon})_{T}$ and
$(\frac{\partial p}{\partial\epsilon})_{s/\rho_{B}}$ at low density. A crutial
point is that the inflections points of $(\frac{\partial
p}{\partial\epsilon})_{T}$ and $(\frac{\partial
p}{\partial\epsilon})_{s/\rho_{B}}$ are different, i.e., the boundary of
vanishing sound speed are different for the two cases.

Figure 9: Contour map of $c^{2}_{T}$ in the $T-\rho_{B}$ plane. The red dashed
line is the boundary of vanishing sound velocity.

Figure 10:  Contour map of $c^{2}_{\mu_{B}}$ in the $T-\rho_{B}$ plane. The
red dashed line is the boundary of vanishing sound velocity.

Both the figure 9 and 10 show that the boundary of zero sound speed at
constant temperature or chemical potential is just the spinodal line
associated with the first-order phase transition. However the boundary of zero
sound velocity at constant $s/\rho_{B}$ is in the interior of the spinodal
structure as shown in Fig. 4. The negative values of $c^{2}_{T}$ and
$c^{2}_{\mu_{B}}$ both appear in the unstable phase of the spinodal structure,
since $(\frac{\partial p}{\partial\epsilon})_{T}$ and $(\frac{\partial
p}{\partial\epsilon})_{\mu_{B}}$ are negative in this region.

From the behavior of $c^{2}_{T}$ derived above, we can approximately deduce
the speed of sound in the quark core of a massive neutron star. Since the
richness of lepton including electron and muon approaches to zero at high
density of a hybrid neutron star, and the richness of $u,d,s$ quark tends to
be equivalent Shao2013 , quite similar to the situation in this study at high
density and low temperature. Therefore, it may be concluded to a certain
degree that the squared speed of sound in the quark core of a massive neutron
star gradually approaches to $1/3$ at high density. A further study on the
speed of sound in neutron star matter with a hadron-quark phase transition
will be conducted with the combination of observation data.

##  4 Summary

In this work, we studied the speed of sound in QCD matter at finite
temperature and density in the PNJL model. We derived the behavior of speed of
sound under different definitions in the $T-\rho_{B}$ phase diagram including
the stable phase, metastable and unstable phases associated with the first-
order phase transition. We systematically discussed the relations between the
speed of sound and QCD phase structure.

The numerical results indicate that the squared speed of sounds under
different definitions are all approaching to 1/3 at high temperature. However,
the behaviors in the phase-transition region are closely related to the non-
perturbative interaction and the phase structure. From the perspective of idea
fluid evolution, more attention are put on the speed of sound under isentropic
condition. The calculation indicates that $c^{2}_{s/\rho_{B}}$ is nonzero at
the CEP under the mean field approximation, and the boundary of vanishing
sound velocity is further derived.

We also obtained the contour maps of $c^{2}_{X}$ $(X=s,\rho_{B},T,$
$\mu_{B})$, and analyzed their relations with the QCD phase transition, as
well as the relations between different definitions. For each definition of
speed of sound, we find that there exist one or several regions with
$(\frac{\partial p}{\partial\epsilon})_{X}<0$, even in the stable phase for
$X=s$ and $\rho_{B}$.

The different definitions of speed of sound are involved in some important
physics problems in nuclear physics, such as the evolution of quark gluon
plasma, the cosmic QCD phase transition, the bulk viscosity of strongly
interacting matter and the equation of state of a hybrid neutron star. The
indepth investigations on these physics issues will be performed in the
future.

Acknowledgments This work is supported by the National Natural Science
Foundation of China under Grant No. 11875213.

Data Availability Data sharing not applicable to this article as no datasets
were generated or analyzed during the current study.

## References

  * (1) H. C. Song, S. A. Bass, U. Heinz, T. Hirano, and C. Shen, 106, 192301 (2011).
  * (2) H. C. Song, S. A. Bass, U. Heinz, Phys. Rev. C 83, 024912 (2011).
  * (3) P. Deb, G. P. Kadam, and H. Mishra, Phys. Rev. D 94, 094002 (2016).
  * (4) F. G. Gardim, G. Giacalone, M. Luzum, and J. Y. Ollitrault, Nat. Phys. 16, 615 (2020).
  * (5) D. Sahu, S. Tripathy, R. Sahoo, and A. R. Dash, Eur. Phys. J. A 56, 187 (2021).
  * (6) D. Biswas, K. Deka, A. Jaiswal, and S. Roy, Phys. Rev. C 102, 014912 (2020).
  * (7) A. Sorensen, D. Oliinychenko, V. Koch, and L. McLerran, Phys. Rev. Lett. 127, 042303 (2021).
  * (8) B. Reed and C. J. Horowitz, Phys. Rev. C 101, 045803 (2020).
  * (9) A. Kanakis-Pegios, P. S. Koliogiannis, and Ch. C. Moustakidis, Phys. Rev. C 102, 055801 (2020).
  * (10) S. Han and M. Prakash, Astrophys. J. 899, 164 (2020).
  * (11) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018).
  * (12) S. K. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, and A. L. Watts, Mon. Not. R. Astron. Soc. 485, 5363 (2019).
  * (13) M. M. Forbes, S. Bose, S. Reddy, D. Zhou, A. Mukherjee, and S. De, Phys. Rev. D 100, 083010 (2019).
  * (14) C. Drischler, S. Han, J. M. Lattimer, M. Prakash, S. Reddy, and T. Zhao, Phys. Rev. C 103, 045808 (2021).
  * (15) R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Phys. Rev. C 102, 055803 (2020).
  * (16) S. Han, M. A. A. Mamun, S. Lalit, C. Constantinou, and M. Prakash, Phys. Rev. D 100, 103022 (2019).
  * (17) T. Kojo, AAPPS Bull. 31, 11 (2021).
  * (18) P. Jaikumar, A. Semposki, M. Prakash, and C. Constantinou, Phys. Rev. D 103, 123009 (2021).
  * (19) O. Philipsen, Prog. Part. Nucl. Phys. 70, 55 (2013).
  * (20) S. Borsányi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 125, 052001 (2020).
  * (21) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, K. K. Szabo, Nature (London) 443, 675 (2006).
  * (22) S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Sabzó, Phys. Lett. B 730, 99 (2014).
  * (23) A. Bazavov et al. (hotQCD Collaboration), Phys. Rev. D. 90, 094503 (2014).
  * (24) M. Motta, R. Stiele, W. M. Alberico, A. Beraudo, Eur. Phys. J. C 80, 770 (2020).
  * (25) S. K. Ghosh, T. K. Mukherjee, M. G. Mustafa, and R. Ray, Phys. Rev. D 73, 114007 (2006).
  * (26) R. Marty, E. Bratkovskaya, W. Cassing, J. Aichelin, and H. Berrehrah, Phys. Rev. C 88, 045204 (2013).
  * (27) K. Saha, S. Ghosh, S. Upadhaya, S. Maity, Phys. Rev. D 97, 116020 (2018).
  * (28) Y. P. Zhao, Phys. Rev. D 101, 096006, (2020).
  * (29) B. J. Schaefer, M. Wagner, and J. Wambach, Phys. Rev. D 81, 074013 (2010).
  * (30) A. Abhishek, H. Mishra, and S. Ghosh, Phys. Rev. D 97, 014005 (2018).
  * (31) R. Venugopalan and M. Prakash, Nucl. Phys. A546, 718 (1992).
  * (32) M. Bluhm, P. Alba, W. Alberico, A. Beraudo, and C. Ratti, Nucl. Phys. A929, 157 (2014).
  * (33) Z. V. Khaidukov, M. S. Lukashov, and Yu. A. Simonov, Phys. Rev. D 98, 074031 (2018).
  * (34) Z. V. Khaidukov and Yu. A. Simonov, Phys. Rev. D 100, 076009 (2019).
  * (35) V. Mykhaylova and C. Sasaki, Phys. Rev. D 103, 014007 (2021).
  * (36) W. B. He, G. Y. Shao, X. Y. Gao, X. R. Yang, and C. L. Xie, Phys. Rev. D 105, 094024 (2022).
  * (37) I. N. Mishustin, Phys. Rev. Lett. 82, 4779 (1999).
  * (38) J. Randrup, Phys. Rev. Lett. 92, 122301 (2004).
  * (39) V. Koch, A. Majumder, and J. Randrup, Phys. Rev. C 72, 064903 (2005).
  * (40) C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. Lett. 99, 232301 (2007).
  * (41) C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. D 77, 034024 (2008).
  * (42) J. Randrup, Phys. Rev. C 79, 054911 (2009); Phys. Rev. C 82, 034902 (2010).
  * (43) J. Steinheimer and J. Randrup, Phys. Rev. Lett. 109, 212301 (2012).
  * (44) J. Steinheimer and J. Randrup, Eur. Phys. J. A 52, 239 (2016).
  * (45) F. Li and C. M. Ko, Phys. Rev. C 93, 035205 (2016).
  * (46) J. Steinheimer and V. Koch, Phys. Rev. C. 96, 034907 (2017).
  * (47) G. Y. Shao, X. Y. Gao, and W. B. He, Eur. Phys. J. A 56, 115 (2020).
  * (48) M. Ahmadvand, K. Bitaghsir Fadafan, Phys. Lett. B 779, 1 (2018).
  * (49) Tuomas V. I. Tenkanen and Jorinde van de Vis, J. High Energ. Phys. 2022, 302 (2022).
  * (50) S. Rößner, C. Ratti, and W. Weise, Phys. Rev. D 75, 034007 (2007).
  * (51) P. Rehberg, S. P. Klevansky, and J. Hüfner, Phys. Rev. C 53, 410 (1996).
  * (52) P. Costa, M. C. Ruivo, C. A. de Sousa, and H. Hansen, Symmetry 2, 1338 (2010).
  * (53) G. Y. Shao, Z. D. Tang, X. Y. Gao, and W. B. He, Eur. Phys. J . C 78, 138 (2018).
  * (54) P. Chomaz, M. Colonna, and J. Randrup, Phys. Rep. 389 263 (2004).
  * (55) G. Y. Shao, X. Y. Gao, and W. B. He, Phys.Rev.D 101, 074029 (2020).
  * (56) G. Y. Shao, M. Colonna, M. Di Toro, Y. X. Liu, and B. Liu, Phys.Rev.D 87, 096012 (2013).

##  Appendix A Derivations of the formulae of speed of sound under different
definitions in the temperature and density space

The general definition of speed of sound is

$c_{X}^{2}=\left(\frac{\partial p}{\partial\epsilon}\right)_{X},$ (20)

where $X$ is a physics quantum fixed in the calculation of sound speed. In
practice, the squared speed of sound $c^{2}_{X}$
$(X=s/\rho_{B},s,\rho_{B},T,\mu_{B})$ under different conditions are taken in
dealing with different physics problems. With the basic definitions of speed
of sound, calculation can only be done along some special paths.

To calculate the speed of sound under different definitions in the whole
$T-\rho_{B}$ space, it is necessary to derive the corresponding formulae in
terms of temperature and density. Using the Jacobian formula in
thermodynamics, we can derive

$\displaystyle c_{X}^{2}\left(T,\rho_{B}\right)$ $\displaystyle=$
$\displaystyle\left(\frac{\partial
p}{\partial\epsilon}\right)_{X}=\frac{\partial(p,X)}{\partial(\epsilon,X)}=\frac{\frac{\partial(p,X)}{\partial\left(T,\rho_{B}\right)}}{\frac{\partial(\epsilon,X)}{\partial\left(T,\rho_{B}\right)}}$
(21) $\displaystyle=$ $\displaystyle\frac{\left(\frac{\partial p}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial
X}{\partial\rho_{B}}\right)_{T}-\left(\frac{\partial
p}{\partial\rho_{B}}\right)_{T}\left(\frac{\partial X}{\partial
T}\right)_{\rho_{B}}}{\left(\frac{\partial\epsilon}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial
X}{\partial\rho_{B}}\right)_{T}-\left(\frac{\partial\epsilon}{\partial\rho_{B}}\right)_{T}\left(\frac{\partial
X}{\partial T}\right)_{\rho_{B}}}$

According to the thermodynamic characteristic function in the giant canonical
ensemble, it is convenient to get the following relations for isospin
symmetric matter

$\left(\frac{\partial p}{\partial
T}\right)_{\rho_{B}}=s+\rho_{B}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}},$ (22) $\left(\frac{\partial
p}{\partial\rho_{B}}\right)_{T}=\rho_{B}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T},$
(23) $\left(\frac{\partial\epsilon}{\partial
T}\right)_{\rho_{B}}=T\left(\frac{\partial s}{\partial T}\right)_{\rho_{B}},$
(24)
$\left(\frac{\partial\epsilon}{\partial\rho_{B}}\right)_{T}=T\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}+\mu_{B},$ (25)
$\left(\frac{\partial{(s/\rho_{B})}}{\partial
T}\right)_{\rho_{B}}=\frac{1}{\rho_{B}}\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}},$ (26)

and

$\left(\frac{\partial{(s/\rho_{B})}}{\partial\rho_{B}}\right)_{T}=\frac{1}{\rho_{B}}\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}-\frac{s}{\rho_{B}^{2}}$ (27)

For the different constraint conditions, $X=s/\rho_{B},s,\rho_{B},T,\mu_{B}$,
we can derived the corresponding formulae of speed of sound as follows

$\\!c^{2}_{s/\rho_{B}}\\!=\\!\frac{\\!s^{2}\\!+\\!\rho_{B}^{2}\\!\left[\\!\left(\frac{\\!\partial\mu_{B}}{\\!\partial\rho_{B}}\\!\right)_{\\!T}\\!\left(\frac{\partial
s}{\\!\partial
T}\\!\right)_{\\!\rho_{\\!B}}\\!-\\!\left(\frac{\partial\mu_{\\!B}}{\partial
T}\\!\right)_{\\!\rho_{\\!B}}\\!\left(\frac{\partial
s}{\partial\rho_{\\!B}}\\!\right)_{\\!T}\\!\right]\\!+\\!s\rho_{\\!B}\\!\left[\\!\left(\frac{\partial\mu_{\\!B}}{\partial
T}\\!\right)_{\\!\rho_{\\!B}}\\!-\\!\left(\frac{\partial
s}{\partial\rho_{\\!B}}\\!\right)_{\\!T}\right]}{\left(Ts+\mu_{B}\rho_{B}\right)\left(\frac{\partial
s}{\partial T}\right)_{\rho_{B}}},$ (28)
$c_{s}^{2}=\frac{\rho_{B}\left[\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}-\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}\right]-s\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}}{\mu_{B}\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}},$ (29)
$c_{\rho_{B}}^{2}=\frac{s+\rho_{B}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}}{T\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}},$ (30)
$c_{T}^{2}=\frac{\rho_{B}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}}{T\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}+\mu_{B}},$ (31)

and

$c_{\mu_{B}}^{2}=\frac{s\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}}{T\left[\left(\frac{\partial
s}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial\mu_{B}}{\partial\rho_{B}}\right)_{T}-\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}\left(\frac{\partial
s}{\partial\rho_{B}}\right)_{T}\right]-\mu_{B}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{\rho_{B}}}.$ (32)


# Pb10-xCux(PO4)6O: a Mott or charge transfer insulator
in need of further doping for (super)conductivity

Liang Sia,b, Markus Wallerbergerb, Andriy Smolyanyukb, Simone di Cataldob, Jan
M. Tomczakb,c, and Karsten Heldb a School of Physics, Northwest University,
Xi’an 710127, China
b Institute of Solid State Physics, TU Wien, 1040 Vienna, Austria
c Department of Physics, King’s College London, Strand, London WC2R 2LS,
United Kingdom

###### Abstract

We briefly review the status quo of research on the putative superconductor
Pb9Cu(PO4)6O also known as LK-99. Further, we provide ab initio derived tight-
binding parameters for a two- and five-band model, and solve these in
dynamical-mean-field theory. The ratio interaction-to-bandwidth makes LK-99 a
Mott or charge transfer insulator. Electron or hole doping (which is different
from substituting Pb by Cu and thus differs from LK-99) is required to make it
metallic and potentially superconducting.

##  I  Introduction

In recent preprints Lee  _et al._ (2023a, b), Lee, Kim, et al. reported the
discovery of a room-temperature superconductor at ambient-pressure:
Pb10-xCux(PO4)6O with $0.9<x<1.1$. They had previously named this material
LK-99 after their initials and the year of the first synthesis. Their more
recent samples show somewhat stronger signatures of superconductivity Lee  _et
al._ (2023c, a, b): (i) a sharp drop in the resistivity Lee  _et al._ (2023c,
b), according to Ref. Lee  _et al._ , 2023a to the order of
$10^{-10}-10^{-11}\,\Omega$cm though in Refs. Lee  _et al._ , 2023c, b a
higher noise level is visible, (ii) a negative (diamagnetic) spin
susceptibility and levitation on a magnet Lee  _et al._ (2023b), and (iii)
sharp voltage jumps at critical currents, with the critical currents vanishing
in approximately a quarter-circle as a function of temperature and magnetic
field Lee  _et al._ (2023c, a, a).

If LK-99 is truly a superconductor at ambient temperature and pressure, it is
arguably one of the most significant physics discoveries of recent history.
However, experimental confirmation is urgently needed: The above experiments,
while indicative of superconductivity, do not unambiguously prove it. (i) The
noise level of the resistivity appears too large for concluding that LK-99 has
zero resistance. (ii) The negative susceptibility and levitation can be caused
by a simple diamagnet. (iii) The voltage jumps might also be caused by contact
issues. One has to admit however that taken everything together, the overall
picture provides quite some indication for superconductivity. If the critical
temperature was 1 K [and if taking into account the first confirmations of (i)
and (ii), see below], the scientific community would now most likely be quite
positive that at least parts of the LK-99 sample are superconducting. But room
temperature superconductivity is an extraordinary claim, and extraordinary
claims rightfully require an extraordinarily solid proof. Such waterproof
evidence has not been given as of the time of submitting this article, neither
has solid evidence against.

Naturally, the results by Lee, Kim et al. led to huge experimental and
theoretical efforts. Let us briefly review the status quo of these subsequent
works as of the day of submission (Aug. 8th, 2023): The levitation (i) has
been reproduced by Wu et al. Wu  _et al._ (2023) and further groups on social
media. A sharp drop in resistivity has been confirmed by Hou et al. Hou  _et
al._ (2023), albeit at 100 K instead of 100∘ Celsius. Hou et al. also report
two strange resistivity jumps above 250 K (which the authors suggest might be
caused by issues with the electrode contacts) as well as an abnormal field
dependence.

In contrast, other experimental groups report an opposite behavior. Liu et al.
Liu  _et al._ (2023) find an increase of the resistivity with decreasing
temperature indicating that LK-99 is a semiconductor or insulator. They also
observe a paramagnetic response instead of a diamagnetic one, and the magnetic
susceptibility as well as the resistivity increases with decreasing
temperature. Kumar et al. Kumar  _et al._ (2023, 2023) also successfully
synthesized Pb10-xCux(PO4)6O in a modified lead apatite structure, and report
a diamagnetic insulator. Abramiam et al. Abramian  _et al._ (2023) conjecture
that the samples by Lee et al. are not pure LK-99, but a coexistence of
superconducting and non-superconducting regions; superconductivity might
emerge from another material. Guo et al. Guo  _et al._ (2023) report a
ferromagnetic hysteresis and half levitation.

As for theory, density functional theory (DFT)Hohenberg and Kohn (1964) is
state-of-the-art for calculating crystal structures and for getting, at the
bare minimum, a first idea of the electronic structure. Five groups Lai  _et
al._ (2023); Griffin (2023); Si and Held (2023); Kurleto  _et al._ (2023);
Cabezas-Escares  _et al._ (2023) independently performed such DFT
calculations, appearing on arXiv within days, and showing similar results (cf.
Kumar  _et al._ (2023); Tao  _et al._ (2023); Sun  _et al._ (2023)): For the
lead apatite crystal structure with one Pb atom replaced by Cu, two very flat
bands cross the Fermi energy. Below these are still flat, but slightly more
dispersive O bands, and another Cu band. Some of the DFT calculations also
analyze possible alternative Cu and O positions Si and Held (2023); Griffin
(2023); Cabezas-Escares  _et al._ (2023); Sun  _et al._ (2023) so far to a
very limited extent; Cabezas-EscaresCabezas-Escares  _et al._ (2023) find an
instability with a simplified frozen phonon calculation. Lai et al. Lai  _et
al._ (2023) suggest that gold-doped lead apatite may have stronger effects
than Cu. Griffin Griffin (2023), Si and Held Si and Held (2023), and Kurleta
et al. Kurleto  _et al._ (2023) argue that the flat bands might boost
electron-phonon mediated superconductivity; Si and Held Si and Held (2023)
also suggest purely electronic flat-band superconductivity Kuroki  _et al._
(2005); Iglovikov  _et al._ (2014); Aoki (2020) as a possible alternative.

Some groups Lai  _et al._ (2023); Griffin (2023); Kurleto  _et al._ (2023);
Cabezas-Escares  _et al._ (2023) consider the flat DFT bands crossing the
Fermi energy as evidence that Cu doping $x\approx 1$ makes insulating lead
apatite metallic, thus explaining the conducting and prospectively
superconducting state of LK-99. While this is suggestive from the DFT results,
Si and Held Si and Held (2023) estimate the interaction-to-bandwidth $U/W$ to
be of ${\cal O}(10)$ and thus conclude that LK-99 must be a Mott (or charge
transfer) insulator, see Fig. 1 for an illustration. They further conjecture
that the accompanying spin-1/2 should show a strong paramagnetic response so
that a diamagnet without superconductivity is difficult to imagine.

Figure 1: Top: Schematics of DFT bandstructure for Pb9Cu(PO4)6O. Middle: Mott-
Hubbard splitting of the Cu $d_{xz,yz}$ orbitals leading to a Mott or charge
transfer insulator. Note here we visualize the case with an orbital symmetry
breaking (ordering) indicating by having different orbitals in the lower and
upper Hubbard band. Interrupted lines indicate a larger energy separation.
Bottom: if doped the Mott or charge transfer insulator becomes metallic. here
hole doping is visualized. For electron doping the quasiparticle band (and
Fermi energy) would be closer to the upper Hubbard band. A charge transfer
insulator for electron doping is not plausible, because the next unoccupied
orbitals (Pb $p$) are too high in energy Si and Held (2023).

Indeed such a MottGebhard (1997) or charge transfer Zaanen  _et al._ (1985)
insulator might explain the simultaneous experimental findings of a
paramagnetic insulator. At the same time, the metallic (and prospectively
superconducting) behavior found in the other experiments Lee  _et al._ (2023c,
a, b); Wu  _et al._ (2023); Hou  _et al._ (2023) is possible if (part of) the
sample is doped, see Fig. 1 (bottom). Indeed high-temperature cuprate
superconductors Bednorz and Müller (1986) are in the same class of a charge
transfer (or Mott) insulator Scalapino (2012), and also have Cu and O orbitals
as the relevant ones. However, this is how far the similarity goes.
Superconducting cuprates have $U/W$ only of ${\cal O}(1)$, and the Cu atoms
form a square lattice in the CuO2 planes, while the lead apatite crystal
structure is hexagonal. To make such a Mott or charge transfer insulator
metallic, one needs electron or hole doping which is not possible by changing
$x$, i.e., the ratio of Cu:Pb. The synthesis procedure Lee  _et al._ (2023b)
and lead apatite crystal structure suggests that electron doping was possible
for $0<y\ll 1$ and $z<0$ and hole doping for $z>0$ in
Pb10-xCux(P1-ySyO4)6O1+z. Note, the nominal oxidation states are: Pb2+, Cu2+,
P5+, S6+ and O2-. Anyhow, this is merely an educated guess.

Also based on the picture of a doped Mott insulator, Baskaran Baskaran (2023)
speculates that Cu atoms cluster in chains or 2D patches with a Cu0 electronic
configuration instead of Cu2+ Lai  _et al._ (2023); Griffin (2023); Si and
Held (2023); Kurleto  _et al._ (2023); Cabezas-Escares  _et al._ (2023), and
thus realize his theory of a broad band Mott localization. First groups also
started doing calculations for two-orbital models on a triangular lattice,
using a Bardeen, Cooper and Schrieffer(BCS) Bardeen  _et al._ (1957) type of
coupling Tavakol and Scaffidi (2023) and slave-boson mean field theory. The
authors find $f$-wave and $s$-wave superconductivity (though at too low
critical temperatures) Oh and Zhang (2023), respectively.

In this paper, we would like to put such tight binding parameters on a more
solid basis. To this end, we do a Wannier function projection and calculate
the tight binding parameters for: (i) a two-band low energy model made up of
the Cu $d_{xz/yz}$ orbitals and (ii) a five-band model also involving the O
$p_{x/y}$ states just below these and, additionally, the next Cu ($d_{z^{2}}$)
orbital, see Fig. 1 (top) [this Fig. does not include the 3rd Cu $d$ orbital
below the O $p_{x/y}$ bands, cf. Fig. 2 below]. These tight-binding models can
be used for subsequent many-body calculations and are listed in 1 and 2,
respectively. We motivate the strength of the Kanamori interaction parameters
on the Cu sites. Finally we solve these models in DMFT and find a Mott or
charge transfer insulator. A similar insulator is also obtained in DFT+$U$,
but only if the crystal symmetry lifts the degeneracy of the Cu $d_{xz/yz}$
orbitals.

Let us put some caveats here regarding the low-energy model. It assumes the
periodic continuation of a unit cell with a single formula unit ($x=1$, $y=0$,
$z=0$) and optimized O and Cu positions Si and Held (2023). This yields among
others, a regular triangular lattice of the Cu sites. Other O and Cu positions
are, however, so close in energy Si and Held (2023) that we must expect a
disordered arrangement of these at room temperature — unless there is a
crystal distortion stabilizing some arrangement. The x-ray diffraction (XRD)
patterns Lee  _et al._ (2023b); Liu  _et al._ (2023); Hou  _et al._ (2023)
clearly indicate an undistorted lead apatite structure without periodic
arrangement of the Cu atoms. For the matter of Mott insulator or not, this is
not relevant, but for (super)conductivity the possible long-range ordering of
the Cu and O atoms or vice versa a disordered arrangement of these is very
relevant.

Orbital-1 | Orbital-2 | $t_{000}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{2}$ | $t_{2}$ | $t_{2}$ | $t_{2}$ | $t_{2}$ | $t_{2}$ | $t_{z}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  | (000) | (100) | (010) | (110) | (-1-10) | (-100) | (0-10) | (120) | (210) | (1-10) | (-1-20) | (-2-10) | (-110) | (001)
$e_{g}$(2) |  $e_{g}$(2) | $-24$ | $6$ | $-5$ | $0$ | $0$ | $6$ | $-5$ | 0 | 0 | 0 | 0 | 0 | 0 | $-10$
$e_{g}$(1) |  $e_{g}$(1) | $-24$ | $-6$ | $5$ | $0$ | $0$ | $-6$ | $5$ | 0 | 0 | 0 | 0 | 0 | 0 | $-10$
$e_{g}$(2) |  $e_{g}$(1) | $0$ | $6$ | $4$ | $-2$ | $15$ | $-11$ | $-13$ | 0 | 0 | 0 | 0 | 0 | 0 | $0$
$e_{g}$(1) |  $e_{g}$(2) | $0$ | $-11$ | $-13$ | $15$ | $-2$ | $6$ | $4$ | 0 | 0 | 0 | 0 | 0 | 0 | $0$
Table 1: Two-band model, comprising two Cu $e_{g}$(1) and $e_{g}$(2) orbitals (corresponding to $d_{yz}$, $d_{yz}$). Hopping parameters up to second nearest from Cu to Cu sites and Orbital-1 to Orbital-2. Here, $t_{1}$ and $t_{2}$ denote 1st and 2nd nearest hoppings. The numbers in brackets indicate the hopping vector in real space, see Fig. 3(f). Hopping terms that are close to or smaller than 1 meV are put as 0 meV. All hoppings are in units of meV; $t_{000}$ is the on-site energy.  Orb-1 | Orb-2 | $t_{000}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{1}$ | $t_{z}$
---|---|---|---|---|---|---|---|---|---
|  | (000) | (100) | (010) | (110) | (-1-10) | (-100) | (0-10) | (001)
$e_{g}$(1) |  $e_{g}$(1) | $-77$ | $0$ | $-4$ | $0$ | $0$ | $0$ | $-4$ | $-4$
$e_{g}$(2) |  $e_{g}$(2) | $-77$ | $-4$ | $0$ | $0$ | $0$ | $-4$ | $0$ | $-4$
$e_{g}$(2) |  $e_{g}$(1) | $0$ | $-4$ | $-4$ | $6$ | $0$ | $0$ | $0$ | $2$
$e_{g}$(1) |  $e_{g}$(2) | $0$ | $0$ | $0$ | $0$ | $6$ | $-4$ | $-4$ | $-2$
$e_{g}$(1) | O-$p$(1) | $0$ | $0$ | $-3$ | $0$ | $0$ | $3$ | $0$ | $0$
$e_{g}$(1) | O-$p$(2) | $-7$ | $0$ | $0$ | $0$ | $-4$ | $0$ | $0$ | $0$
$e_{g}$(2) | O $p$(1) | $0$ | $0$ | $6$ | $2$ | $0$ | $5$ | $0$ | $0$
$e_{g}$(2) | O $p$(2) | $0$ | $0$ | $3$ | $2$ | $0$ | $-4$ | $0$ | $0$
O $p$(1) | O $p$(1) | $-366$ | $10$ | $0$ | $-11$ | $-11$ | $10$ | $0$ | $-85$
O $p$(2) | O $p$(2) | $-366$ | $-10$ | $0$ | $11$ | $11$ | $-10$ | $0$ | $-85$
O $p$(2) | O $p$(1) | $0$ | $-4$ | $-22$ | $15$ | $-6$ | $18$ | $0$ | $-3$
O $p$(1) | O $p$(2) | $0$ | $18$ | $0$ | $-6$ | $15$ | $-4$ | $-22$ | $3$
$a_{1g}$ | $a_{1g}$ | $-715$ | $6$ | $6$ | $6$ | $6$ | $6$ | $6$ | $12$
Table 2: Five-band model, comprising besides the two Cu $e_{g}$ orbitals of
the two-band model additionally two O $p$(1) and $p$(2) orbitals
(corresponding to $p_{x}$ and $p_{y}$) as well as one further Cu $a_{1g}$
($d_{z^{2}}$) orbital. The numbers in brackets indicate the hopping vector in
real space, see Fig. 3 (f). Hopping parameters up to 1st nearest neighbor
between all Cu and O sites and orbitals are given. As before, $t_{1}$
indicates 1st nearest hoppings, the numbers in brackets indicate hopping
vectors in real space; hopping terms obviously smaller than 1 meV are
approximated as 0 meV. The energy unit is meV.  Figure 2: (a) DFT-relaxed
structure of Pb9Cu(PO4)6O; (b) schematic figure of energy band splitting from
the octahedral Cu-coordination of the CuO6 motif; Wannier projections for the
two bands model (c) and five bands model (d), respectively. The Wannier bands
(dots) are virtually identical to the DFT bands (lines). Figure 3: Wannier
bands character of (a) $e_{g}$(1), (b) $e_{g}$(2), (c) $a_{1g}$, (d) O $p$(1),
and (e) O $p$(2). (f) Schematic hopping terms of Cu site, the red and green
arrows indicate first and second nearest hopping, the numbers in bracket
indicate the real space hopping vectors. Figure 4: (Two-band) tight-binding
model fermi surface of Pb9Cu(PO4)6O.

##  II Computational methods

DFT-level structural relaxations and static calculations are performed by
employing Vasp Kresse and Hafner (1993); Kresse and Furthmüller (1996) and
Wien2K Blaha  _et al._ (2001); Schwarz and Blaha (2003) code with the Perdew-
Burke-Ernzerhof version for solids of the generalized gradient approximation
(GGA-PBESol) Perdew  _et al._ (2008) and a dense 9$\times$9$\times$13 $k$-mesh
for the unit cell of Pb9Cu(PO4)6O. The relaxed ground state crystal structure
is shown in Fig. 2(a). The two low-energy effective Hamiltonians are generated
by projecting the derived DFT bands, now computed by WIEN2K, around the Fermi
level onto Wannier functions Wannier (1937); Marzari  _et al._ (2012) using
WIEN2WANNIER Mostofi  _et al._ (2008); Kuneš  _et al._ (2010). The real-space
Wannier Hamiltonian is then transformed to momentum space using a
$\mathbf{k}$-mesh with 18,125 reducible points.

For the DMFT calculations, this DFT-derived one-particle Hamiltonian is
supplemented by a local Kanamori interaction on the Cu sites, see Section III,
and we employ the fully localized limit as double counting correction scheme
Anisimov  _et al._ (1991). O-$p$ orbitals are considered as non-interacting.
We solve the resulting many-body Hamiltonian at room temperature (298 K,
$\beta=1/(k_{B}T)=39\,\mathrm{eV}^{-1}$) within DMFT employing a continuous-
time quantum Monte Carlo (QMC) solver in the hybridization expansions Gull
_et al._ (2011) using W2dynamics Parragh  _et al._ (2012); Wallerberger  _et
al._ (2019). Real-frequency spectra are obtained with the ana_cont code
Kaufmann and Held (2021) via analytic continuation using the maximum entropy
method (MaxEnt) Gubernatis  _et al._ (1991); Sandvik (1998).

Further, the rotationally invariant DFT+$U$ scheme Liechtenstein  _et al._
(1995) with $U=3$ eV and $J=0.7$ eV on top of the regular PBE Perdew  _et al._
(1996, 1997) functional as implemented in the Vasp package was employed for
DOS and bandstructure calculations for two representative relaxed crystal
structures. The plane wave cutoff is 600 eV, 4$\times$4$\times$5 and
5$\times$4$\times$4 $k$-mesh was used for _P3_ (143) and _Pm_ (6) structures.
AFLOW-SYM Hicks  _et al._ (2018) was used for the symmetry analysis.

##  III Tight binding models

The two-band and five-band tight binding model consist of $m=1..2$ and
$m=1..5$ orbitals in the unit cell, respectively. For the two-band model these
are the $d_{xz}$ and $d_{xz}$ orbitals of the Cu-site, for the five-band model
there are two additional O $p_{x}$ and $p_{y}$ orbitals from oxygen sites and
one more $d_{z^{2}}$ orbital from the Cu-site. This motif is periodically
extended on a hexagonal $k$-grid. This tight-binding Hamiltonian ${\cal
H}^{0}$ is supplemented by a local Coulomb interaction term ${\cal H}_{int}$
on the Cu sites.

$\displaystyle{\cal H}={\cal H}^{0}+{\cal H}_{int}.$ (1)

For the non-interacting part

$\displaystyle{\cal H}^{0}=\sum_{\mathbf{k},\sigma,m,n}{\cal
H}^{0}_{m,n}(\mathbf{k}),$ (2)

we set up a tight-binding parametrization

$\displaystyle{\cal
H}^{0}_{m,n}(\mathbf{k})=\sum_{i,j}t_{im,jn}e^{i\mathbf{k}(\mathbf{R}_{i}-\mathbf{R}_{j})}c^{{\dagger}}_{i,m\sigma}c_{j,n\sigma}^{\phantom{{\dagger}}},$
(3)

where $c^{{\dagger}}_{i,m\sigma}(c_{i,m\sigma}^{\phantom{{\dagger}}})$ is the
creation (annihilation) operator, and $i$, $j$ indicate unit cells
$\mathbf{R}_{i}$, $\mathbf{R}_{j}$, while $m$, $n$ are orbital indices, and
$\sigma$ the spin index. For the interaction part,

$\displaystyle{\cal H}_{int}=\sum_{i,\sigma,\sigma^{\prime}}{\cal
H}_{int}^{\sigma,\sigma^{\prime}}(i)$ (4)

we use the Kanamori form

$\displaystyle{\cal H}_{int}^{\sigma,\sigma^{\prime}}(i)$ $\displaystyle=$
$\displaystyle U\sum_{m}n_{m,\sigma}n_{m,\sigma^{\prime}}$ $\displaystyle+$
$\displaystyle\sum_{m>n}[U^{\prime}n_{m,\sigma}n_{n,\sigma^{\prime}}+(U^{\prime}-J)n_{m,\sigma}n_{n,\sigma^{\prime}}]$
$\displaystyle-$ $\displaystyle\sum_{m,n}^{m\neq
n}J(d^{{\dagger}}_{m,\sigma}d^{{\dagger}}_{n,\sigma^{\prime}}d_{n,\sigma}d_{m,\sigma^{\prime}}+d^{{\dagger}}_{m,\sigma}d^{{\dagger}}_{m,\sigma^{\prime}}d_{n,\sigma}d_{n,\sigma^{\prime}}),$

where all number operators, $n_{m,\sigma}$, act on the same Cu site (i.e., are
in the same unit cell $i$), and the $m,n$ orbitals are restricted to the two
and three Cu orbitals for the respective two- and five-band model, defined
below.

###  III.1  Two-band tight-binding model

As only two bands cross the Fermi energy, our initial objective is to
establish a two-band model for Pb9Cu(PO4)6O. As depicted in Fig. 2(a,b), the
bands intersecting at the Fermi level (Ef) primarily arise from the Cu
$d_{yz}$ and $d_{xz}$ orbitals (in the coordinate system with $z$ aligned with
the $c$ unit cell vector). Thus, a minimal low-energy model with only these
two orbitals appears possible. It can facilitate subsequent calculations that
extend beyond the scope of DFT, as fewer orbitals require less computational
resources for complicated many-body calculations. The energy range chosen for
this two-band Wannier projections is -0.1 eV to 0.1 eV. Fig. 2(c) shows the
excellent fit of the bands in the Wannier gauge to the DFT. Truncating the
hopping amplitudes at the second nearest neighbors, yields the hopping
elements collected in Table 1.

It is worth noting that, due to the presence of a trigonal distortion, the
orthogonality between $d_{yz}$ and $d_{xz}$ bands, which is preserved in an
undistorted CuO6 octahedron, is lifted. The distortion itself can be best seen
in Fig. 2(b). It leads to the emergence of non-zero hopping terms such as the
nearest neighbor ($t_{1}$) hopping from $e_{g}$(1) to $e_{g}$(2) along the
(110) direction. For all second nearest neighbor hoppings, the predicted
values are close to zero (less than 1 meV), indicating that hopping between Cu
ions in Pb9Cu(PO4)6O can be safely restricted to first nearest neighbors.
These nearest Cu neighbors are separated in space by approximately 10 Å.

The hopping energy along the $z$-direction, $t_{z}$, is -10 meV. Even this
hopping is notably smaller than the corresponding $t_{z}$ in infinite-layer
nickelates and cuprates ($\sim$-36 meV; where the in-plane hopping is from
-370 to -450 meV) Kitatani  _et al._ (2020). That is, in contrast to cuprates
and nickelates, in-plane and out-of-plane hoppings are quite comparable; also
the Cu-Cu distance in the $z$-direction is with $\sim$7.4 Å similar to that
in-plane ones. Altogether, we can hence conclude that LK-99 has a three
dimensional electronic structure.

Furthermore, we illustrate the Fermi surface using the two-band model in Fig.
4. (The Fermi surface of the five-band model is also exactly the same as that
of the original DFT). Intriguingly, the Fermi surface of LK-99 exhibits
striking resemblances to that of UPt3 De Visser  _et al._ (1987); Sauls
(1994); Li  _et al._ (2022). In UPt3, the prevailing consensus attributes the
emergent superconductivity to heavy fermions, rather than electron-phonon
coupling. This parallel suggests that the presence of super flat bands and
correlations with $U/W$ of the order of ${\cal O}(10)$ might play a pivotal
role in driving a transition from a normal to a superconducting state.

###  III.2 Five-band tight-binding model

As the O orbitals may play an important role if LK-99 is doped with electrons
or holes, see Fig. 1, we further construct a five-band model. Here, we also
include the Cu $d_{z^{2}}$ orbital to be on the safe side. As long as this
lowest band remains firmly below the Fermi surface (i.e., is fully occupied up
to hybridizations/orbital admixing), it need not be considered in subsequent
many-body calculations. But it can be included with a simple Hartree shift,
given by the occupations (and spin polarizations) of the two Cu $e_{g}$
orbitals. If on the other hand, this Cu $d_{z^{2}}$ orbital accumulates holes,
this is indicative that further Cu orbitals besides the two $e_{g}$ orbitals
need to be included in the calculation.

For this extended five-band model, the hopping parameters are detailed in
Table 2. In congruence with our initial analysis in Fig.2(b) and the energy
ordering displayed in Fig.2(d), the on-site energies ($t_{000}$) of Cu
$e_{g}$, O $p$, and Cu $a_{1g}$ ($d_{z^{2}}$) orbitals amount to -77 meV, -366
meV, and -715 meV, respectively. Notably, focusing on the Cu $a_{1g}$ bands,
we observe isotropic hoppings of the order of $\sim$6 meV along all in-plane
directions, while a more substantial hopping is evident along the
$z$-direction. This consistently aligns with the inherent symmetry of the
$d_{z^{2}}$ orbital [Fig.2(b)]. Furthermore, as depicted in Fig.3(c), the Cu
$a_{1g}$ ($d_{z^{2}}$) band has a remarkably flat dispersion and exhibits only
minimal hybridization with all other bands, a posteriori justifying their
exclusion from a tight binding model of LK-99.

Turning to the O $p$ orbitals in Table 2, the $p(1)$ to $p(2)$ intra-orbital
hoppings can reach up to 11 meV, while inter-orbital hoppings extend to -22
meV. A particularly noteworthy observation is the strong intra-orbital hopping
between O $p$ orbitals along the $z$-direction ($t_{z}$), amounting to -85
meV. This concurs with the pronounced dispersion of O-$p$ orbitals along the
$\Gamma$-$A$ path, see Fig. 2(d). Considering the distance of approximately
$\sim$7.4 Å, such a substantial hopping is unexpected and raises the
possibility of inducing anisotropic, quasi one-dimensional transport in LK-99
(if it is a charge transfer insulator, and if LK-99 is hole doped).

We can further utilize the tight-binding model for visualizing the band
hybridization, by plotting the individual contribution of all five orbitals to
the electronic structure. This is shown in Fig. 3. While the two Cu $e_{g}$
and two O $p$ Wannier orbitals strongly admix among themselves individually,
the mixing between these two types of orbitals as well as between them and the
Cu $d_{z^{2}}$ orbital is weak. A notable exception is the quite strong
intermixing (hybridization) between Cu $d_{z^{2}}$ and O $p$ orbitals around
the H k-point. Despite the very small hopping (hybridization) between these
orbitals (e.g., the Cu-O inter-orbital hopping is only $t_{1}(-100)=3$ meV in
Table 2), the energy difference between the Cu and O bands is similarly small
at momentum H. For this reason they still strongly hybridize. For example, in
perturbation theory the admixture is $t_{1}(-100)/\Delta_{H}$, where
$\Delta_{H}$ is the energy difference at H.

###  III.3 Interaction parameters

Constrained random phase approximation (cRPA) Miyake and Aryasetiawan (2008)
calculations for other Cu-$d^{9}$ based materials Jang  _et al._ (2016); Di
Sante  _et al._ (2017); Ivashko  _et al._ (2019); Nilsson  _et al._ (2019)
suggest an intra-orbital Hubbard interaction $U=\mathcal{O}(2-3)$ eV for the
two (three) Cu $d_{xz/yz}$(+$d_{z^{2}}$) orbitals in the 2-band (5-band)
calculation. Note that the (minimal) Cu-O distance in LK-99 (2.08Å) is longer
by about 10% compared to typical cuprates, suggesting weaker hybridizations
and screening. On the other hand, the bare charge-transfer energy between
Cu-$d$ and O-$p$ orbitals, $E_{d}-E_{p}$, is much smaller than in cuprates,
possibly enhancing screening effects. Further, to account for additional
spectral-weight transfers from retarded processes, one typically uses a static
interaction $U$ that is enhanced with respect to the cRPA-value. Therefore, we
advocate $U=2.5-3.5$ eV, a Hund’s exchange $J=0.7$ eV, and an inter-orbital
interaction $U^{\prime}=U-2J$. Given the flatness of the relevant copper
bands, we expect the parent compound to be insulating for any reasonable
interaction strength. The precise value will, however, be important in
determining the insulating nature (Mott or charge transfer).

##  IV DMFT

In order to study effects of electronic correlations on top of the tight-
binding description, we perform a DMFT calculation. DMFT describes the local
dynamics of electrons, by monitoring the charge (and spin) fluctuations on a
given lattice site in the presence of the local Hubbard interactions. The
latter, in particular, penalizes occupying a site with more than one electron
per orbital and leads to a renormalization of the quasi-particle band-
structure. Or, if the interaction is strong enough, to a Mott-Hubbard
splitting of the DFT bands.

Fig. 5 shows the analytically-continued (to real frequencies) DFT+DMFT
spectral function for both (i) the two-band and (ii) the five-band model. The
(nominal) filling for the models is (i) $n=3$ and (ii) $n=9$ electrons per Cu
site. For both models we consider two different inter-orbital interactions
$U^{\prime}=2$ eV and $U^{\prime}=3$ eV, respectively, and a fixed Hund’s
exchange $J=0.7$ eV; for the intra-orbital (Hubbard) interaction we use
$U=U^{\prime}+2J$. The differences (uncertainties) between $e_{g}(1)$ and
$e_{g}(2)$ spectra emerge since we did not enforce orbital and spin symmetry.

Figure 5: Total and orbitally-resolved DMFT spectral function $A(\omega)$ for
(a) the two-band tight-binding model at $U^{\prime}=2$ eV and (b)
$U^{\prime}=3$eV, as well as (c) the five-band tight-binding model at
$U^{\prime}=2$ eV and (d) $U^{\prime}=3$eV, all at $T=298$K.

For both models, we see a clear gap at the Fermi energy. Both Cu $e_{g}$
(i.e., $d_{xz}$ and $d_{yz}$) orbitals split into a lower and an upper Hubbard
band. Note that despite this Mott-Hubbard splitting, the orbital degeneracy
remains. Because the system is not particle-hole symmetric, the weights of the
upper and lower Hubbard bands are not symmetric.

Since the lower Hubbard band describes transitions from 3 electrons on the Cu
site to 2 electrons and since 2 electrons have a singlet-triplet splitting,
the lower Hubbard band has to show a multiplet splitting. Such a splitting is
also seen in Fig. 5. In contrast, the transition from 3 electrons to 4
electrons, i.e., the upper Hubbard band, must not show such a splitting.

The gaps for the two-band model are larger than for the five-band models. The
reason for this is that in the two-band model three more orbitals contribute
to the screening, which reduces the effective interaction in a cRPA
calculation. That is, we should, on the very limited $U$ grid available and
bare an actual cRPA calculation, rather compare $U^{\prime}=2$ eV for the two-
band model with $U^{\prime}=3$ eV for the five band model.

As for the question of Mott or charge transfer insulator we need to look at
the five-band model in Fig. 5 (c,d): At $U^{\prime}=3$ eV in Fig. 5 we clearly
have a charge transfer insulator, with the first band below the Fermi energy
stemming from the oxygen $p$ orbitals. At $U^{\prime}=2$ eV we have a very
strong admixture between Cu and O bands in the first spectral contribution
below the Fermi energy. The system is neither a clear charge transfer nor a
clear Mott insulator: We are right at the crossover between a charge-transfer
to a Mott insulator. For even smaller $U$’s, given that this tendency
continues, the system might be a Mott insulator. As we do not have interaction
parameters from cRPA calculations yet, there is a substantial uncertainty in
$U$; and we cannot give a definite answer whether LK-99 is a Mott or a charge-
transfer insulator. As cRPA estimations of $U$ also have a substantial error,
it might also be beyond present-day theoretical tools to decide this question.

Let us also note the small peak at $-0.2$ eV for $U^{\prime}=3$ eV in Fig. 5
(c). This might be akin to a Zhang-Rice singlet, see e.g. Hansmann  _et al._
(2014) for a similar peak in DFT+DMFT calculations for cuprates where such a
Zhang-Rice singlet peak is visible. A difference is that, here, for LK-99,
this peak is at the upper edge of the oxygen bands, not above these. Also, we
should note that this small peak might be a maximum entropy artifact.

Figure 6: Imaginary part of the local DMFT Greens function $G$ vs. Matsubara
frequency $\nu_{n}$. Parameters as in Fig. 5.

To avoid the maximum entropy uncertainties, we further present in Fig. 6 the
local ($\mathbf{k}$-integrated) DMFT Green’s function, specifically its
imaginary part and now for (imaginary) Matsubara frequencies $\nu_{n}$. Here
$G_{\nu_{n}}\rightarrow 0$ for $\nu_{n}\rightarrow 0$ signals that there are
no states at the Fermi energy, i.e., we have an insulator. Both, for $n=9$
electrons in the five-band and $n=3$ for the three-band model,
$G_{\nu_{n}}\rightarrow 0$. LK-99 is an insulator.

##  V DFT+$U$

Previous DFT+$U$ calculations Lai  _et al._ (2023); Kurleto  _et al._ (2023)
did not find an insulating state for the lattice symmetry and the unit cell
also employed here. The degeneracy of the $d_{yz}$ and $d_{xz}$ orbital
prevents a splitting since in DFT+$U$ spin or orbital polarization is needed
for splitting off Hubbard bands. Our DFT+$U$ result of Pb9Cu(PO4)6O in Fig.
7(a) show the same.

Figure 7: DFT+$U$ band structure (left) and DOS (right) for two structures of
Pb9Cu(PO4)6O. Spin-up and spin-down contributions to bands (fatbands) and to
DOS are highlighted in blue and orange. Panel (a) shows the calculation for
the structure presented in Fig. 2(a), while panel (b) is for the structure
with a different arrangement of Pb and Cu atoms (crystal structures are given
in the insets with Cu and Pb atoms in blue and gray, respectively).

However, introducing structural distortion via various positions of Cu-doping,
like at Pb(2) in Fig. 2(a), further reduced system symmetry from _P3_ (143) of
Fig. 2(a) to _Pm_ (6) after DFT structural relaxations; cf. insets of Fig. 7.
In this lower _Pm_ symmetry, each Cu atom acquired four surrounding O ligands,
leading to a crystal field of strongly distorted tetrahedral shape. This
breaks the orbital degeneracy of the _P3_ (143) structure. In DFT+$U$, see
Fig. 7(b), static electronic correlations then strongly enhance the orbital
splitting and additionally yield a spin polarization and a corresponding spin
splitting. This ultimately gives rise to an insulating state with a similar
gap as in DMFT, but here with a full spin and orbital polarization. Note that
in DFT+$U$ $E_{f}$ happens to be just below the upper Hubbard band in Fig.
7(b).

##  VI Conclusion

To sum up, on the basis of density-functional theory, we performed a Wannier
function projection for analyzing the electronic structure of Pb9Cu(PO4)6O
(LK-99). This allows for subsequent many-body calculations. Specifically, we
constructed two distinct low-energy models: (i) a minimal two-band model
including only the Cu $d_{xz}$ and $d_{yz}$ orbitals and (ii) a five-band
model that additionally encompasses the O $p_{x}$ and $p_{y}$ orbitals as well
as the Cu $d_{z^{2}}$ orbital. These tight binding models shed some light on
LK-99: It has a three-dimensional electronic structure of the two low-energy
Cu $d_{xz}$ and $d_{yz}$ bands with next nearest-neighbor hoppings up to
distances of 10 Å. This is contrasted with a more one-dimensional
bandstructure of the O $p_{x}$ and $p_{y}$ orbitals; the latter orbitals can
become relevant if LK-99 is a charge transfer insulator. For cuprate
superconductors, higher Tc’s have been associated with smaller charge-transfer
energies Weber  _et al._ (2012), which, there, are still of the order
$\mathcal{O}(1-3)$ eV. This empirical trend for cuprates nicely aligns with an
allegedly much higher Tc of LK-99, here, since the bare Cu-O charge-transfer
energy in our five-orbital model is an order of magnitude smaller.

The ratio of interaction-to-bandwidth $U/W$ is an order of magnitude larger
than what is needed to turn a metal into a Mott insulator Gebhard (1997).
Hence, we see in DMFT a splitting of the Cu $d_{xz}$/$d_{yz}$ into a lower and
an upper Hubbard band so that LK-99 becomes either a Mott insulator or charge
transfer insulator. Given the uncertainty in $U$ we cannot say, at the moment,
which version of these closely related insulators is realized. This splitting
is a dynamical and purely electronic mechanism. In DMFT we do not observe
indications for orbital or spin symmetry breaking. In contrast to DMFT, to
find an insulator in DFT+$U$, a symmetry breaking of the Cu $e_{g}$
($d_{xz}$/$d_{yz}$) orbitals is required. For the crystal structure of Fig.
2(a) these are, however, degenerate. In our DFT+$U$ we therefore also
considered a structure with Cu on the Pb(1) site in Fig. 2(a). In this case,
the $e_{g}$ symmetry is broken and hence the Cu bands can split into Hubbard
bands, realizing an insulator also in DFT+$U$.

Such a dynamic splittings into Hubbard bands and dynamical orbital
reoccupations can match with likewise dynamic but slower Jahn-Teller phonons
(lattice distortions). For example, in manganites electronic and lattice modes
mutually support each other and thus localize charge carriers without symmetry
breaking Yang  _et al._ (2017). While phonons have not yet been calculated, a
similar or different interplay between electron and phonon dynamics might play
an important role for LK-99.

#### Acknowledgments

We thank Robert Svagera and Neven Barišić for discussions and acknowledge
funding through the Austrian Science Fund (FWF) projects I 5398, P 36213, SFB
Q-M&S (FWF project ID F86), the joint project I 6142 of FWF and the French
National Research Agency (ANR), and Research Unit QUAST by the Deutsche
Foschungsgemeinschaft (DFG; project ID FOR5249) and FWF (project ID I 5868).
L. S. is thankful for the starting funds from Northwest University.
Calculations have been mainly done on the Vienna Scientific Cluster (VSC). For
the purpose of open access, the authors have applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising from this
submission.

## References

  * Lee  _et al._ (2023a) S. Lee, J.-H. Kim, and Y.-W. Kwon, arXiv:2307.12008 (2023a).
  * Lee  _et al._ (2023b) S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, arXiv:2307.12037 (2023b).
  * Lee  _et al._ (2023c) S. Lee, J. Kim, S. Im, S. An, Y.-W. Kwon, and K. H. Auh, J. Korean Cryst. Growth Cryst. Technol. 33, 61 (2023c).
  * Wu  _et al._ (2023) H. Wu, L. Yang, B. Xiao, and H. Chang, arXiv:2308.01516 (2023).
  * Hou  _et al._ (2023) Q. Hou, W. Wei, X. Zhou, Y. Sun, and Z. Shi, arXiv:2308.01192 (2023).
  * Liu  _et al._ (2023) L. Liu, Z. Meng, X. Wang, H. Chen, Z. Duan, X. Zhou, H. Yan, P. Qin, and Z. Liu, arXiv:2307.16802 (2023).
  * Kumar  _et al._ (2023) K. Kumar, N. K. Karn, and V. P. S. Awana, arXiv:2307.16402 (2023).
  * Kumar  _et al._ (2023) K. Kumar, N. K. Karn, Y. Kumar, and V. P. S. Awana, arXiv:2308.03544 (2023).
  * Abramian  _et al._ (2023) P. Abramian, A. Kuzanyan, V. Nikoghosyan, S. Teknowijoyo, and A. Gulian, arXiv:2308.01723 (2023).
  * Guo  _et al._ (2023) K. Guo, Y. Li, and S. Jia, arXiv:2308.03110 (2023).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Lai  _et al._ (2023) J. Lai, J. Li, P. Liu, Y. Sun, and X.-Q. Chen, arXiv:2307.16040 (2023).
  * Griffin (2023) S. M. Griffin, arXiv:2307.16892 (2023).
  * Si and Held (2023) L. Si and K. Held, arXiv:2308.00676 (2023).
  * Kurleto  _et al._ (2023) R. Kurleto, S. Lany, D. Pashov, S. Acharya, M. van Schilfgaarde, and D. S. Dessau, arXiv:2308.00698 (2023).
  * Cabezas-Escares  _et al._ (2023) J. Cabezas-Escares, N. F. Barrera, C. Cardenas, and F. Munoz, arXiv:2308.01135 (2023).
  * Tao  _et al._ (2023) K. Tao, R. Chen, L. Yang, J. Gao, D. Xue, and C. Jia, arXiv:2308.03218 (2023).
  * Sun  _et al._ (2023) Y. Sun, K.-M. Ho, and V. Antropov, arXiv:2308.03454 (2023).
  * Kuroki  _et al._ (2005) K. Kuroki, T. Higashida, and R. Arita, Phys. Rev. B 72, 212509 (2005).
  * Iglovikov  _et al._ (2014) V. I. Iglovikov, F. Hébert, B. Grémaud, G. G. Batrouni, and R. T. Scalettar, Phys. Rev. B 90, 094506 (2014).
  * Aoki (2020) H. Aoki, J. Supercond. Nov. Magn. 33, 2341 (2020).
  * Gebhard (1997) F. Gebhard,  _The Mott Metal-insulator transition_ (Springer-Verlag (Berlin), 1997).
  * Zaanen  _et al._ (1985) J. Zaanen, G. A. Sawatzky, and J. W. Allen, Phys. Rev. Lett. 55, 418 (1985).
  * Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Scalapino (2012) D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Baskaran (2023) G. Baskaran, arXiv:2308.01307 (2023).
  * Bardeen  _et al._ (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 106, 162 (1957).
  * Tavakol and Scaffidi (2023) O. Tavakol and T. Scaffidi, arXiv:2308.01315 (2023).
  * Oh and Zhang (2023) H. Oh and Y.-H. Zhang, arXiv:2308.02469 (2023).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Computational materials science 6, 15 (1996).
  * Blaha  _et al._ (2001) P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, and J. Luitz, An augmented plane wave + local orbitals program for calculating crystal properties (2001).
  * Schwarz and Blaha (2003) K. Schwarz and P. Blaha, Computational Materials Science 28, 259 (2003), Proceedings of the Symposium on Software Development for Process and Materials Design.
  * 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, Phys. Rev. Lett. 100, 136406 (2008).
  * Wannier (1937) G. H. Wannier, Phys. Rev. 52, 191 (1937).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Mostofi  _et al._ (2008) A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer physics communications 178, 685 (2008).
  * Kuneš  _et al._ (2010) J. Kuneš, R. Arita, P. Wissgott, A. Toschi, H. Ikeda, and K. Held, Computer Physics Communications 181, 1888 (2010).
  * Anisimov  _et al._ (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  * Gull  _et al._ (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Parragh  _et al._ (2012) N. Parragh, A. Toschi, K. Held, and G. Sangiovanni, Phys. Rev. B 86, 155158 (2012).
  * Wallerberger  _et al._ (2019) M. Wallerberger, A. Hausoel, P. Gunacker, A. Kowalski, N. Parragh, F. Goth, K. Held, and G. Sangiovanni, Computer Physics Communications 235, 388 (2019).
  * Kaufmann and Held (2021) J. Kaufmann and K. Held, arXiv:2105.11211 (2021).
  * Gubernatis  _et al._ (1991) J. E. Gubernatis, M. Jarrell, R. N. Silver, and D. S. Sivia, Phys. Rev. B 44, 6011 (1991).
  * Sandvik (1998) A. W. Sandvik, Phys. Rev. B 57, 10287 (1998).
  * Liechtenstein  _et al._ (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996), publisher: American Physical Society.
  * Perdew  _et al._ (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 78, 1396 (1997), publisher: American Physical Society.
  * Hicks  _et al._ (2018) D. Hicks, C. Oses, E. Gossett, G. Gomez, R. H. Taylor, C. Toher, M. J. Mehl, O. Levy, and S. Curtarolo, Acta Crystallographica Section A 74, 184 (2018).
  * Kitatani  _et al._ (2020) M. Kitatani, L. Si, O. Janson, R. Arita, Z. Zhong, and K. Held, npj Quantum Materials 5, 59 (2020).
  * De Visser  _et al._ (1987) A. De Visser, A. Menovsky, and J. Franse, Physica B+C 147, 81 (1987).
  * Sauls (1994) J. Sauls, Advances in Physics 43, 113 (1994).
  * Li  _et al._ (2022) S. Li, Z. Tong, C. Shao, H. Bao, T. Frauenheim, and X. Liu, The Journal of Physical Chemistry Letters 13, 4289 (2022).
  * Miyake and Aryasetiawan (2008) T. Miyake and F. Aryasetiawan, Phys. Rev. B 77, 085122 (2008).
  * Jang  _et al._ (2016) S. W. Jang, H. Sakakibara, H. Kino, T. Kotani, K. Kuroki, and M. J. Han, Scientific Reports 6, 33397 (2016).
  * Di Sante  _et al._ (2017) D. Di Sante, A. Hausoel, P. Barone, J. M. Tomczak, G. Sangiovanni, and R. Thomale, Phys. Rev. B 96, 121106 (2017).
  * Ivashko  _et al._ (2019) O. Ivashko, M. Horio, W. Wan, N. B. Christensen, D. E. McNally, E. Paris, Y. Tseng, N. E. Shaik, H. M. Rønnow, H. I. Wei, C. Adamo, C. Lichtensteiger, M. Gibert, M. R. Beasley, K. M. Shen, J. M. Tomczak, T. Schmitt, and J. Chang, Nature Communications 10, 786 (2019).
  * Nilsson  _et al._ (2019) F. Nilsson, K. Karlsson, and F. Aryasetiawan, Phys. Rev. B 99, 075135 (2019).
  * Hansmann  _et al._ (2014) P. Hansmann, N. Parragh, A. Toschi, G. Sangiovanni, and K. Held, New J. Phys. 16, 033009 (2014).
  * Weber  _et al._ (2012) C. Weber, C. Yee, K. Haule, and G. Kotliar, EPL (Europhysics Letters) 100, 37001 (2012).
  * Yang  _et al._ (2017) S.-L. Yang, J. A. Sobota, Y. He, Y. Wang, D. Leuenberger, H. Soifer, M. Hashimoto, D. H. Lu, H. Eisaki, B. Moritz, T. P. Devereaux, P. S. Kirchmann, and Z.-X. Shen, Phys. Rev. B 96, 245112 (2017).