# Revisiting holographic model for thermal and dense QCD with a critical point

Qingxuan Fu<EMAIL_ADDRESS>Song He<EMAIL_ADDRESS>Li
Li 3,4,5333<EMAIL_ADDRESS>and Zhibin Li<EMAIL_ADDRESS>
###### Abstract

To quantitatively provide reliable predictions for the hot and dense QCD
matter, a holographic model should be adjusted to describe first-principles
lattice results available at vanishing baryon chemical potential. The equation
of state from two well-known lattice groups, the HotQCD collaboration and the
Wuppertal-Budapest (WB) collaboration, shows visible differences at high
temperatures. We revisit the Einstein-Maxwell-dilaton (EMD) holographic model
for hot QCD with 2+1 flavors and physical quark masses by fitting lattice QCD
data from the WB collaboration. Using the parameterization for the scalar
potential and gauge coupling proposed in our work [Phys.Rev.D 106 (2022) 12,
L121902], the equation of state, the higher order baryon number
susceptibilities, and the chiral condensates are in quantitative agreement
with state-of-the-art lattice results. We find that the critical endpoint
(CEP) obtained from fitting the WB collaboration data is nearly identical to
the one from the HotQCD collaboration, suggesting the robustness of the
location of the CEP. Moreover, our holographic prediction for the CEP location
is in accord with more recent Bayesian analysis on a large number of
holographic EMD models and an effective potential approach of QCD from gap
equations.

1 Center for Theoretical Physics and College of Physics, Jilin University,
Changchun 130012, China
2 Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Muhlenberg 1, 14476 Golm, Germany
3 CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China
4 School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute
for Advanced Study, University of Chinese Academy of Sciences, Hangzhou
310024, China
5 Peng Huanwu Collaborative Center for Research and Education, Beihang
University, Beijing 100191, China
6 Institute for Astrophysics, School of Physics, Zhengzhou University,
Zhengzhou 450001, China

###### Contents

  1. 1 Introduction
  2. 2 Holographic setup and parameter fixing
  3. 3 Thermodynamics and CEP at finite $\mu_{B}$
  4. 4 Critical phenomena near the CEP
  5. 5 Chiral condensate
  6. 6 Conclusion and discussion
  7. 7 Acknowledgments

##  1 Introduction

Studying the properties of the quark-gluon matter is essential to
understanding the theory of the strong interactions known as Quantum
Chromodynamics (QCD), which is the fundamental force describing the
interactions between quarks and gluons. It is the main goal of heavy ion
collision experiments at the Relativistic Heavy Ion Collider (RHIC) and the
Large Hadron Collider (LHC). However, obtaining a quantitative understanding
of the QCD phase structure at finite temperature $T$ and baryon chemical
potential $\mu_{B}$ remains challenging due to the strongly coupled nature of
the system under extreme conditions.

Several non-perturbative approaches have been proposed to study the QCD phase
diagram under various conditions. Lattice QCD, formulated on a grid of points
in space and time, provides a first-principle computation at zero baryon
chemical potential and gives reliable information at small $\mu_{B}$ by
extrapolating the lattice data. Lattice QCD suggests that the chiral and
confinement/deconfinement phase transitions occur as an analytic crossover for
small chemical potentials [1, 2, 3]. On the other hand, effective theories
such as the Dyson-Schwinger equation (DSE) [4, 5, 6, 7, 8, 9], the Nambu-Jona-
Lasinio (NJL) model [10, 11, 12, 13] and the functional renormalization group
(FRG) [14, 15, 16] have suggested that the crossover will become a first-order
phase transition as $\mu_{B}$ increases. The critical point at the endpoint of
the line of the first-order QCD phase transitions is known as the QCD critical
endpoint (CEP). Despite decades of theoretical and experimental efforts,
neither an exact location nor the properties of the CEP are well known.
Nevertheless, currently lattice QCD results disfavor the existence of the CEP
for $\mu_{B}/T\leq 3$ and $\mu_{B}<300~{}\text{MeV}$ [17, 18, 19, 20, 21, 22].

An alternative non-perturbative approach is the so-called holographic QCD by
the gauge gravity duality, which maps the strongly coupled non-Abelian gauge
theories into a weakly coupled gravitational system with one higher dimension.
It aims to capture the essential characteristics of realistic QCD and confront
lattice QCD data at a quantitative level. Of particular interest is the
Einstein-Maxwell-Dilaton (EMD) theory that was originally introduced in [23,
24]. The holographic EMD theory has been used to study various issues related
to the quark matter in the literature, see [25, 26, 27] for recent reviews. In
recent studies, efforts have been made to improve the holographic QCD models
and achieve a better quantitative description for hot QCD with 2+1 flavors and
physical quark masses, see _e.g._ [28, 29, 30, 31, 32]. Non-perturbative
effects and flavor dynamics are effectively incorporated into the model
parameters by matching with lattice QCD data. These improved models have shown
qualitative consistency with the expected QCD phase diagram, including the
existence of a CEP. Chiral condensation and gluon dynamics have also been
studied, and the results are in quantitative agreement with lattice QCD
simulations.

For holographic EMD models, the biggest difficulty is to find the scalar
potential $V(\phi)$ and the coupling between scalar and gauge field $Z(\phi)$.
As an effective field theory, the model parameters of the bulk gravitational
theory should be fixed by matching with lattice QCD results, for which it is
crucial to use up-to-date lattice simulation to make reliable predictions at
finite $\mu_{B}$. In our recent work [30], we provided a parameterization for
$V(\phi)$ and $Z(\phi)$ with the five parameters fixed completely through
lattice simulations at $\mu_{B}=0$: the equation of state and the second-order
baryon number susceptibility. By adopting the lattice QCD data from the HotQCD
collaboration for 2+1 quark flavors, we have obtained quantitative
descriptions for QCD matter [33, 34]. In particular, we predicted the exact
location of the CEP with $(T_{CEP},\mu_{CEP})=(105,555)~{}\text{MeV}$. A
constraint on the holographic QCD phase transition from pulsar timing array
observations can be found in [35]. Surprisingly, our parameterization is able
to also capture hot and dense QCD matter for 2 flavors [36] and pure gluon
[37].

Besides the HotQCD collaboration, another mainstream lattice QCD group is the
Wuppertal-Budapest collaboration (WB collaboration). The WB collaboration
adopts a different setting compared to HotQCD, resulting in a different
equation of state [38, 2]. In particular, the equation of state at high $T$
starts to disagree even when considering the error bars. It is therefore
natural to ask if the position of CEP is sensitive to the setting from both
lattice QCD groups. In this study, we will use the parameterization in our
previous work [30] to fit the lattice data of WB collaboration at zero density
[2]. We will compare our results for the higher order baryon number
susceptibilities with the lattice computation from the WB collaboration [38],
as well as the equation of state at finite but small $\mu_{B}/T$ [20]. Both
will give good quantitative agreement. Moreover, the position of the CEP
obtained from fitting the WB collaboration data is nearly identical to the one
we obtained by using the HotQCD data [30]. It suggests that the location of
CEP is pretty robust, not sensitive to the difference in the equation of
state. Moreover, we will also show our holographic computation of the chiral
condensation for the light quarks and the strange quark.

The subsequent sections of this work are organized as follows: In Section 2,
we provide a brief introduction to our holographic setup and compare our
holographic results to the lattice QCD data from the WB collaboration at zero
baryon chemical potential. Section 3 is dedicated to presenting the QCD phase
structure derived from our holographic model, along with detailing the
location of the CEP. In Section 4, we delve into the computation of critical
exponents in the vicinity of the CEP. We present our holographic computation
of the chiral condensates for the light and strange quarks in Section 5.
Finally, we provide a summary and discussion in Section 6.

##  2 Holographic setup and parameter fixing

The gravitational action of the five-dimensional EMD theory is given by

$\displaystyle\begin{gathered}S=\frac{1}{2\kappa_{N}^{2}}\int
d^{5}x\sqrt{-g}\left[\mathcal{R}-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi-\frac{Z(\phi)}{4}F_{\mu\nu}F^{\mu\nu}-V(\phi)\right]+S_{\partial}\,.\end{gathered}$
(2.2)

where $\kappa_{N}^{2}$ is the effective Newton constant and $S_{\partial}$ is
the boundary term. The dilaton $\phi$ is responsible for breaking the
conformal invariance of boundary field theory and $A_{\mu}$ introduces the
baryon chemical potential. Non-perturbative effects are effectively adopted
into $V(\phi)$ and $Z(\phi)$ which are parameterized as [30]

$\displaystyle
V(\phi)=-12\cosh\left[c_{1}\phi\right]+\left(6c_{1}^{2}-\frac{3}{2}\right)\phi^{2}+c_{2}\phi^{6},$
(2.3) $\displaystyle
Z(\phi)=\frac{1}{1+c_{3}}\operatorname{sech}\left[c_{4}\phi^{3}\right]+\frac{c_{3}}{1+c_{3}}e^{-c_{5}\phi}\,,$
(2.4)

with five free parameters $c_{1},c_{2},c_{3},c_{4},c_{5}$.

The charged hairy planer black bole takes

$\displaystyle\begin{gathered}ds^{2}=-f(r)e^{-\eta(r)}dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}(dx^{2}+dy^{2}+dz^{2})\,,\\\
\phi=\phi(r),\quad A=A_{t}(r)dt\,,\end{gathered}$ (2.7)

with $r$ the holographic radial coordinate. The black hole event horizon is
given by the largest root of $f(r_{h})=0$, and the asymptotically AdS boundary
is located at $r\rightarrow\infty$. The asymptotic expansion near the AdS
boundary reads

$\begin{split}\phi&=\frac{\phi_{s}}{r}+\cdots+\frac{\phi_{v}}{r^{3}}+\cdots,\quad
A_{t}=\mu_{B}-\frac{\kappa_{N}^{2}n_{B}}{r^{2}}+\cdots\,,\\\
f&=r^{2}+\cdots+\frac{f_{v}}{r^{2}},\qquad\qquad\;\eta=0+\cdots\,,\end{split}$
(2.8)

with $\phi_{s},\phi_{v},\mu_{B},n_{B},f_{v}$ constants that will be fixed by
solving the equations of motion. The source term $\phi_{s}$ essentially breaks
the conformal symmetry and plays the role of the energy scale.

In our coordinate system, the boundary term for a well-defined Dirichlet
variational principle and for removing divergence is given by

$\displaystyle\begin{aligned} S_{\partial}&=S_{\text{Gibbons-
Hawking}}+S_{\text{counter}}\\\
&=\frac{1}{2\kappa_{N}^{2}}\int_{r\rightarrow\infty}dx^{4}\sqrt{-h}\left[2K-6-\frac{1}{2}\phi^{2}-\frac{6c_{1}^{4}-1}{12}\phi^{4}\ln[r]-b\phi^{4}+\frac{1}{4}F_{\rho\lambda}F^{\rho\lambda}\ln[r]\right].\end{aligned}$
(2.9)

Here $h$ is the determinant of the induced metric at the AdS boundary and $K$
is the trace of the extrinsic curvature defined by the outward pointing normal
vector to the boundary. The parameter $b$ is determined by requiring the
pressure $P(T=0)=0$ at vanishing chemical potential. Via the standard
holographic dictionary, we can then obtain many thermodynamic observables from
the UV data (2.8). In particular, $\mu_{B}$ and $n_{B}$ are baryon chemical
potential and baryon number density, respectively. The free energy density,
energy density, and pressure are given by

$\displaystyle\begin{aligned}
\Omega=&\frac{1}{2\kappa_{N}^{2}}\left(f_{v}-\phi_{s}\phi_{v}-\frac{3-48b-8c_{1}^{4}}{48}\phi_{s}^{4}\right)\,,\\\
\epsilon=&\frac{1}{2\kappa_{N}^{2}}\left(-3f_{v}+\phi_{s}\phi_{v}+\frac{1+48b}{48}\phi_{s}^{4}\right)\,,\\\
P=&\frac{1}{2\kappa_{N}^{2}}\left(-f_{v}+\phi_{s}\phi_{v}+\frac{3-48b-8c_{1}^{4}}{48}\phi_{s}^{4}\right)\,.\end{aligned}$
(2.10)

Then one obtains the trace anomaly $I=\epsilon-3P$. Meanwhile, the temperature
and entropy density can be obtained at the event horizon $r=r_{h}$.

$\displaystyle
T=\frac{1}{4\pi}f^{\prime}\left(r_{h}\right)e^{-\eta\left(r_{h}\right)/2},\hskip
22.76219pts=\frac{2\pi}{\kappa_{N}^{2}}r_{h}^{3}\,.$ (2.11)

The squared sound speed and the $n$-th order baryon number susceptibility are
given by

$\displaystyle c_{s}^{2}=\left.\frac{dP}{d\epsilon}\right|_{\mu_{B}},\hskip
22.76219pt\chi_{n}^{B}(T,\mu_{B})=\left.\frac{\partial^{n}}{\partial(\mu_{B}/T)^{n}}\frac{P}{T^{4}}\right|_{T}\,.$
(2.12)

The baryon number susceptibilities are closely related to various cumulants of
the baryon number distribution measured in heavy-ion collision experiments. In
particular, at $\mu_{B}=0$ one has

$\displaystyle\chi_{2}^{B}(T,0)=\left.\frac{\partial^{2}}{\partial(\mu_{B}/T)^{2}}\frac{P}{T^{4}}\right|_{T}=\lim_{\mu_{B}\rightarrow
0}\frac{1}{T^{2}}\frac{n_{B}}{\mu_{B}}\,.$ (2.13)

The hairy black hole solutions will be obtained numerically after further
requiring the regularity at the event horizon. Please refer to [30, 39] for
more details about the thermodynamics of hairy black holes.

As mentioned above, we fix all free parameters of the EMD model by matching
the equation of state and the second-order baryon number susceptibility
evaluated at $\mu_{B}=0$ with the corresponding lattice QCD data of WB
collaboration [38, 2]. We obtain that

$\displaystyle\begin{aligned} &c_{1}=0.7055,\quad c_{2}=0.0038,\quad
c_{3}=1.7877,\quad c_{4}=0.1050,\quad c_{5}=30.2054\,,\\\
&\phi_{s}=1080~{}\text{MeV},\quad\kappa_{N}{}^{2}=(2\pi)1.76,\quad
b=-0.2554\,.\end{aligned}$ (2.14)

These values are close to the one by fitting the lattice data of HotQCD
collaboration [30], which is expected since the results from both lattice
groups show small but visible differences at high temperatures. The fitting
results are depicted in Fig. 1. The left panel involves the comparison of the
normalized entropy density, trace anomaly, and pressure obtained from our
holographic model and the corresponding lattice data [2]. The right panel
shows the squared speed of sound. One finds a good agreement with the lattice
data over a considerable temperature range from $T=100~{}\text{MeV}$ to
$500~{}\text{MeV}$.

Figure 1: (Left) Dimensionless entropy density $s/T^{3}$, trace anomaly
$(\epsilon-3P)/T^{4}$ and pressure $P/T^{4}$ with (Right) squared speed of
sound $c_{s}^{2}$ at zero baryon density. Solid curves represent results from
our holographic model, and data with error bars are lattice QCD data from WB
collaboration [2].

Once obtaining the equation of state, one can compute the baryon number
susceptibilities $\chi_{n}^{B}$ using (2.12). They are the coefficients in the
Taylor expansion of $P$ around $\mu_{B}=0$, thus can reflect the information
of the equation of state at small $\mu_{B}$. The baryon number
susceptibilities at $\mu_{B}=0$ serve as an important probe to test the
accuracy of our model. In Fig. 2, we present a direct comparison between
$\chi_{n}^{B}$ obtained from our holographic model and the corresponding
lattice QCD computation from WB collaboration. The top left panel of Fig. 2
shows $\chi_{2}^{B}$ and $Td\chi_{2}^{B}/dT$, which we have used to determine
the parameters $c_{3},c_{4},c_{5}$ in the coupling function $Z(\phi)$. The
other three panels depict the temperature dependence of higher-order
susceptibilities. In all cases, we observe good quantitative consistency,
further demonstrating the accuracy of our holographic model.

Figure 2: Baryon number susceptibilities $\chi^{B}_{2}$ and
$Td\chi^{B}_{2}/dT$ (top left), $\chi^{B}_{3}/\chi^{B}_{1}$ (top right),
$\chi^{B}_{4}/\chi^{B}_{2}$ (bottom left) and $\chi^{B}_{6}$ (bottom right) at
$\mu_{B}=0$ compared with lattice data [40, 38, 20]. The light-gray band
denotes the region of continuous limit from lattice QCD simulation.

##  3 Thermodynamics and CEP at finite $\mu_{B}$

Using a novel expansion scheme, the most up-to-date lattice data for the QCD
equation of state is available up to $\mu_{B}/T=3.5$ [20]. To further validate
the accuracy of our holographic model, we compare the holographic prediction
at finite $\mu_{B}$ with the one recently reported by [20]. In Fig 3, we show
the normalized entropy density, pressure, energy density, and baryon number
density as a function of temperature for various $\mu_{B}/T$. Our holographic
predictions for the entropy density, pressure, and energy density are in
quantitative agreement with the lattice results all the way up to
$\mu_{B}/T=3.5$. Regarding the baryon number density, there shows also
quantitative agreement for most of the values of $T$ and $\mu_{B}$, although
the holographic model overestimates the lattice results for $n_{B}$ at high
temperatures $T\geq 195~{}\text{MeV}$ when $\mu_{B}/T\geq 3.0$. Nevertheless,
this further supports the reliability and effectiveness of our holographic
model in describing the thermodynamic properties of quark matter at nonzero
$\mu_{B}$.

Figure 3: Dimensionless entropy density $s/T^{3}$, pressure $P/T^{4}$, baryon
density $\rho/T^{3}$, and energy density $\epsilon/T^{4}$ as functions of
temperature $T$ at various $\mu_{B}/T$. The solid curves obtained from our
holographic model are compared with the lattice QCD data [20], represented by
dots with error bars.

After solving the bulk equations of motion, we can obtain all extrema of the
free energy that corresponds to the coexistence region of not only
thermodynamically stable minima but also metastable and unstable saddle
points. The behavior of free energy density $\Omega$ versus temperature at
fixed baryon chemical potential is depicted in Fig. 4. At low $\mu_{B}$,
$\Omega$ is a single-valued function of $T$, corresponding to a crossover
without a clear phase boundary. In contrast, one observes the characteristic
multivalued swallow-shape for $\Omega$, yielding a first-order phase
transition. By identifying the physical state corresponding to the minimum
free energy density, we can locate the boundary of the first-order phase
transitions.

Figure 4: Free energy density $\Omega$ as a function of temperature $T$ at
different baryon chemical potential $\mu_{B}$. As $\mu_{B}$ increases, the
smooth crossover becomes a first-order one.

The $T$-$\mu_{B}$ phase diagram obtained from our holographic approach is
shown in Fig. 5. The black solid line represents the line of the first-order
phase transitions and the red dashed one denotes the inflection point of
$\chi_{2}^{B}$ that characterizes the crossover region. The CEP at the end of
the first-order phase transition line is indicated by a red dot. More
precisely, the location of the CEP is at
$(T_{CEP},\mu_{CEP})=(109,552)~{}\text{MeV}$. Remarkably, the CEP position
obtained by fitting the WB lattice data is very close to the one we reported
by fitting the hotQCD data, $(T_{CEP},\mu_{CEP})=(105,555)~{}\text{MeV}$ [30].
In contrast, an early study [29] fitting the WB lattice data gave a
significantly different location of the CEP,
$(T_{CEP},\mu_{CEP})=(89,724)~{}\text{MeV}$. Although a relatively complicated
form of $V$ and $Z$ was adopted in [29] (see [26] for a review), our model
improves the quality of the fitting for the same lattice data. For example,
one can compare the trace anomaly in Fig. 2 of [41] with our case in the left
panel of Fig. 1 and the baryon number density in Fig. 12 of [41] with our case
in Fig 3. We also note that the scaling dimension of the gauge theory operator
dual to the dilaton field is $\Delta\approx 2.73$ of [41], while in our
approach $\Delta=3$, which is much easier to extract the UV date from (2.8)
numerically. A more recent study [31] incorporated the potential
reconstruction method together with machine learning. The location of CEP for
2+1 flavors was found to be at $(T_{CEP},\mu_{CEP})=(94,740)~{}\text{MeV}$
which is significantly from our results. However, as visible from Fig. 4 of
[31], this holographic model overestimates the lattice results for the trace
anomaly for $T\leq 200\text{MeV}$. Meanwhile, the baryon number susceptibility
$\chi^{B}_{2}$ displayed in Fig. 5 of [31] yields a poor agreement with the
lattice data.

Figure 5: Holographic prediction for the QCD phase diagram in $T$-$\mu_{B}$
plane. The red dashed line represents the crossover boundary, determined by
the inflection point of $\chi_{2}^{B}$, and the black solid curve is the
first-order transition line obtained by minimum free energy density. Between
them there exists a critical endpoint(bold red dot) with
$T_{CEP}=109~{}\text{MeV}$, $\mu_{CEP}=552~{}\text{MeV}$. CEP predicted by
other QCD models is also presented, including D1-D6 from Dyson-Schwinger
equation(DSE) [4, 5, 6, 7, 8, 9], F1-F2 from functional renormalization
group(FRG) [14, 15], Polyakov-Nambu-Jona-Lasinio model [12], N1 from Nambu-
Jona-Lasinio [10], E1 for effective potential approach [42], as well as H1-H3
from different holographic QCD models [23, 28, 29]. The pink and gray areas
provide $68\%$ and $95\%$ confidence regions for the CEP location under
different parametric ansatz [43]. Surprisingly, the position of our CEP is
completely within the $95\%$ confidence region.

Nevertheless, it is worth noting that a recent study [43] shows that the
predictions for the CEP location in different realizations of the model from a
Bayesian analysis overlap at one sigma, with $T_{CEP}=101-108\ \text{MeV}$ and
$\mu_{CEP}=560-625\ \text{MeV}$.555Note that the location of the CEP reported
by [31] is outside the one sigma region of [43]. A better choice of function
configurations for the deformed factor and gauge coupling function in [31]
could give potential improvement on the fitting quality with lattice QCD.
Meanwhile, the CEP from an effective potential approach is found to be located
at $(T_{CEP},\mu_{CEP})=(101.3,558)~{}\text{MeV}$ [42]. Both agree with our
CEP, suggesting that our prediction for the location of the CEP is robust.

##  4 Critical phenomena near the CEP

The critical exponents are vital parameters used to describe critical points
and are regarded as universal physical quantities. They transcend the
specifics of a physical system and depend solely on the system’s degrees of
freedom and correlation length. As one approaches the critical endpoint (CEP),
thermodynamic quantities often exhibit power-law scaling with temperature or
chemical potential. Consequently, along different directions, a set of
critical exponents can be defined to provide a more nuanced characterization
of the system’s critical behavior.

The exponent $\alpha$ and $\gamma$ are defined along the first-order axis
which is the tangent of the first-order line near the CEP, with the power law
of the specific heat and baryon number susceptibility, respectively.

$C_{\rho}\equiv T\left(\frac{\partial s}{\partial
T}\right)_{\rho_{B}}=-T\left(\frac{\partial^{2}\Omega}{\partial
T^{2}}-\frac{(\partial^{2}\Omega/\partial
T\partial\mu)^{2}}{(\partial^{2}\Omega/\partial\mu^{2})}\right)\sim|T-T_{\text{CEP}}|^{-\alpha}.$
(4.1) $\chi^{B}_{2}=\frac{1}{T^{2}}\left(\frac{\partial
n_{B}}{\partial\mu_{B}}\right)_{T}\sim|T-T_{\text{CEP}}|^{-\gamma}.$ (4.2)

On either side of the first-order phase transition line, which is determined
by the minimum free energy, thermodynamic quantities display two distinct
branches corresponding to the high-temperature and low-temperature phases.
These branches typically exhibit discontinuities at the phase boundary,
enabling us to observe significant changes in various thermodynamic quantities
during the first-order phase transition. The critical exponent $\beta$ is
related to the power-law relationship of the discontinuous change in entropy
density $s$ along the first-order phase transition line and the temperature
$T$.

$\Delta s=s_{>}-s_{<}\sim(T_{\text{CEP}}-T)^{\beta}.$ (4.3)

Where $s_{>}$ and $s_{<}$ represent entropy density from the high-temperature
and low-temperature branches, respectively. And the critical exponent $\delta$
is defined follow the relation of $s-s_{\text{CEP}}$ and
$\mu_{B}-\mu_{\text{CEP}}$ at the critical isotherm $T=T_{\text{CEP}}$.

$s-s_{\text{CEP}}\sim|\mu_{B}-\mu_{\text{CEP}}|^{1/\delta}\,.$ (4.4)

Fig. 6 shows the behavior of different thermodynamic quantities near the
critical endpoint. And critical exponents are obtained by fitting the slope of
different log plots.

Figure 6: The critical exponents $\alpha$, $\beta$, $\gamma$ and $\delta$ from
our holographic model.

The values of the four critical exponents are shown in Table 1. We also
compared the values of the critical exponents near CEP given by our
holographic model and the non-QCD fluids experiment data, the results of 3D
Ising model, the mean-field theory, and the holographic QCD model (the DGR
model) [44, 23]. Surprisingly, our results are closest to the results of the
mean-field theory, which is consistent with the results of the holographic QCD
model [23].

| Experiment | 3D Ising | Mean field | DGR model | HQCD
---|---|---|---|---|---
$\alpha$ | 0.110-0.116 | 0.110(5) | 0 | 0 | 1.9739 $\times$ $10^{-7}$(1603)
$\beta$ | 0.316-0.327 | 0.325$\pm$0.0015 | 1/2 | 0.482 | 0.5098(67)
$\gamma$ | 1.23-1.25 | 1.2405$\pm$0.0015 | 1 | 0.942 | 1.0159(62)
$\delta$ | 4.6-4.9 | 4.82(4) | 3 | 3.035 | 3.0399(843)
Table 1: Critical exponents from experiments in non-QCD fluids, the full
quantum 3D Ising model, mean-field (van der Waals) theory, the DGR model, and
our 2+1-flavor holographic model.

##  5 Chiral condensate

To incorporate flavor dynamics into our model, we adopt techniques akin to
those in the KKSS model [45]. This involves introducing a holographic probe
action to describe the chiral condensates. With the hairy bulk geometry (2.7)
for the (2+1)-flavor system, we adopt the following effective action for the
light quakes ($u$ and $d$ quarks) and the strange quark ($s$ quark).

$\displaystyle\begin{gathered}S_{X_{q}}=\frac{1}{2\kappa_{N}^{2}}\int
d^{5}x\sqrt{-g}Z_{q}(\phi)\left[-\frac{1}{2}\nabla_{\mu}X_{q}\nabla^{\mu}X_{q}-V\left(X_{q}\right)\right]+S_{X_{q},\partial}\,,\end{gathered}$
(5.2)

with $q=l$ for the light quarks and $q=s$ for the $s$ quark. The boundary term
$S_{X_{q},\partial}$ is essential to obtain the physical condensation. The
bulk scalar field $X_{q}$ is dual to the chiral operator
$\bar{\psi}_{q}\psi_{q}$ with the scaling dimension $\Delta_{c}=3$. Following
our previous work [30], the coupling functions $Z_{q}(\phi)$ and $V(X_{q})$
take

$\displaystyle\begin{gathered}Z_{q}(\phi)=a_{1}e^{a_{q}\phi^{2}},\quad
V\left(X_{q}\right)=-\frac{3}{2}X_{q}^{2}+a_{0}X_{q}^{4},\quad
q=l,s\,,\end{gathered}$ (5.4)

where $a_{0},a_{1},a_{q}$ are parameters that will be fixed by fitting lattice
QCD data.

The equation of motion on top of (2.7) reads,

$\displaystyle
X_{q}^{\prime\prime}+\left(\frac{f^{\prime}}{f}-\frac{\eta^{\prime}}{2}+\frac{3}{r}\right)X_{q}^{\prime}+\frac{\partial_{\phi}Z_{q}X_{q}{}^{\prime}\phi^{\prime}}{Z_{q}}-\frac{1}{f}\partial_{X_{q}}V=0\,,$
(5.5)

and the behavior of $X_{q}$ near the AdS boundary is given by

$\displaystyle
X_{q}(r)=\frac{m_{q}}{r}+\cdots+\frac{\sigma_{q}}{r^{3}}+\cdots.$ (5.6)

Here $\sigma_{q}$ is a constant and $m_{q}$ is nothing but the mass of each
quark. Matching lattice QCD data, we choose $m_{l}=5.1\mathrm{MeV}$ and
$m_{s}=102\mathrm{MeV}$. The chiral action suffers from divergence and should
be renormalized with a counter term that takes

$\displaystyle\begin{aligned}
S_{X_{q},\partial}=\frac{1}{2\kappa_{N}^{2}}\int_{r\rightarrow\infty}dx^{4}\sqrt{-h}\left[-\frac{1}{2}a_{1}X_{q}^{2}+a_{1}a_{0}X_{q}^{4}\ln[r]+\frac{a_{1}\left(1-6a_{q}\right)}{6}X_{q}^{2}\phi^{2}\ln[r]\right]\,.\end{aligned}$
(5.7)

Then we can obtain the chiral condensate

$\displaystyle\begin{aligned}
\langle\bar{\psi}\psi\rangle_{q,T}=\frac{\delta\left(S_{X_{q}}+S_{X_{q},\partial}\right)}{\delta
m_{q}}=\frac{a_{1}}{2\kappa_{N}^{2}}\left[2\sigma_{q}+2a_{0}m_{q}^{3}+\frac{1}{4}m_{q}\phi_{s}^{2}\right].\end{aligned}$
(5.8)

To compare directly with the W-B lattice simulation, we define the
renormalized chiral condensate and subtracted chiral condensate as [46]

$\displaystyle\begin{aligned}
\langle\bar{\psi}\psi\rangle_{R}=-\frac{m_{l}}{X^{4}}\left[\langle\bar{\psi}\psi\rangle_{l,T}-\langle\bar{\psi}\psi\rangle_{l,0}\right]~{}\text{and}~{}~{}\Delta_{l,s}=\frac{\langle\bar{\psi}\psi\rangle_{l,T}-\frac{m_{l}}{m_{s}}\langle\bar{\psi}\psi\rangle_{s,T}}{\langle\bar{\psi}\psi\rangle_{l,0}-\frac{m_{l}}{m_{s}}\langle\bar{\psi}\psi\rangle_{s,0}}.\end{aligned}$
(5.9)

with $\langle\bar{\psi}\psi\rangle_{q,0}$ the chairl condensate at
$T=0$.666$\langle\bar{\psi}\psi\rangle_{q,0}$ can be obtained either by
treated as a free parameter to be adjusted to fit lattice data or by
evaluating $\langle\bar{\psi}\psi\rangle_{q,T}$ at sufficiently low
temperatures. Both yield identical results,
$\langle\bar{\psi}\psi\rangle_{l,0}=1.81\times 10^{7}$,
$\langle\bar{\psi}\psi\rangle_{s,0}=-1.033\times 10^{8}$. The temperature
dependence of the renormalized chiral condensate and subtracted chiral
condensate are depicted in Fig. 7 with
$a_{0}=13,a_{1}=2.76,a_{l}=0.58,a_{s}=0.8$ and $X=m_{\pi}=135$ MeV. One finds
a remarkably quantitative agreement with lattice data [46]. It could serve as
a strong support for our holographic setup.

Figure 7: Renormalized chiral condensate and subtracted chiral condensate
compare with lattice data [46].

Our holographic predictions for the renormalized chiral condensate
$\langle\bar{\psi}\psi\rangle_{R}$ and subtracted chiral condensate
$\Delta_{l,s}$ as a function of temperature are displayed in Fig. 8 along
lines of constant $\mu_{B}$. One finds that both the renormalized chiral
condensate $\langle\bar{\psi}\psi\rangle_{R}$ and subtracted chiral condensate
$\Delta_{l,s}$ for different $\mu_{B}$ tend to a fixed value at high
temperatures. $\langle\bar{\psi}\psi\rangle_{R}$ is promoted while
$\Delta_{l,s}$ is suppressed as the baryon chemical potential is increased.
For a given $\mu_{B}$, $\Delta_{l,s}$ is larger at low temperatures and
decreases as the temperature increases. At high temperatures, the renormalized
chiral condensates $\langle\bar{\psi}\psi\rangle_{R}$ all decrease to negative
values. One also notices the formation of a characteristic multivalued S-shape
for both chiral condensates when $\mu_{B}$ is larger than the critical value
$\mu_{CEP}=552~{}\text{MeV}$, suggesting the development of a first-order
phase transition.

Figure 8: Temperature dependence of the renormalized chiral condensate and
subtracted chiral condensate for different values of $\mu_{B}$. The
characteristic multivalued S-shape develops when $\mu_{B}$ is beyond the
critical value $\mu_{CEP}=552~{}\text{MeV}$.

##  6 Conclusion and discussion

The motivation for our study arises from a unique situation in lattice QCD:
there exist only two sets of lattice configurations that provide the equation
of state (EoS) for QCD matter at zero chemical potential and finite
temperature. While these lattice QCD results are quantitatively consistent in
the low-temperature regime, they exhibit visible discrepancies at high
temperatures. The discrepancies stem from differences in the choice of fermion
propagators and initial configurations employed by the two lattice groups,
resulting in markedly distinct EoS for hot QCD matter.

Currently, there is no definitive criterion to discern the discrepancies
between these two sets of lattice data. To address this issue, we have
previously developed a holographic QCD model that successfully reproduces the
EoS consistent with one of these lattice groups, known as ”hot-QCD.” Given the
inconsistent behavior of lattice data in the high-temperature region, it
becomes imperative to establish another holographic QCD model that
quantitatively predicts the critical point in the QCD phase diagram based on
the data from the other lattice group, referred to as ”WB”.

The precise location of this critical point is of paramount importance, as it
exhibits high sensitivity to the model parameters. Thus, the primary
motivation behind this work is to thoroughly investigate and determine the
exact position of the critical point, addressing the discrepancies in the
lattice data and advancing our understanding of the QCD phase diagram.

In this work, we build up a bottom-up holographic model to confront the most
recent lattice results for EOS from the different lattice QCD collaboration
and offer a reliable first-order transition line and CEP in the QCD phase
diagram. By computing thermodynamic quantities via holographic
renormalization, the EOS is found to be quantitatively matched with the latest
lattice QCD simulation. The QCD first-order transition line is fixed from the
free energy and the corresponding CEP in the $T-\mu_{B}$ plane is predicted at
$(T_{C}=109\text{MeV},\mu_{C}=552\text{MeV})$. The CEP in the 2+1 flavor
holographic QCD model, calibrated using data from the Wuppertal-Budapest
lattice QCD group, aligns closely with predictions from a separate holographic
model established by the HotQCD collaboration [30]. This proximity is also
mirrored in the Bayesian analysis results on the location of the QCD critical
point from a holographic perspective [43], indicating that the location of the
CEP exhibits low sensitivity to variations in lattice QCD data inputs.
Additionally, the associated critical exponents approximate those expected
from mean-field theory. Contrarily, as demonstrated by [36], the holographic
2-flavor CEP diverges significantly from mean-field predictions, underscoring
the distinctiveness of our holographic approach compared to conventional large
N QCD models. Our development of diverse holographic models, tailored to mimic
specific QCD systems (e.g., 2+1 flavor[30], pure gluon [37] and 2-flavor [36])
based on corresponding lattice QCD data, enables precise microscopic
descriptions. The similarity in critical exponents suggests that the
2+1-flavor model may belong to the same critical class as the mean-field
theory. Finally, we investigate the behavior of renormalized chiral condensate
and subtracted chiral condensate which is consistent with the WB lattice QCD
simulation [46].

Dedicating to a precise characterization of the properties and differences of
the two phases along the first-order phase transition is an interesting
direction for further study. In the current investigation, we have set up the
preliminary hQCD model to quantitatively study the phase transition in the QCD
phase diagram. The current model captures the characteristic confining
properties and many other relevant physical quantities should be taken into
account to complete the phase diagram, including the chiral condensation and
the transport properties in QGP and hadron gas. In addition, the present
results, in particular those regarding the CEP, should be embedded into the
framework of a more general and multidimensional view of the QCD phase
diagram, including an external magnetic field, an isospin chemical potential,
and a rotation. It will be interesting to consider the real-time dynamics in
our hQCD model far from equilibrium.

##  7 Acknowledgments

This work is supported by the NSFC Grants No.12075101 and No.12235016,
No.12122513, No.12075298, and No.12047503. S.H. would also appreciate
financial support from the Fundamental Research Funds for the Central
Universities and Max Planck Partner Group. Z.L. acknowledges support from the
National Key Program for Science and Technology Research Development
(2023YFB3002500) and the Natural Science Foundation of Henan Province of China
Grant No. 242300420235.

## References

  * [1] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, “The QCD equation of state with dynamical quarks,” JHEP 11 (2010) 077, arXiv:1007.2580 [hep-lat].
  * [2] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, “Full result for the QCD equation of state with 2+1 flavors,” Phys. Lett. B 730 (2014) 99–104, arXiv:1309.5258 [hep-lat].
  * [3] HotQCD Collaboration, A. Bazavov et al., “Equation of state in ( 2+1 )-flavor QCD,” Phys. Rev. D 90 (2014) 094503, arXiv:1407.6387 [hep-lat].
  * [4] X.-y. Xin, S.-x. Qin, and Y.-x. Liu, “Quark number fluctuations at finite temperature and finite chemical potential via the Dyson-Schwinger equation approach,” Phys. Rev. D 90 no. 7, (2014) 076006, arXiv:2109.09935 [hep-ph].
  * [5] F. Gao and Y.-x. Liu, “QCD phase transitions via a refined truncation of Dyson-Schwinger equations,” Phys. Rev. D 94 no. 7, (2016) 076009, arXiv:1607.01675 [hep-ph].
  * [6] S.-x. Qin, L. Chang, H. Chen, Y.-x. Liu, and C. D. Roberts, “Phase diagram and critical endpoint for strongly-interacting quarks,” Phys. Rev. Lett. 106 (2011) 172301, arXiv:1011.2876 [nucl-th].
  * [7] C. Shi, Y.-L. Wang, Y. Jiang, Z.-F. Cui, and H.-S. Zong, “Locate QCD Critical End Point in a Continuum Model Study,” JHEP 07 (2014) 014, arXiv:1403.3797 [hep-ph].
  * [8] C. S. Fischer, J. Luecker, and C. A. Welzbacher, “Phase structure of three and four flavor QCD,” Phys. Rev. D 90 no. 3, (2014) 034022, arXiv:1405.4762 [hep-ph].
  * [9] F. Gao and J. M. Pawlowski, “QCD phase structure from functional methods,” Phys. Rev. D 102 no. 3, (2020) 034027, arXiv:2002.07500 [hep-ph].
  * [10] M. Asakawa and K. Yazaki, “Chiral Restoration at Finite Density and Temperature,” Nucl. Phys. A 504 (1989) 668–684.
  * [11] T. M. Schwarz, S. P. Klevansky, and G. Papp, “The Phase diagram and bulk thermodynamical quantities in the NJL model at finite temperature and density,” Phys. Rev. C 60 (1999) 055205, arXiv:nucl-th/9903048.
  * [12] Z. Li, K. Xu, X. Wang, and M. Huang, “The kurtosis of net baryon number fluctuations from a realistic Polyakov–Nambu–Jona-Lasinio model along the experimental freeze-out line,” Eur. Phys. J. C 79 no. 3, (2019) 245, arXiv:1801.09215 [hep-ph].
  * [13] P. Zhuang, M. Huang, and Z. Yang, “Density effect on hadronization of a quark plasma,” Phys. Rev. C 62 (2000) 054901, arXiv:nucl-th/0008043.
  * [14] W.-j. Fu, J. M. Pawlowski, and F. Rennecke, “QCD phase structure at finite temperature and density,” Phys. Rev. D 101 no. 5, (2020) 054032, arXiv:1909.02991 [hep-ph].
  * [15] H. Zhang, D. Hou, T. Kojo, and B. Qin, “Functional renormalization group study of the quark-meson model with $\omega$ meson,” Phys. Rev. D 96 no. 11, (2017) 114029, arXiv:1709.05654 [hep-ph].
  * [16] W.-j. Fu, X. Luo, J. M. Pawlowski, F. Rennecke, R. Wen, and S. Yin, “Hyper-order baryon number fluctuations at finite temperature and density,” Phys. Rev. D 104 no. 9, (2021) 094047, arXiv:2101.06035 [hep-ph].
  * [17] V. Vovchenko, J. Steinheimer, O. Philipsen, and H. Stoecker, “Cluster Expansion Model for QCD Baryon Number Fluctuations: No Phase Transition at $\mu_{B}/T<\pi$,” Phys. Rev. D 97 no. 11, (2018) 114030, arXiv:1711.01261 [hep-ph].
  * [18] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, “QCD Crossover at Finite Chemical Potential from Lattice Simulations,” Phys. Rev. Lett. 125 no. 5, (2020) 052001, arXiv:2002.02821 [hep-lat].
  * [19] A. Bazavov et al., “Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data,” Phys. Rev. D 101 no. 7, (2020) 074502, arXiv:2001.08530 [hep-lat].
  * [20] S. Borsányi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pásztor, C. Ratti, and K. K. Szabó, “Lattice QCD equation of state at finite chemical potential from an alternative expansion scheme,” Phys. Rev. Lett. 126 no. 23, (2021) 232001, arXiv:2102.06660 [hep-lat].
  * [21] D. Bollweg, D. A. Clarke, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee, P. Petreczky, C. Schmidt, and S. Sharma, “Equation of state and speed of sound of (2+1)-flavor QCD in strangeness-neutral matter at non-vanishing net baryon-number density,” arXiv:2212.09043 [hep-lat].
  * [22] O. Philipsen, “Lattice Constraints on the QCD Chiral Phase Transition at Finite Temperature and Baryon Density,” Symmetry 13 no. 11, (2021) 2079, arXiv:2111.03590 [hep-lat].
  * [23] O. DeWolfe, S. S. Gubser, and C. Rosen, “A holographic critical point,” Phys. Rev. D 83 (2011) 086005, arXiv:1012.1864 [hep-th].
  * [24] O. DeWolfe, S. S. Gubser, and C. Rosen, “Dynamic critical phenomena at a holographic critical point,” Phys. Rev. D 84 (2011) 126014, arXiv:1108.2029 [hep-th].
  * [25] Y. Chen, D. Li, and M. Huang, “The dynamical holographic QCD method for hadron physics and QCD matter,” Commun. Theor. Phys. 74 no. 9, (2022) 097201, arXiv:2206.00917 [hep-ph].
  * [26] R. Rougemont, J. Grefa, M. Hippert, J. Noronha, J. Noronha-Hostler, I. Portillo, and C. Ratti, “Hot QCD Phase Diagram From Holographic Einstein-Maxwell-Dilaton Models,” arXiv:2307.03885 [nucl-th].
  * [27] M. Jarvinen, “Holographic dense QCD in the Veneziano limit,” EPJ Web Conf. 274 (2022) 08006, arXiv:2211.10005 [hep-ph].
  * [28] J. Knaute, R. Yaresko, and B. Kämpfer, “Holographic QCD phase diagram with critical point from Einstein–Maxwell-dilaton dynamics,” Phys. Lett. B 778 (2018) 419–425, arXiv:1702.06731 [hep-ph].
  * [29] R. Critelli, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, “Critical point in the phase diagram of primordial quark-gluon matter from black hole physics,” Phys. Rev. D 96 no. 9, (2017) 096026, arXiv:1706.00455 [nucl-th].
  * [30] R.-G. Cai, S. He, L. Li, and Y.-X. Wang, “Probing QCD critical point and induced gravitational wave by black hole physics,” Phys. Rev. D 106 no. 12, (2022) L121902, arXiv:2201.02004 [hep-th].
  * [31] X. Chen and M. Huang, “Machine learning holographic black hole from lattice QCD equation of state,” arXiv:2401.06417 [hep-ph].
  * [32] X.-Y. Liu, X.-C. Peng, Y.-L. Wu, and Z. Fang, “Holographic study on QCD phase transition and phase diagram with two flavors,” Phys. Rev. D 109 no. 5, (2024) 054032, arXiv:2312.01346 [hep-ph].
  * [33] Y.-Q. Zhao, S. He, D. Hou, L. Li, and Z. Li, “Phase diagram of holographic thermal dense QCD matter with rotation,” JHEP 04 (2023) 115, arXiv:2212.14662 [hep-ph].
  * [34] Z. Li, J. Liang, S. He, and L. Li, “Holographic study of higher-order baryon number susceptibilities at finite temperature and density,” Phys. Rev. D 108 no. 4, (2023) 046008, arXiv:2305.13874 [hep-ph].
  * [35] S. He, L. Li, S. Wang, and S.-J. Wang, “Constraints on holographic QCD phase transitions from PTA observations,” arXiv:2308.07257 [hep-ph].
  * [36] Y.-Q. Zhao, S. He, D. Hou, L. Li, and Z. Li, “Phase structure and critical phenomena in 2-flavor QCD by holography,” arXiv:2310.13432 [hep-ph].
  * [37] S. He, L. Li, Z. Li, and S.-J. Wang, “Gravitational waves and primordial black hole productions from gluodynamics by holography,” Sci. China Phys. Mech. Astron. 67 no. 4, (2024) 240411, arXiv:2210.14094 [hep-ph].
  * [38] S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo, and C. Ratti, “Higher order fluctuations and correlations of conserved charges from lattice QCD,” JHEP 10 (2018) 205, arXiv:1805.04445 [hep-lat].
  * [39] L. Li, “On Thermodynamics of AdS Black Holes with Scalar Hair,” Phys. Lett. B 815 (2021) 136123, arXiv:2008.05597 [gr-qc].
  * [40] S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, “Freeze-out parameters: lattice meets experiment,” Phys. Rev. Lett. 111 (2013) 062005, arXiv:1305.5161 [hep-lat].
  * [41] J. Grefa, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, “Hot and dense quark-gluon plasma thermodynamics from holographic black holes,” Phys. Rev. D 104 no. 3, (2021) 034002, arXiv:2102.12042 [nucl-th].
  * [42] H.-w. Zheng, Y. Lu, F. Gao, S.-x. Qin, and Y.-x. Liu, “The effective potential of composite operator in the first order region of QCD phase transition,” arXiv:2312.00382 [hep-ph].
  * [43] M. Hippert, J. Grefa, T. A. Manning, J. Noronha, J. Noronha-Hostler, I. Portillo Vazquez, C. Ratti, R. Rougemont, and M. Trujillo, “Bayesian location of the QCD critical point from a holographic perspective,” arXiv:2309.00579 [nucl-th].
  * [44] N. Goldenfeld, Lectures on phase transitions and the renormalization group.  1992\.
  * [45] A. Karch, E. Katz, D. T. Son, and M. A. Stephanov, “Linear confinement and AdS/QCD,” Phys. Rev. D 74 (2006) 015005, arXiv:hep-ph/0602229.
  * [46] Wuppertal-Budapest Collaboration, S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, “Is there still any $T_{c}$ mystery in lattice QCD? Results with physical masses in the continuum limit III,” JHEP 09 (2010) 073, arXiv:1005.3508 [hep-lat].


# First-principles Calculations of Raman and Infrared Spectroscopy For Phase
Identification and Strain Calibration of Hafnia.

Aldo Raeliarijaona<EMAIL_ADDRESS>R. E. Cohen
<EMAIL_ADDRESS>Extreme Materials Initiative, Earth and Planets
Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW,
Washington, DC 20015, USA

###### Abstract

Using density functional perturbation theory (DFPT) we computed the phonon
frequencies, Raman and IR activities of hafnia polymorphs (P42nmc , Pca21 ,
Pmn21 , Pbca OI, brookite, and baddeleyite) for phase identification. We
investigated the evolution of Raman and IR activities with respect to
epitaxial strain and provide plots of frequency differences as a function of
strain for experimental calibration and identification of the strain state of
the sample. We found Raman signatures of different hafnia polymorphs:
$\omega(A_{1g})=300$ cm-1 for P42nmc , $\omega(A_{1})=343$ cm-1 for Pca21 ,
$\omega(B_{2})=693$ cm-1 for Pmn21 , $\omega(A_{g})=513$ cm-1 for Pbca (OI),
$\omega(A_{g})=384$ cm-1 for brookite, and $\omega(A_{g})=496$ cm-1 for
baddeleyite. We also identified the Raman $B_{1g}$ mode, an anti-phase
vibration of dipole moments, ( $\omega(B_{1g})=758$ cm-1 for OI,
$\omega(B_{1g})=784$ cm-1 for brookite) as the Raman signature of antipolar
Pbca structures. We calculated a large splitting between longitudinal optical
(LO) and transverse optical (TO) modes ($\Delta{\omega_{\text{LO-
TO}}(A^{z}_{1})}=255$ cm-1 in Pca21 , and $\Delta{\omega_{\text{LO-
TO}}(A_{1})}=263$ cm-1 in Pmn21 ) to the same order as those observed in
perovskite ferroelectrics, and related them to the anomalously large Born
effective charges of Hf atoms ($Z^{*}(\text{Hf})=5.54$).

Hafnia is of great interest because of its compatibility with silicon and its
robust ferroelectricity at the nanoscale. It has a high permittivity and has
been proposed as a dielectric material for dynamic random access memories
(DRAM) and logic devices[1, 2]. At ambient pressure and temperature it
stabilizes in the monoclinic crystal structure (P21/c) [3], but can also have
different structures depending on temperature and pressure. At ambient
pressure there are two high-temperature polymorphs: the tetragonal (P42nmc ),
when 1700∘C $<$ T $<$ 2600∘C, and the cubic ($Fm\bar{3}m$) when T
$>$2600∘C[4]. Other phases have been stabilized under different pressure, such
as Pbca, when 5.9 GPa $<$ P $<$ 15 GPa, and Pnma when under P $>$ 15 GPa[4,
5], epitaxial strain[6], or doping[7, 8].

Recent investigations have highlighted the possibility and importance of
ferroelectricity in hafnia[9, 10, 11, 12, 13, 14], and studies of yttrium-
doped hafnia revealed the Pca21 polymorph as the polar phase responsible for
ferroelectricity[15]. Antiferroelectricity has also been observed in doped
hafnia [16] and further theoretical studies have shown that the orthorhombic
Pbca structure, which contains two Pca21 cells with opposite oxygen
displacement relative to the hafnia atoms, as the antiferroelectric
orthorhombic structure[11]. Hafnia is also of fundamental interest due to its
fluorite-based structure, different from the wide class of well-studied
perovskite ferroelectrics. One hindrance to the rapid development of hafnia
ferroelectrics has been a lack of basic understanding of ferroelectricity.

Structural studies are usually based on X-ray diffraction (XRD) methods.
Complementary to XRD, inelastic light scattering such as Raman or IR
spectroscopy is also a tool that can be used to determine the structure and
structural response of different hafnia phases. The latter approach reveals
information about interatomic interaction, hence atomic arrangements, through
the vibrational spectra. Epitaxial strain is important in thin films and it
can affect the Raman frequencies and intensities. One of the advantages of
Raman spectroscopy is that it can provide information about the strain on the
sample once calibrated. Only Raman shifts and peak intensities are accessible
experimentally thus the need of first-principles studies of atomic vibrations
for symmetry assignment of Raman peaks. Although Raman and IR data are
available from previous calculations of hafnia polymorphs[17, 18, 19] they do
not cover all the relevant phases, namely both polar orthorhombic phases or
both antipolar orthorhombic phases. Furthermore, experimentally it is normal
to observe mixture of different phases in hafnia such as baddeleyite and Pca21
[20], which makes the symmetry determination using Raman or IR spectroscopy
crucial.

Our aim in this paper is to provide tools for experiments to identify the
structure of hafnia and determine the epitaxial strain to which the sample is
subject. Using DFPT we computed the phonon frequencies and their activities
(Raman and IR) at different strain values
($\eta=[-3.0\%,-1.5\%,0.0\%,1.5\%,3.0\%]$) for different hafnia polymorphs,
namely baddeleyite, P42nmc , Pmn21 , Pca21 , Pbca (OI), and brookite. We
tracked the evolution of phonon mode frequencies under strain, and provided
means to identify the strain through identification of strain-induced
frequency shifts of signature phonons.

We performed first-principles calculations using Quantum Espresso [21, 22,
23], with optimized norm-conserving Vanderbilt (ONCV) pseudopotentials [24].
The cell parameters and atomic positions were optimized using the local
density approximation (LDA) with the Perdew-Zunger (PZ) exchange-correlation
functional [25], generated using the oncvpsp code [24]. The plane-wave
expansion is truncated using a cutoff energy of E${}_{\textbf{cutoff}}$ =
1306.6 eV, and the Brillouin zone was sampled using an 8$\times$8$\times$8
Monkhorst-Pack grid [26]. We calculated the Brillouin zone center phonon
frequencies using density functional perturbation theory (DFPT) [27]
implemented in the QE/PH package [21, 22, 23]. We considered the powder Raman
spectra, with the Raman activity computed using Placzek approximation[28, 29,
30]. We plotted the spectra by convoluting with a Gaussian shape function
whose width is 5 cm-1 estimated from experiments. Similar widths have been
used in previous studies[30].

We set a coordinate system similar to the ones defined in Ref.6. We consider
the epitaxy of the tetragonal case with a square substrate such as Yttria-
stabilized zirconia (YSZ), with z normal to the epitaxial plane, and a and b
in registry with the epitaxial surface. Strain $\eta$ was applied by setting:

$a_{\eta}=(1+\eta)a_{0},$ (1)

where $a_{0}$ is the strain-free ground state lattice constant of the
tetragonal phase. The out-of-plane lattice constant ($c$) and the atomic
positions were relaxed at each strain.

To validate the accuracy of our DFPT calculations, we compare the calculated
Raman spectra of the strain-free baddeleyite and brookite to the measured
Raman spectra[31] (Fig.1). The Raman spectrum of brookite[31] was measured
from a sample under hydrostatic pressure of P = 5.9 GPa. So, for fair
comparison, the cell parameters and the atomic positions for our calculations
for brookite were relaxed under the same pressure. The calculated spectra of
baddeleyite and brookite match the experimentally measured spectra well
(Fig.1). For baddeleyite, the high intensity $A_{g}$ peaks
($\omega(A_{g})=495$ cm-1 ), the shoulder peak $B_{g}$ ($\omega(B_{g})=510$
cm-1 ), and the double peaks $B_{g}$ ($\omega(B_{g})=637$ cm-1 ) and $A_{g}$
($\omega(A_{g})=670$ cm-1 ) all closely match with experimental peaks; for
brookite we can point to the high intensity $A_{g}$ peak ($\omega(A_{g})=384$
cm-1 ), the $B_{1g}$ peak ($\omega(B_{1g})=513$ cm-1 ), the double peaks
$B_{3g}$ ($\omega(B_{3g})=593$ cm-1 ) and $B_{1g}$ ($\omega(B_{g})=604$ cm-1
). Our calculated spectra for baddeleyite, tetragonal and brookite also
compare well with other calculations[17, 18] (see Supplementary Material
Fig.S1).

Figure 1: Raman spectra of (a) baddeleyite,(b) brookite compared to
experimental and theoretical Raman. Gray and green spectra are experimental
data from Ref.31 and the red curve is a theoretical spectrum from Ref.17. We
convoluted the spectra with a Gaussian shape function of width 5 cm-1
estimated from experiments. The squares indicate modes chosen to compute the
frequency differences (Fig.4). For brookite Raman spectra (b), the blue and
green curves were obtained under hydrostatic pressure (P = 5.9 GPa).

We plot the Raman spectra of the different strain-free hafnia polymorphs
(Fig.1,Fig.2) to help with the identification of the crystal structures. The
high intensity peaks of the Raman spectra can help in the identification of
the structure symmetry. For example, the high intensity mode at
$\omega(A_{g})=496$ cm-1 is found only in baddeleyite. The Pmn21 structure can
be identified using the high frequency $B_{2}$ ($\omega(B_{2})=693$ cm-1 ) and
$A_{1}$ ($\omega(A_{1})=722$ cm-1 ) peaks. The Pca21 structure can be
identified using the $A_{1}$ peak ($\omega=340$ cm-1 ) and the strong $A_{1}$
peak ($\omega=670$ cm-1 ). The presence of these peaks can be used to identify
the Pca21 structure. This is an advantage of the Raman (or IR) spectroscopy
over the XRD given that phase identification using XRD can be harder due to
the monoclinic and Pca21 phase mixtures and the overlap of XRD peaks [20],
whereas with Raman spectroscopy the difference can be resolved easily. In Pbca
(OI) the signature modes are the strong $A_{g}$ peaks at $\omega(A_{g})=513$
cm-1 , and the double $A_{g}$, and $B_{2g}$ peaks at $\omega(A_{g})=658$ cm-1
and $\omega(B_{2g})=668$ cm-1 , respectively. Finally, in brookite the $A_{g}$
mode ($\omega(A_{g})=384$ cm-1 ) can be used as its fingerprint.

Figure 2: Raman spectra of (a) P42nmc , (b) Pmn21 , (c) Pca21 , and (d) Pbca
(OI). We convoluted the spectra using a Gaussian shape function of width 5
cm-1 estimated from experiments. The squares indicate modes chosen to compute
the frequency differences (Fig.4).

Additionally, we find a signature mode for antipolar polymorphs. This is the
high frequency B1g peak ($\omega(B_{1g})=811$ cm-1 for brookite, or
$\omega(B_{1g})=791$ cm-1 for Pbca OI). This mode is absent in the FE (Pca21
and Pmn21 ) structures and it is an anti-phase vibration of the O atoms in
neighboring cells of the Pbca (Fig.3). This anti-phase vibration of the O
atoms can be thought as an anti-phase collective motion of the dipole moments.

Figure 3: The atomic displacements corresponding to the high frequency
$B_{1g}$ mode in Pbca AFE structures. The larger spheres represent the Hf
atoms and the smaller spheres the O atoms. The blue arrows on the atoms
indicate the direction of the atomic vibrations. The arrows on the bottom
right corner show the crystallographic directions.

The splitting between longitudinal optical (LO) and transverse optical (TO)
phonon modes is a subject of interest in ferroelectrics. Indeed, perovskite
ferroelectrics are known to have phonon modes that exhibit giant LO-TO
splittings[32]. Using the rigorous definition of LO-TO splitting[33], we show
that these large LO-TO splittings also occur in hafnia. In polar Pmn21 phase,
the calculated LO-TO splitting is $\Delta{\omega_{\text{LO-TO}}}(A_{1})=263$
cm-1 for the high frequency $A_{1}$ mode ($\omega(A_{1})=722$ cm-1 ); for the
Pca21 phase we computed the LO-TO splitting to be $\Delta{\omega_{\text{LO-
TO}}(A^{z}_{1})}=255$ cm-1 for $\omega(A^{z}_{1})=676$ cm-1 and
$\Delta{\omega_{\text{LO-TO}}(B^{x}_{2})}=146$ cm-1 for the $B^{x}_{2}$ mode
at $\omega(B^{x}_{2})=722$ cm-1 (Fig.2). These large LO-TO splittings can be
attributed to the unusually large Born effective charge ($Z^{*}$) of Hf atoms
in hafnia (see Supplementary Material). Although the born effective charges
$Z^{*}(O)$ are close to the nominal charge of O atoms (-2), with few
exceptions such as: $Z^{*}_{yy}(O1)=Z^{*}_{xx}(O2)=-3.42$ in P42nmc ,
$Z^{*}_{xx}(O1)=-3.15$ and $Z^{*}_{zz}(O2)=-3.08$ in Pmn21 , the Hf atoms have
anomalously large dynamical charges up to 5.54. Other calculations also
reported similar values[34, 19]. The three modes with large LO-TO splitting
mentioned earlier all have in common the antiphase motion between the Hf and
O.

The constraint on the lattice constant imposed by the substrate affects the
interatomic distances in the sample, which not only shifts the phonon
frequencies but also changes the Raman and IR intensities. The shift in Raman
or IR frequencies thus indicates the strains in a hafnia film. The evolution
of phonon modes with respect to strain can be tracked using correspondence
between the phonon eigenvectors at each strain following the approach in
Ref.33. Firstly, the phonon eigenvectors at strain $\eta_{2}$ were projected
to the eigenvectors at strain $\eta_{1}$:

$\epsilon_{n}^{i\alpha}(\eta_{2})=\sum\limits_{m}a_{mn}\epsilon_{m}^{i\alpha}(\eta_{1}),$
(2)

where $m$ and $n$ are the phonon mode indices, $i$ and $\alpha$ are the atoms
and the direction indices respectively; the coefficient
$a_{mn}=\sum\limits_{i\alpha}\langle\epsilon_{m}^{i\alpha}(\eta_{1})|\epsilon_{n}^{i\alpha}(\eta_{2})\rangle$
indicates the correlation or projection of phonon eigenmode $n$ at strain
$\eta_{2}$ to the eigenmodes at strain $\eta_{1}$. Then by choosing the mode
with the maximum coefficient $a_{mn}$ (Eq. 2), i.e. the mode with the highest
correspondence probability to mode $n$ denoted by $\omega^{m^{\prime}}_{\eta
1}$, we assign $\omega^{n}_{\eta 2}\rightarrow\omega^{m^{\prime}}_{\eta 1}$.

The calculation of the frequency difference is straightforward once equipped
with the strain-to-strain mode correspondence (Eq.2):

$\Delta{\omega_{l,n}(\eta)}=\omega_{l^{\prime}}(\eta)-\omega_{n^{\prime}}(\eta),$
(3)

where it should be understood that $\omega_{l^{\prime}}(\eta)$, and
$\omega_{n^{\prime}}(\eta)$ are the phonon modes that correspond to
$\omega_{l}(0)$, and $\omega_{n^{\prime}}(0)$ at zero strain, respectively.
$\Delta{\omega_{l,n}}(\eta)$ changing sign means that the order of the modes
considered switched.

Figure 4: Evolution of frequency difference with respect to epitaxial strain
for selected modes of (a) baddeleyite, (b) P42nmc , (c) Pmn21 , (d) Pca21 ,
(e) Pbca OI, and f brookite. The star, diamond, and triangle symbols in (a),
(b), and (f) are frequency differences using data from Ref.31, Ref.17, and
Ref.35, respectively. The blue and red curves are calibration using IR and
Raman active modes, respectively. The Raman modes used for calibrations are
marked by the empty squares on the Raman spectra in Fig1 and the IR modes are
marked as blue modes in Table S3–S8 (see Supplementary Material).

The evolution of the phonon modes is not straightforward, and choices had to
be made on the modes used to compute the frequency difference (Fig.4). The
modes chosen (Fig. 1, and Fig.2) were modes with non-negligible Raman or IR
intensities and the ones showing monotonic behavior if possible, within a
range of 200 cm-1 for plotting purposes. In Fig.4 we used results from studies
with samples that are either under hydrostatic pressure[31] or doped
hafnia[35] for comparison. For the case where doping helps with the
stabilization of a particular polymorph, peaks associated to defects or
dopants were identified[35] so the comparison with the remaining peaks,
probably from the tetragonal polymorph, with our calculation results can be
made.

For the comparison of Pbca brookite spectra, the Raman spectra in Ref.31 were
obtained under hydrostatic pressure. We performed two different first-
principles calculations for the brookite phase, one at ambient pressure and
one under 5.9 GPa. We find that pressure hardens the individual Raman modes,
as expected. But the difference between two modes ($\Delta{\omega_{l,n}}$) for
the system at ambient pressure was similar to $\Delta{\omega_{l,n}}$ for the
system under hydrostatic pressure. It is thus fair to compare
$\Delta{\omega_{l,n}}$ from the experiment under pressure with
$\Delta{\omega_{l,n}}$ from DFPT at ambient pressure (Fig.4f).

For most cases, the frequency difference between two peaks changes linearly
with strain, such as $\Delta{\omega_{B_{1g},B_{2g}}}(\eta)=80+4.2\eta$ for the
Raman active modes chosen in brookite or
$\Delta{\omega_{A_{2},B_{2}}}(\eta)=89+5.9\eta$ for the Raman active modes
chosen in Pmn21 . On the other hand, for the Pca21 and baddeleyite the
frequency difference is better described by a quadratic polynomial of the
strain $\eta$:

$\Delta{\omega_{A_{g},B_{g}}}(\eta)=142-13.7\eta+2\eta^{2},$ (4)

for the Raman active modes chosen in baddeleyite, and

$\Delta{\omega_{A_{1},B_{1}}}(\eta)=15+14.1\eta+3.7\eta^{2}$ (5)

for Pca21 .

Now we briefly discuss an interesting comparison between hafnia and zirconia,
two isomorphous compounds[3]. Recent studies on ZrO2 combined experimental
measurements and computations of Raman spectra to identify the zirconia
polymorphs using Raman[36]. The Raman spectra of monoclinic hafnia and
zirconia are similar, e.g. the $A_{g}$ peak at $\omega(A_{g})=495$ cm-1 in
hafnia and the $A_{g}$ peak at $\omega(A_{g})=470$ cm-1 in zirconia, or the
$B_{g}$ peak at $\omega(B_{g})=637$ cm-1 in hafnia compared to
$\omega(B_{g})=598$ cm-1 in zirconia. The quantitative mismatch between the
spectra in hafnia and zirconia (25 cm-1 for $A_{g}$ and 30 cm-1 for Bg) can be
attributed to differences in the interatomic force constant (IFC) resulting
from the difference in the cation. Interestingly, despite these quantitative
discrepancies in the Raman frequencies of hafnia and zirconia, the frequency
differences $\Delta{\omega_{l,n}}$ are comparable. In monoclinic hafnia for
example, the frequency difference $\Delta{\omega_{A_{g},B_{g}}}=142$ cm-1
compares well with $\Delta{\omega_{A_{g},B_{g}}}=140$ cm-1 in zirconia, and in
tetragonal the frequency difference $\Delta{\omega_{E_{g},B_{1g}}}=130$ cm-1
in hafnia is a close match to $\Delta{\omega_{E_{g},B_{1g}}}=135$ cm-1 in
zirconia.

We computed the phonon frequencies of P42nmc , Pmn21 , Pca21 , Pbca (OI),
brookite, and baddeleyite using DFPT at different values of epitaxial strains.
We found that Raman spectra can be used to identify the symmetry or phases of
hafnia, namely $\omega(A_{1g})=300$ cm-1 for P42nmc , $\omega(A_{1})=343$ cm-1
for Pca21 , $\omega(B_{2})=693$ cm-1 for Pmn21 , $\omega(A_{g})=513$ cm-1 for
Pbca (OI), $\omega(A_{g})=384$ cm-1 for brookite, and $\omega(A_{g})=496$ cm-1
for baddeleyite. We also identified the Raman signature of AFE structures as
the $\omega(B_{1g})=784$ cm-1 peak for brookite or the $\omega(B_{1g})=758$
cm-1 peak for Pbca (OI). The presence (or absence) of this Raman signal can be
used to distinguish between AFE and FE orthorhombic hafnia. Further, we showed
the evolution of frequency differences between selected normal mode
frequencies with respect to strain for calibration purposes and identification
of strain state of the hafnia crystal.

## Supplementary Material

Comparison of bookite Raman spectra at ambient pressure(Fig.S1), the
calculated Born effective charges (Table S1–S2), the infrared frequencies and
IR activities of different hafnia polymorphs (Table S3–S8).

###### Acknowledgements.

The authors thank Pavlo Zubko for helpful discussions. This work is supported
by U. S. Office of Naval Research Grants No. N00014-17-1-2768 and
N00014-20-1-2699, and the Carnegie Institution for Science. Computations were
supported by DOD HPC, Carnegie computational resources, and REC gratefully
acknowledges the Gauss Centre for Supercomputing e.V. (www.gauss-center.eu)
for funding this project by providing computing time on the GCS Supercomputer
SuperMUC-NG at Leibniz Supercomputing Centre (LRZ, www.lrz.de). Some of the
computing for this project was performed on the Memex cluster. We would like
to thank Carnegie Institution for Science and the Carnegie Sci-Comp Committee
for providing computational resources and support that contributed to these
research results.

## References

  * Kingon, Maria, and Streiffer [2000] A. I. Kingon, J.-P. Maria, and S. K. Streiffer, “Alternative dielectrics to silicon dioxide for memory and logic devices,” Nature 406, 1032–1038 (2000).
  * Wilk, Wallace, and Anthony [2001] G. D. Wilk, R. M. Wallace, and J. M. Anthony, “High- gate dielectrics: Current status and materials properties considerations,” Journal of Applied Physics 89, 5243–5275 (2001).
  * Ruh and Corfield [1970] R. Ruh and P. W. R. Corfield, “Crystal Structure of Monoclinic Hafnia and Comparison with Monoclinic Zirconia,” Journal of the American Ceramic Society 53, 126–129 (1970).
  * Ohtaka  _et al._ [2001] O. Ohtaka, H. Fukui, T. Kunisada, T. Fujisawa, K. Funakoshi, W. Utsumi, T. Irifune, K. Kuroda, and T. Kikegawa, “Phase relations and volume changes of hafnia under high pressure and high temperature,” Journal of the American Ceramic Society 84, 1369–1373 (2001).
  * Huan  _et al._ [2014] T. D. Huan, V. Sharma, G. A. Rossetti, and R. Ramprasad, “Pathways Towards Ferroelectricity in Hafnia,” Phys. Rev. B 90, 064111 (2014).
  * Liu and Hanrahan [2019] S. Liu and B. M. Hanrahan, “Effects of Growth Orientations and Epitaxial Strains on Phase Stability of HfO2 Thin Films,” Physical Review Materials 3, 054404 (2019).
  * Park  _et al._ [2017] M. H. Park, T. Schenk, C. M. Fancher, E. D. Grimley, C. Zhou, C. Richter, J. M. LeBeau, J. L. Jones, T. Mikolajick, and U. Schroeder, “A comprehensive study on the structural evolution of HfO2 thin films doped with various dopants,” J. Mater. Chem. C 5, 4677–4690 (2017).
  * Batra  _et al._ [2017] R. Batra, T. D. Huan, G. A. Rossetti, and R. Ramprasad, “Dopants Promoting Ferroelectricity in Hafnia: Insights from a comprehensive Chemical Space Exploration,” Chemistry of Materials 29, 9102–9109 (2017).
  * Böscke  _et al._ [2011] T. S. Böscke, J. Müller, D. Bräuhaus, U. Schröder, and U. Böttger, “Ferroelectricity in Hafnium Oxide Thin Films,” Applied Physics Letters 99, 102903 (2011).
  * Schroeder, Hwang, and Funakubo [2019] U. Schroeder, C. Hwang, and H. Funakubo,  _Ferroelectricity in Doped Hafnium Oxide: Materials, Properties and Devices_, Woodhead Publishing Series in Electronic and Optical Materials (Elsevier Science, 2019).
  * Lee  _et al._ [2020] H.-J. Lee, M. Lee, K. Lee, J. Jo, H. Yang, Y. Kim, S. C. Chae, U. Waghmare, and J. H. Lee, “Scale-free Ferroelectricity Induced by Flat Phonon Bands in HfO2,” Science 369, 1343–1347 (2020).
  * Qi  _et al._ [2020] Y. Qi, S. Singh, C. Lau, F.-T. Huang, X. Xu, F. J. Walker, C. H. Ahn, S.-W. Cheong, and K. M. Rabe, “Stabilization of Competing Ferroelectric Phases of HfO2 Under Epitaxial Strain,” Phys. Rev. Lett. 125, 257603 (2020).
  * Xu  _et al._ [2021] X. Xu, F.-T. Huang, Y. Qi, S. Singh, K. M. Rabe, D. Obeysekera, J. Yang, M.-W. Chu, and S.-W. Cheong, “Kinetically Stabilized Ferroelectricity in Bulk Single-Crystalline HfO2:Y,” Nature Materials 20, 826–832 (2021).
  * Mikolajick and Schroeder [2021] T. Mikolajick and U. Schroeder, “Ferroelectricity in Bulk Hafnia,” Nature Materials 20, 718–719 (2021).
  * Müller  _et al._ [2011] J. Müller, U. Schröder, T. S. Böscke, I. Müller, U. Böttger, L. Wilde, J. Sundqvist, M. Lemberger, P. Kücher, T. Mikolajick, and L. Frey, “Ferroelectricity in yttrium-doped hafnium oxide,” Journal of Applied Physics 110, 114113 (2011).
  * Park  _et al._ [2015] M. H. Park, Y. H. Lee, H. J. Kim, Y. J. Kim, T. Moon, K. D. Kim, J. Müller, A. Kersch, U. Schroeder, T. Mikolajick, and C. S. Hwang, “Ferroelectricity and antiferroelectricity of doped thin hfo2-based films,” Advanced Materials 27, 1811–1831 (2015).
  * Zhou  _et al._ [2014] B. Zhou, H. Shi, X. D. Zhang, Q. Su, and Z. Y. Jiang, “The Simulated Vibrational Spectra of HfO2 Polymorphs,” Journal of Physics D: Applied Physics 47, 115502 (2014).
  * Gunst, Stradi, and Blom [2020] T. Gunst, D. Stradi, and A. Blom, “Identification of zirconia and hafnia crystalline phases by optical spectroscopy from first-principles,” in  _Nanoengineering: Fabrication, Properties, Optics, Thin Films, and Devices XVII_, Vol. 11467, edited by B. Panchapakesan, A.-J. Attias, and W. Park, International Society for Optics and Photonics (SPIE, 2020) pp. 34 – 47\.
  * [19] S. Fan, S. Singh, X. Xu, K. Park, Y. Qi, a. V. S.W. Cheong, K. M. Rabe, and J. L. Musfeldt, “Vibrational fingerprints of ferroelectric hafnia,” arXiv:2201.12643 .
  * Mukundan  _et al._ [2021] V. Mukundan, S. Consiglio, D. H. Triyoso, K. Tapily, M. E. McBriarty, S. Schujman, K. Beckmann, J. Hazra, V. Kaushik, N. Cady, R. D. Clark, G. J. Leusink, and A. C. Diebold, “Ferroelectric Phase Content in 7 nm Hf(1−x)ZrxO2 Thin Films Determined by X-Ray-Based Methods,” physica status solidi (a) 218, 2100024 (2021).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, “QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials,” Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, “Advanced Capabilities for Materials Modelling With QUANTUM ESPRESSO,” Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, “Quantum ESPRESSO Toward the Exascale,” The Journal of Chemical Physics 152, 154105 (2020).
  * Hamann [2013] D. R. Hamann, “Optimized Norm-Conserving Vanderbilt Pseudopotentials,” Phys. Rev. B 88, 085117 (2013).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, “Self-Interaction Correction to Density-Functional Approximations for Many-Electron Systems,” Phys. Rev. B 23, 5048–5079 (1981).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, “Special Points for Brillouin-Zone Integrations,” Phys. Rev. B 13, 5188–5192 (1976).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, “Phonons and Related Crystal Properties from Density-Functional Perturbation Theory,” Rev. Mod. Phys. 73, 515–562 (2001).
  * Porezag and Pederson [1996] D. Porezag and M. R. Pederson, “Infrared intensities and raman-scattering activities within density-functional theory,” Phys. Rev. B 54, 7830–7836 (1996).
  * Lazzeri and Mauri [2003] M. Lazzeri and F. Mauri, “First-Principles Calculation of Vibrational Raman Spectra in Large Systems: Signature of Small Rings in Crystalline ${\mathrm{s}\mathrm{i}\mathrm{o}}_{2}$,” Phys. Rev. Lett. 90, 036401 (2003).
  * Caracas and Cohen [2006] R. Caracas and R. E. Cohen, “Theoretical determination of the Raman spectra of MgSiO3perovskite and post-perovskite at high pressure,” Geophysical Research Letters 33 (2006), https://doi.org/10.1029/2006GL025736.
  * Jayaraman  _et al._ [1993] A. Jayaraman, S. Y. Wang, S. K. Sharma, and L. C. Ming, “Pressure-induced Phase Transformations in HfO2 to 50 GPa studied by Raman spectroscopy,” Physical Review B 48, 9205–9211 (1993).
  * Zhong, King-Smith, and Vanderbilt [1994] W. Zhong, R. D. King-Smith, and D. Vanderbilt, “Giant lo-to splittings in perovskite ferroelectrics,” Phys. Rev. Lett. 72, 3618–3621 (1994).
  * Raeliarijaona and Fu [2015] A. Raeliarijaona and H. Fu, “ModeSequence, Frequency Change of Nonsoft Phonons, and LO-TO Splitting in strained tetragonal ${\mathrm{batio}}_{3}$,” Phys. Rev. B 92, 094303 (2015).
  * Neal  _et al._ [2021] S. N. Neal, S. Li, T. Birol, and J. L. Musfeldt, “Chemical bonding and Born charge in 1T-HfS2,” npj 2D Materials and Applications 5, 45 (2021).
  * Masatomo  _et al._ [1996] Y. Masatomo, T. H., O. Katsuya, H. Teruo, K. Masato, A. Haruo, I. Yasuro, S. Yasuo, and Y. Masahiro, “Formation of Metastable Forms by Quenching of the HfO2RO1.5 Melts (R = Gd, Y and Yb),” Journal of Physics and Chemistry of Solids 57, 289–295 (1996).
  * Materano  _et al._ [2022] M. Materano, P. Reinig, A. Kersch, M. Popov, M. Deluca, T. Mikolajick, U. Boettger, and U. Schroeder, “Raman Spectroscopy as a Key Method to Distinguish the Ferroelectric Orthorhombic Phase in Thin ZrO2-Based Films,” Physica Status Solidi (RRL) – Rapid Research Letters 16, 2100589 (2022).