# Project excited camel

# Project hiding frog, the hidinger stuff.

./Jul_2021_OY

The dielectric response of rock-salt crystals at finite temperatures from
first principles
Supplementary Material

##  SI Experimental setup

All Raman measurements were performed on our home-built system. Figure S1
shows a schematic of the system. The beam path begins with a 488nm solid-state
(Coherent Sapphire SF 488-100 CDRH) laser, which will excite the Raman
scattering in the sample. The laser is filtered for any amplified spontaneous
emission by a volume holographic ASE (Ondax) filters. Next, the beam shape is
optimized and culminated by a pinhole spatial filter. Control over incident
and scattered polarization is realized by a set of two polarizers and two
half-wave plates. First, linear polarization is insured by a (Thorlabds)
calcite polarizer. A monochromatic half-wave plate (HW0 in Fig. S1) is placed
before the polarizer and rotated for optimal incident intensity, correcting
for the original laser polarization. The beam then passes through a
(NoiseBlock 90/10 Ondax) beam splitter (BS) with the experimental polarization
orientation controlled by another monochromatic half-wave plate (HW1 in Fig.
S1). The polarized beam then enters a microscope and focused on the sample by
a 10X (Zeiss) objective. Excitation powers of 7-30 mW were used (depending on
temperature), measured just before the microscope. The sample (see SIII) is
placed inside a (TS1000 HiTemp) Linkam stage under Argon purge. Purging with a
mono-atomic gas proved crucial, since under vacuum local laser heating would
create defects, whisle nitrogen would dominate the Raman signal with its
rotational excitations. A mono-atomic gas is also the best option for
effective thermal coupling to the Linkam hot plate. Inside the Linkam the
magic happens, and laser light is scattered off the thermally excited crystal.
The back-scattered light is collected by the same objective and passes back
through HW1 and through the BS, where 90% of the Rayleigh is eliminated and
Raman signal transmitted. Note that the beam’s polarization at this point is
not known, since we can not assume the scattered light, especially the Raman
signal, scattered back into the same incident polarization, as is described in
Fig. S1 (b). Indeed, there is no reason to assume the beam is still linearly
polarized. We may, however, as always, describe the beam’s polarization as the
sum of two linear components, one parallel to the incident beam polarization,
and the other perpendicular. After passing again through HW1, all light
scattered in parallel to the incident beam polarization will now parallel the
original polarization determined by the first polarizer. All light scattered
perpendicular to the incident beam, will now be also polarized perpendicular
to the first polarizer. After the BS, the scattered beam (orange line in Fig.
S1) passes through an achromatic half-wave plate (HW2 in Fig. S1) which either
leaves or rotates the signal polarization by 900. Another identical polarizer
(designated ”analyzer” in Fig. S1) filters the polarization of the signal. We
thus separate the two scattered polarization components. With HW2 at 00, only
the parallel signal will make it past the analyzer, with HW2 at 450, only the
perpendicular signal is transmitted. The system achieved an extinction ratio
$(I_{\perp}/I_{\parallel})$ between 1/500 and 1/100 depending on the
orientation of HW1. After the analyzer, the remaining Rayleigh component of
the signal is attenuated by two (Ondax SureBlock ultra-narrowband) notch
filters and the signal is directed into a one meter (Horiba FHR 1000)
spectrometer and detected by a (Horiba Synapse) CCD.

Figure S1: (a) Schematic of experimental setup. Laser light comes from top.
Scattered signal goes right for detection. (b) Diagram for the different
polarizations involved in the measurement.
$\overset{\rightarrow}{\boldsymbol{E}_{i}}$ is the incident linear
polarization. $\overset{\rightarrow}{\boldsymbol{E}_{s}}$ is the polarization
for the scattered signal and need not be linearly polarized. Parallel
measurement (HW2 at 00) probes the scattered field’s projection on
$\hat{\boldsymbol{e}}_{\parallel}$. The perpendicular measurement (HW2 at 450)
probes the scattered field’s projection on $\hat{\boldsymbol{e}}_{\perp}$.
Figure S2: NaCl Raman PO color map for measurement and TDEP simulation in 400
K. Top Row is for experimental measurement, bottom row for calculation. Left
column is for parallel configuration and right is for perpendicular. Figure
S3: NaCl Raman PO color map for measurement and TDEP simulation in 500 K. Top
Row is for experimental measurement, bottom row for calculation. Left column
is for parallel configuration and right is for perpendicular. Figure S4: NaCl
Raman PO color map for measurement and TDEP simulation in 600 K. Top Row is
for experimental measurement, bottom row for calculation. Left column is for
parallel configuration and right is for perpendicular. Figure S5: NaCl Raman
PO color map for measurement and TDEP simulation in 700 K. Top Row is for
experimental measurement, bottom row for calculation. Left column is for
parallel configuration and right is for perpendicular.

##  SII Measurement and analysis

Here We give the detailed description of all data acquisition processing and
present all PO maps for all measured temperatures.

To accumulate a Raman crystallography data-set, we rotate the polarization
orientation of the incident laser beam $(\overrightarrow{\boldsymbol{E}}_{i})$
in Fig. S1 (b) in $10^{0}$ steps covering $\theta=0-360^{0}$ in a plane
parallel to the NaCl (001) crystal face. At each $\theta$ step, a Raman
spectrum is collected for both HW2 configurations (see
$\hat{\boldsymbol{e}}_{\parallel}$ and $\hat{\boldsymbol{e}}_{\perp}$ in Fig.
S1).

The signal of 2nd order Raman is very weak, so long exposure times and
multiple accumulations are necessary. For each measured spectrum, that is, for
a given temperature, at a given incident polarization angle at given detection
configuration, approximately 6 accumulations of 60 seconds were collected.

Some residual baseline persisted in all measured spectra, resembling a
Lorentzian tail towards lower frequencies. This was cancelled in all
unpolarized figures via standard baseline reduction. All PO maps are given in
figs. S2, S3, S4 and S5. These are presented with no augmentation other than a
low-frequency cut-off, where the signal is governed by the notch filters and
not the actual scattering. For each temperature, the corresponding TDEP
calculation (see SVII) is also presented.

##  SIII Single Crystal Synthesis

The NaCl single crystal used in Raman measurements was home grown. Double-
deionized water (DDW) was heated slowly to around $40-50^{0}C$. High purity
$(\geq 99\%)$ Sodium Chloride powder was added and stirred into the DDW until
saturation. The saturated solution was poured into another wide flask with a
few fresh crystals for crystallization sites. The flask was left covered with
a filter paper (to avoid dust penetration and allow evaporation) in a chemical
hood for several days until the crystals reached a size of $\sim$1cm. Finally,
the crystals were taken out of the dried flask and placed on a clean wipe to
adsorb all residual water. A picture of the crystal used is given in figure
S6.

Figure S6: One of the NaCl single crystals synthesized for the experiment..

##  SIV Factor group analysis for 2nd order Raman

The full factor group analysis for the rock-salt structure has already been
worked out by Burnstein [1], following methods established by Birman [2] and
Elliot [3]. These ostensibly provide a complete tool-set to interpret 2nd
order Raman spectra by assigning each spectral feature to a known combination
of vibrational modes. As mentioned in the main text, this treatment may prove
useful for mode assignment of a known phononic dispersion relations to an
observed spectrum. However, even in the high symmetry case of the cubic
$O_{h}$ point group, selection rules for 2nd order Raman are too relaxed to
offer complete predictive power, even with reliable dispersion relations at
hand.

To demonstrate this, we perform the factor group analysis for for one high
symmetry point in the Birllouin zone, specifically, for the
$\mathrm{X^{\prime}_{5}(TO)},\>\mathrm{X^{\prime}_{5}(TA)}$ phonon pair in the
reciprocal space X-point in a rock-salt structure.

The basic selection rule governing Raman scattering is [4]:

$\Gamma_{f}\subseteq\Gamma_{15}\otimes\Gamma_{15}$ (S1)

with $\Gamma_{f}$ the representation of the excited phonon state, and
$\Gamma_{15}$ the representation for the momentum operator in BSW notation
(Raman scattering can be viewed as the product o two dipole transitions). The
crux of the matter lies in the fact that in 2nd order Raman this final phonon
state is generally composed of two off $\Gamma$-point phonons, with its
representation given by the Kronecker product of the two single-phonon states.
For example:

$\displaystyle\textrm{Quantum
state:}\>\>\left|\mathrm{X}_{2phonon}\right\rangle=\left|\mathrm{X^{\prime}_{5}(TO),\>X^{\prime}_{5}(TA)}\right\rangle\>\to\>$
$\displaystyle\textrm{Representation:}\>\>\>\mathrm{X}_{2phonon}=\mathrm{X^{\prime}_{5}(TO)\otimes
X^{\prime}_{5}(TA)=X_{1}\oplus X_{2}\oplus X_{3}\oplus X_{4}}$ (S2)

with all the irreducible representations in the direct sum being Raman active
according to Eq. (S1). This leads to the symmetry allowed Raman tensor form
[5]:

$R_{\mathrm{X^{\prime}}_{5}\times\mathrm{X^{\prime}}_{5}}=R_{\mathrm{X}_{1}}+R_{\mathrm{X}_{2}}+R_{\mathrm{X}_{3}}+R_{\mathrm{X}_{4}}=\\\
\begin{pmatrix}a&d&0\\\ d&b&0\\\ 0&0&c\end{pmatrix}$ (S3)

which is a very general Raman tensor indeed. Now, consider that this tensor is
responsible for the contribution to the scattered intensity inside some
frequency interval $I(\Omega)\mathrm{d}\Omega$ from a single phonon
combination. In practice, equation $(11)$ allows a multitude of combinations
to contribute to the same interval, such that accurately predicting the
observed intensity, let alone the PO behavior of the signal, becomes
intractable. The situation becomes even worse for lower symmetry structures,
where more phonon combinations will be Raman active and have more general
Raman tensors.

##  SV Numerical determination of interaction tensors

Here we briefly reiterate the TDEP procedure for consistency of notation and
to clarify the analogy between how dipole and polarizibility interactions are
determined and how interatomic force constants are determined. The starting
point is a set of displacements $\boldsymbol{u}$ and forces $\boldsymbol{f}$
from a set of $N_{c}$ supercells sampled from a canonical ensemble at
temperature $T$. Although omitted in the notation, all interaction tensors
depend on temperature as well as the volume (or more generally strain). The
interatomic force constants are determined as follows: consider a supercell
with $N_{a}$ atoms and forces and displacements given by the $3N_{a}\times 1$
vectors $\boldsymbol{u}$ and $\boldsymbol{f}$ (for clarity we will make note
of the dimensions of the matrices in each step):

$\displaystyle\underbrace{\boldsymbol{f}}_{3N_{a}\times 1}$
$\displaystyle=-\underbrace{\boldsymbol{\Phi}}_{3N_{a}\times
3N_{a}}\underbrace{\boldsymbol{u}}_{3N_{a}\times 1}$ (S4) this can
equivalently be written as
$\displaystyle\underbrace{\boldsymbol{f}}_{3N_{a}\times 1}$
$\displaystyle=-\underbrace{(\boldsymbol{I}\otimes\boldsymbol{u}^{T})}_{3N_{a}\times(3N_{a})^{2}}\underbrace{\boldsymbol{\Phi}_{v}}_{(3N_{a})^{2}\times
1}$ (S5) with the Kronecker product $\otimes$. To express the forces from
$N_{c}$ supercells the matrices are stacked on top of each other:
$\displaystyle\underbrace{\begin{pmatrix}\boldsymbol{f}_{1}\\\ \vdots\\\
\boldsymbol{f}_{N_{c}}\end{pmatrix}}_{3N_{a}N_{c}\times 1}$
$\displaystyle=-\underbrace{\begin{pmatrix}\boldsymbol{I}\otimes\boldsymbol{u}_{1}^{T}\\\
\vdots\\\ \boldsymbol{I}\otimes\boldsymbol{u}_{N_{c}}^{T}\\\
\end{pmatrix}}_{3N_{a}N_{c}\times(3N_{a})^{2}}\underbrace{\boldsymbol{\Phi}_{v}}_{(3N_{a})^{2}\times
1}$ (S6) I can without any loss of generality express the second order force
constants as a linear combination of $N_{x}$ irreducible components
$\boldsymbol{x}^{\Phi\textrm{II}}$, where $N_{x}\leq(3N_{a})^{2}$, via a
matrix $\boldsymbol{C}^{\Phi\textrm{II}}$:
$\displaystyle\underbrace{\begin{pmatrix}\boldsymbol{f}_{1}\\\ \vdots\\\
\boldsymbol{f}_{N_{c}}\end{pmatrix}}_{3N_{a}N_{c}\times 1}$
$\displaystyle=-\underbrace{\underbrace{\begin{pmatrix}\boldsymbol{I}\otimes\boldsymbol{u}_{1}^{T}\\\
\vdots\\\ \boldsymbol{I}\otimes\boldsymbol{u}_{N_{c}}^{T}\\\
\end{pmatrix}}_{3N_{a}N_{c}\times(3N_{a})^{2}}\underbrace{\boldsymbol{C}^{\Phi\textrm{II}}}_{(3N_{a})^{2}\times
N_{x}}}_{=\boldsymbol{A}^{\Phi\textrm{II}}\,,3N_{a}N_{c}\times
N_{x}}\underbrace{\boldsymbol{x}^{\Phi\textrm{II}}}_{N_{x}\times 1}$ (S7)

The set of symmetry relations and invariances determines $\boldsymbol{C}$ and
$N_{x}$, [6, 7, 8] and in general $N_{x}\ll(3N_{a})^{2}$ which considerably
simplifies the problem of numerically determining $\boldsymbol{x}$. In
practice, only the matrix $\boldsymbol{A}$ is determined and stored. The
nominally very large matrices that go into the construction of
$\boldsymbol{A}$ are quite sparse and construction presents negligible
computational cost. The second order force constants imply matrices that could
possibly be treated naively, but the generalization to higher order quickly
becomes intractable with dense matrix storage. The third order force
constants, for example, comes down to

$\underbrace{\begin{pmatrix}\boldsymbol{f}_{1}\\\ \vdots\\\
\boldsymbol{f}_{N_{c}}\end{pmatrix}}_{3N_{a}N_{c}\times
1}=-\underbrace{\underbrace{\begin{pmatrix}\boldsymbol{I}\otimes\boldsymbol{u}_{1}^{T}\otimes\boldsymbol{u}_{1}^{T}\\\
\vdots\\\
\boldsymbol{I}\otimes\boldsymbol{u}_{N_{c}}^{T}\otimes\boldsymbol{u}_{N_{c}}^{T}\\\
\end{pmatrix}}_{3N_{a}N_{c}\times(3N_{a})^{3}}\underbrace{\boldsymbol{C}^{\Phi\textrm{III}}}_{(3N_{a})^{3}\times
N_{x}}}_{=\boldsymbol{A}^{\Phi\textrm{III}}\,,3N_{a}N_{c}\times
N_{x}}\underbrace{\boldsymbol{x}^{\Phi\textrm{III}}}_{N_{x}\times 1}$ (S8)

where the contracted matrix $\boldsymbol{A}$ is many orders of magnitude
smaller than the matrices it is built from. The interatomic force constants
are determined in succession:

$\displaystyle\boldsymbol{A}^{\Phi\textrm{II}}\boldsymbol{x}^{\Phi\textrm{II}}$
$\displaystyle=\boldsymbol{f}$ (S9)
$\displaystyle\boldsymbol{A}^{\Phi\textrm{III}}\boldsymbol{x}^{\Phi\textrm{III}}$
$\displaystyle=\boldsymbol{f}-\boldsymbol{A}^{\Phi\textrm{II}}\boldsymbol{x}^{\Phi\textrm{II}}$
(S10)
$\displaystyle\boldsymbol{A}^{\Phi\textrm{IV}}\boldsymbol{x}^{\Phi\textrm{IV}}$
$\displaystyle=\boldsymbol{f}-\boldsymbol{A}^{\Phi\textrm{II}}\boldsymbol{x}^{\Phi\textrm{II}}-\boldsymbol{A}^{\Phi\textrm{III}}\boldsymbol{x}^{\Phi\textrm{III}}$
(S11)

Where $\boldsymbol{f}$ denotes the reference forces to be reproduced.
Equations (S9)–(S11) are overdetermined and solved as a least squares problem.
This ensures that the baseline harmonic part becomes as large as possible, and
that the higher order terms become smaller and smaller.

In a similar manner we search for a way to determine $M$ and $P$ (the
expansion of the dipole moment and polarizability with respect to position) to
all orders. The input is again a supercell, where in addition to forces we
also have determined $\partial^{2}U/\partial E^{2}$ as well as
$\partial^{2}U/\partial E\partial R$, that is the polarizability and Born
charges. In a manner analogous to the force constants we determine the
irreducible components of the interaction tensors and express the
polarizability for a supercell as:

$\begin{split}\underbrace{\partial_{E}\partial_{E}U(\boldsymbol{u})}_{9\times
1}&=P_{0}+\underbrace{\boldsymbol{I}\otimes\boldsymbol{u}^{T}}_{9\times
27N_{a}}\underbrace{P}_{27N_{a}\times
1}+\underbrace{\boldsymbol{I}\otimes\boldsymbol{u}^{T}\otimes\boldsymbol{u}^{T}}_{9\times
9(3N_{a})^{2}}\underbrace{P}_{9(3N_{a})^{2}\times 1}+\ldots\\\
&=P_{0}+\boldsymbol{A}^{P\textrm{I}}\boldsymbol{x}^{P\textrm{I}}+\boldsymbol{A}^{P\textrm{II}}\boldsymbol{x}^{P\textrm{II}}+\ldots\end{split}$
(S12)

The terms in the polarizability expansion are determined slightly differently
than the interatomic force constants. For a set of $N_{c}$ supercells I build
the set of differences in polarizability between supercell $i$ and $j$:

$\displaystyle\Delta^{0}P_{ij}$
$\displaystyle=\partial_{E}\partial_{E}U(\boldsymbol{u}_{i})-\partial_{E}\partial_{E}U(\boldsymbol{u}_{j})$
(S13) and solve for $\boldsymbol{x}^{P\textrm{I}}$ from
$\displaystyle\left(\boldsymbol{A}^{P\textrm{I}}_{i}-\boldsymbol{A}^{P\textrm{I}}_{j}\right)\boldsymbol{x}^{P\textrm{I}}$
$\displaystyle=\Delta^{0}P_{ij}$ (S14) where the set of all pairs of supercell
$ij$ are used. The fluctuations from the first order are removed to construct
$\displaystyle\Delta^{1}P_{ij}$
$\displaystyle=\left(\partial_{E}\partial_{E}U(\boldsymbol{u}_{i})-\boldsymbol{A}^{P\textrm{I}}_{i}\boldsymbol{x}^{P\textrm{I}}\right)-\left(\partial_{E}\partial_{E}U(\boldsymbol{u}_{j})-\boldsymbol{A}^{P\textrm{I}}_{j}\boldsymbol{x}^{P\textrm{I}}\right)$
(S15)
$\displaystyle\left(\boldsymbol{A}^{P\textrm{II}}_{i}-\boldsymbol{A}^{P\textrm{II}}_{j}\right)\boldsymbol{x}^{P\textrm{II}}$
$\displaystyle=\Delta^{1}P_{ij}$ (S16)

And so on for all orders considered. The baseline polarizability is then
determined via

$P_{0}+P^{\textrm{th}}=\left\langle\partial_{E}\partial_{E}U-\boldsymbol{A}^{P\textrm{I}}\boldsymbol{x}^{P\textrm{I}}-\boldsymbol{A}^{P\textrm{II}}\boldsymbol{x}^{P\textrm{II}}-\ldots\right\rangle$
(S17)

Where $P^{\textrm{th}}$ is the contribution to the susceptibility from the
second order terms in the polarizability expansion [9]:

$P^{\textrm{th}}=\sum_{\lambda}P_{\lambda\lambda}(2n_{\lambda}+1)$ (S18)

where the matrix elements are defined below. For a perfectly harmonic system
(or any harmonic sampling of phase space) the average of the first order term
disappears.

The interaction parameters for the dipole moment are determined from the set
of Born charges calculated in the supercell. We deliberately do not work with
polarization directly. The quantized nature of polarization makes finding the
polarization on the same branch for a large number of supercells tedious if
not outright problematic. The starting point is again the idea of expressing
the Born charges for a supercell in terms of the irreducible components:

$\begin{split}\underbrace{\partial_{E}\partial_{R}U}_{9N_{a}\times
1}(\boldsymbol{u})&=\underbrace{\boldsymbol{M}^{\textrm{I}}}_{9N_{a}\times
1}+\underbrace{\boldsymbol{I}\otimes\boldsymbol{u}^{T}}_{9N_{a}\times(9N_{a})(3N_{a})}\underbrace{M^{\textrm{II}}}_{(9N_{a})(3N_{a})\times
1}+\underbrace{\boldsymbol{I}\otimes\boldsymbol{u}^{T}\otimes\boldsymbol{u}^{T}}_{9N_{a}\times(9N_{a})(3N_{a})^{2}}\underbrace{M^{\textrm{III}}}_{(9N_{a})(3N_{a})^{2}\times
1}+\ldots\\\
&=\boldsymbol{A}^{\textrm{MI}}\boldsymbol{x}^{\textrm{MI}}+\boldsymbol{A}^{\textrm{MII}}\boldsymbol{x}^{\textrm{MII}}+\boldsymbol{A}^{\textrm{MIII}}\boldsymbol{x}^{\textrm{MIII}}+\end{split}$
(S19)

Where the terms are determined in succession:

$\displaystyle\boldsymbol{A}^{\textrm{MI}}\boldsymbol{x}^{\textrm{MI}}$
$\displaystyle=\partial_{E}\partial_{R}U$
$\displaystyle\boldsymbol{A}^{\textrm{MII}}\boldsymbol{x}^{\textrm{MII}}$
$\displaystyle=\partial_{E}\partial_{R}U-\boldsymbol{A}^{\textrm{MI}}\boldsymbol{x}^{\textrm{MI}}$
$\displaystyle\boldsymbol{A}^{\textrm{MIII}}\boldsymbol{x}^{\textrm{MIII}}$
$\displaystyle=\partial_{E}\partial_{R}U-\boldsymbol{A}^{\textrm{MI}}\boldsymbol{x}^{\textrm{MI}}-\boldsymbol{A}^{\textrm{MII}}\boldsymbol{x}^{\textrm{MII}}$

##  SVI Long-ranged electrostatics in polar materials

Long-ranged interactions are not captured appropriately by the interatomic
force constants and have to be treated with care. The formalism to deal with
this in quasiharmonic calculations is well established, see e.g. Gonze and Lee
[10]. For finite temperature calculations of interactions beyond pair
interaction the formalism requires a slight adjustment.

We will first discuss the problem in quite general terms. Our starting point
is a general long-ranged (polar) pair interaction $f(r)$ that we range-
separate in some way:

$f(r)=f^{\textrm{sr}}(r)+f^{\textrm{lr}}(r).$ (S20)

If we set the range-separation such that the short-ranged $f^{\textrm{sr}}$
part fits snugly in the simulation cell, that part of the interaction is taken
care of by the usual force constants and requires no additional consideration.
The long-ranged tail $f^{\textrm{lr}}$ is what remains to be treated. This is
done by calculating the Hessian of the long-ranged part:

$\Phi^{\textrm{lr}}_{ij\alpha\beta}=\left.\frac{\partial^{2}f^{\textrm{lr}}(r)}{\partial
u_{i}^{\alpha}\partial u_{j}^{\beta}}\right|_{r=r_{ij}}$ (S21)

after which the forces from the long-ranged interactions can be calculated

$f_{i\alpha}^{\textrm{lr}}=-\sum_{j\beta}\Phi^{\textrm{lr}}_{ij\alpha\beta}u_{j\beta}$
(S22)

these long-ranged forces are removed from the calculated forces, such that the
procedure for determining the interatomic force constants becomes:

$\displaystyle\boldsymbol{A}^{\Phi\textrm{II}}\boldsymbol{x}^{\Phi\textrm{II}}$
$\displaystyle=\boldsymbol{f}-\boldsymbol{f}^{\textrm{lr}}$ (S23)
$\displaystyle\boldsymbol{A}^{\Phi\textrm{III}}\boldsymbol{x}^{\Phi\textrm{III}}$
$\displaystyle=\boldsymbol{f}-\boldsymbol{f}^{\textrm{lr}}-\boldsymbol{A}^{\Phi\textrm{II}}\boldsymbol{x}^{\Phi\textrm{II}}$
(S24)
$\displaystyle\boldsymbol{A}^{\Phi\textrm{IV}}\boldsymbol{x}^{\Phi\textrm{IV}}$
$\displaystyle=\boldsymbol{f}-\boldsymbol{f}^{\textrm{lr}}-\boldsymbol{A}^{\Phi\textrm{II}}\boldsymbol{x}^{\Phi\textrm{II}}-\boldsymbol{A}^{\Phi\textrm{III}}\boldsymbol{x}^{\Phi\textrm{III}}$
(S25)

When determining the dynamical matrix, the long-ranged pair interactions are
added back again:

$\boldsymbol{\Phi}(\boldsymbol{q})=\boldsymbol{\Phi}^{\textrm{sr}}(\boldsymbol{q})+\boldsymbol{\Phi}^{\textrm{lr}}(\boldsymbol{q})$
(S26)

This procedure is beneficial in several ways. By subtracting forces (instead
of dynamical matrices as in Gonze and Lee [10]) numerical issues due to
aliasing (finite size) errors are greatly reduced. More importantly, it also
ensures there is no multiple counting: the long-ranged interactions are
subtracted from interactions of all orders, ensuring the short-sightedness of
higher order interactions. What we just described is the general procedure,
for practical calculations one needs to assume an explicit form for the long-
ranged interactions. We use the (now temperature-dependent) Born charges and
dielectric constant, and the range-separatation is done using the standard
Ewald technique. The resulting real-space force constants are then [10]:

$\displaystyle\Phi^{\textrm{dd}}|^{\alpha\beta}_{ij}$
$\displaystyle=\sum_{\gamma\delta}Z_{i}^{\alpha\gamma}Z_{j}^{\beta\delta}\widetilde{\Phi}^{\gamma\delta}_{ij}$
(S27) $\displaystyle\widetilde{\Phi}^{\alpha\beta}_{ij}$
$\displaystyle=\frac{1}{4\pi\epsilon_{0}}\frac{1}{\sqrt{\det\epsilon}}\left(\frac{\epsilon^{-1}_{\alpha\beta}}{|\Delta_{ij}|_{\epsilon}^{3}}-3\frac{\Delta_{ij}^{\alpha}\Delta_{ij}^{\beta}}{|\Delta_{ij}|_{\epsilon}^{5}}\right)$
(S28)

Here $\boldsymbol{\epsilon}$ is the dielectric tensor and $\Delta$ realspace
distances using the dielectric tensor as a metric:

$\displaystyle\mathbf{r}_{ij}$
$\displaystyle=\mathbf{R}_{j}+\mathbf{\tau}_{j}-\mathbf{\tau}_{i}$ (S29)
$\displaystyle\mathbf{\Delta}_{ij}$
$\displaystyle=\boldsymbol{\epsilon}^{-1}\mathbf{r}_{ij}$ (S30)
$\displaystyle|\Delta_{ij}|_{\epsilon}$
$\displaystyle=\sqrt{\mathbf{\Delta}_{ij}\cdot\mathbf{r}_{ij}}$ (S31)

This poses an issue when calculating the dynamical matrix: the interactions
die off as $~{}1/r^{3}$ which makes it necessary to extend the sum over
lattice vectors to infinity. This is remedied with the usual Ewald technique:

$\widetilde{\mathbf{\Phi}}_{ij}(\mathbf{q})=\widetilde{\mathbf{\Phi}}^{\textrm{r}}+\widetilde{\mathbf{\Phi}}^{\textrm{q}}+\widetilde{\mathbf{\Phi}}^{\textrm{c}}$
(S32)

Dividing the sum into a realspace part, a reciprocal part and a connecting
part. The realspace part is given by

$\displaystyle\widetilde{\mathbf{\Phi}}^{\textrm{r}}_{ij}$
$\displaystyle=-\frac{\Lambda^{3}}{4\pi\epsilon_{0}\sqrt{\det\epsilon}}\sum_{\mathbf{R}}\mathbf{H}(\Lambda\Delta_{ij},\Lambda|\Delta_{ij}|_{\epsilon})e^{i\mathbf{q}\cdot\mathbf{R}}$
(S33)
$\displaystyle\frac{\partial\widetilde{\mathbf{\Phi}}^{\textrm{r}}_{ij}}{\partial
q_{\alpha}}$
$\displaystyle=-\frac{\Lambda^{3}}{4\pi\epsilon_{0}\sqrt{\det\epsilon}}\sum_{\mathbf{R}}iR_{\alpha}\mathbf{H}(\Lambda\Delta_{ij},\Lambda|\Delta_{ij}|_{\epsilon})e^{i\mathbf{q}\cdot\mathbf{R}}$
(S34)

where

$H_{\alpha\beta}(\mathbf{x},y)=\frac{x_{\alpha}x_{\beta}}{y^{2}}\left[\frac{3\,\textrm{erfc}\,y}{y^{3}}+\frac{2e^{-y^{2}}}{\sqrt{\pi}}\left(\frac{3}{y^{2}}+2\right)\right]-\epsilon^{-1}_{\alpha\beta}\left[\frac{\textrm{erfc}\,y}{y^{3}}+\frac{2e^{-y}}{\sqrt{\pi}y^{2}}\right]\,.$
(S35)

In reciprocal space we have

$\displaystyle\widetilde{\Phi}_{\textrm{q}}|_{ij}$
$\displaystyle=\sum_{\mathbf{K}=\mathbf{G}+\mathbf{q}}\chi_{ij}(\mathbf{K},\Lambda)\left(\mathbf{K}\otimes\mathbf{K}\right)$
(S36) $\displaystyle\frac{\partial}{\partial
q_{x}}\widetilde{\Phi}_{\textrm{q}}|_{ij}$
$\displaystyle=\sum_{\mathbf{K}=\mathbf{G}+\mathbf{q}}\chi_{ij}(\mathbf{K},\Lambda)\left(\mathbf{K}\otimes\mathbf{K}\right)\left(i\tau^{x}_{ij}-\left[\sum_{\alpha}K_{\alpha}\epsilon_{\alpha
x}\right]\left[\frac{1}{\|\mathbf{K}\|_{\epsilon}}+\frac{1}{4\Lambda^{2}}\right]\right)+\chi_{ij}(\mathbf{K},\Lambda)\begin{pmatrix}2K_{x}&K_{y}&K_{z}\\\
K_{y}&0&0\\\ K_{z}&0&0\end{pmatrix}$ (S37)
$\displaystyle\frac{\partial}{\partial
q_{y}}\widetilde{\Phi}_{\textrm{q}}|_{ij}$
$\displaystyle=\sum_{\mathbf{K}=\mathbf{G}+\mathbf{q}}\chi_{ij}(\mathbf{K},\Lambda)\left(\mathbf{K}\otimes\mathbf{K}\right)\left(i\tau^{y}_{ij}-\left[\sum_{\alpha}K_{\alpha}\epsilon_{\alpha
y}\right]\left[\frac{1}{\|\mathbf{K}\|_{\epsilon}}+\frac{1}{4\Lambda^{2}}\right]\right)+\chi_{ij}(\mathbf{K},\Lambda)\begin{pmatrix}0&K_{x}&0\\\
K_{x}&2K_{y}&K_{z}\\\ 0&K_{z}&0\end{pmatrix}$ (S38)
$\displaystyle\frac{\partial}{\partial
q_{z}}\widetilde{\Phi}_{\textrm{q}}|_{ij}$
$\displaystyle=\sum_{\mathbf{K}=\mathbf{G}+\mathbf{q}}\chi_{ij}(\mathbf{K},\Lambda)\left(\mathbf{K}\otimes\mathbf{K}\right)\left(i\tau^{z}_{ij}-\left[\sum_{\alpha}K_{\alpha}\epsilon_{\alpha
z}\right]\left[\frac{1}{\|\mathbf{K}\|_{\epsilon}}+\frac{1}{4\Lambda^{2}}\right]\right)+\chi_{ij}(\mathbf{K},\Lambda)\begin{pmatrix}0&0&K_{x}\\\
0&0&K_{y}\\\ K_{x}&K_{y}&2K_{z}\end{pmatrix}$ (S39)

where

$\displaystyle\chi_{ij}(\mathbf{K},\Lambda)$
$\displaystyle=\frac{1}{\Omega\epsilon_{0}}\frac{\exp\left(i\mathbf{q}\cdot\mathbf{\tau}_{ij}\right)\exp\left(-\frac{\|\mathbf{K}\|_{\epsilon}}{4\Lambda^{2}}\right)}{\|\mathbf{K}\|_{\epsilon}}$
(S40) $\displaystyle\|\mathbf{K}\|_{\epsilon}$
$\displaystyle=\sum_{\alpha\beta}\epsilon_{\alpha\beta}K_{\alpha}K_{\beta}$
(S41)

and finally a connecting part given by

$\widetilde{\mathbf{\Phi}}_{\textrm{c}}|_{ij}=\delta_{ij}\frac{\Lambda^{3}}{3\epsilon_{0}\pi^{3/2}\sqrt{\det\boldsymbol{\epsilon}}}\widetilde{\boldsymbol{\epsilon}}$
(S42)

This expression constitutes the dipole-dipole interactions at all distances.
As explained in the beginning of the section we only want to subtract the
long-ranged components from the forces. This separation is realized by
choosing the Ewald parameter $\Lambda$ such that $\widetilde{\Phi}^{r}$
disappears at the supercell boundary, and keeping only the reciprocal space
part of the Ewald sum. So, to summarize, the procedure is as follows:

  * •

Determine Ewald parameter $\Lambda$ such that the short-ranged interactions
die off at the supercell boundary.

  * •

Calculate $\boldsymbol{\Phi}^{\textrm{lr}}(\boldsymbol{q}=0)$ for the
simulation supercell via Eq. (S36).

  * •

Use $\boldsymbol{\Phi}^{\textrm{lr}}(\boldsymbol{q}=0)$ to calculate the
forces $\boldsymbol{f}^{\textrm{lr}}$ originating from the long-ranged tails
of the dipole-dipole interactions.

  * •

Use equations (S23)–(S25) to determine interatomic force constants.

  * •

When determining the dynamical matrix at any $\boldsymbol{q}$, add
$\boldsymbol{\Phi}^{\textrm{lr}}(\boldsymbol{q})$ back again.

For reference we also specified the gradient of the dynamical matrix in case
the group velocities are needed. To conclude this section we note that there
is no principal difference between this approach and what is regularly used,
but in practice the presented approach is easier to implement, suffers from
less numerical artifacts and makes avoiding double-counting easier. The
double-counting we refer to arises because if we only correct the dynamical
matrices, some long-ranged polar interactions get picked up by the third and
fourth order force constants. The non-analytical behavior of Eq. (S28) is
already included in the the long-ranged reciprocal space part and independent
of Ewald coupling parameter.

##  SVII Determining the Raman spectra

In this section we outline the procedure to determine Raman spectra of an
anharmonic crystal. We will begin by sketching the procedure to obtain the
phonon-phonon correlation functions in the presence of three-phonon
anharmonicity via the equation of motion approach, since intermediate results
are needed to obtain the polarizability-polarizability correlation functions.

If we express the displacements in terms of creation and annihilation
operators we can identify the reciprocal space coefficients needed for the
perturbative expansion, repeated here for clarity and consistency of notation.
Start with the standard creation and annihilation operators for a phonon with
momentum $\boldsymbol{q}$ and polarization $s$ (when there is no loss of
clarity we use the compound index $\lambda$ to denote $\boldsymbol{q}s$):

$\displaystyle a_{\boldsymbol{q}s}=$
$\displaystyle\frac{1}{\sqrt{2N\hbar}}\sum_{i\alpha}\epsilon_{\boldsymbol{q}s}^{i\alpha}\left(\sqrt{m_{i}\omega_{\boldsymbol{q}s}}u_{i\alpha}-i\frac{p_{i\alpha}}{\sqrt{m_{i}\omega_{\boldsymbol{q}s}}}\right)e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}$
(S43) $\displaystyle a_{\boldsymbol{q}s}^{\dagger}=$
$\displaystyle\frac{1}{\sqrt{2N\hbar}}\sum_{i\alpha}\epsilon_{\boldsymbol{q}s}^{i\alpha\,\dagger}\left(\sqrt{m_{i}\omega_{\boldsymbol{q}s}}u_{i\alpha}+i\frac{p_{i\alpha}}{\sqrt{m_{i}\omega_{\boldsymbol{q}s}}}\right)e^{i\mathbf{q}\cdot\mathbf{r}_{i}}$
(S44) $\displaystyle
A_{\boldsymbol{q}s}=a_{\boldsymbol{q}s}+a_{\bar{\boldsymbol{q}}s}^{\dagger}=$
$\displaystyle\frac{2}{\sqrt{2N\hbar}}\sum_{i\alpha}\epsilon_{\boldsymbol{q}s}^{i\alpha}\sqrt{m_{i}\omega_{\boldsymbol{q}s}}u_{i\alpha}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}$
(S45) $\displaystyle
B_{\boldsymbol{q}s}=a_{\boldsymbol{q}s}-a_{\bar{\boldsymbol{q}}s}^{\dagger}=$
$\displaystyle-\frac{2i}{\sqrt{2N\hbar}}\sum_{i\alpha}\epsilon_{\boldsymbol{q}s}^{i\alpha}\frac{p_{i\alpha}}{\sqrt{m_{i}\omega_{\boldsymbol{q}s}}}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}$
(S46)

Here $\epsilon$ are phonon eigenvectors
($\boldsymbol{\epsilon}_{i}\cdot\boldsymbol{\epsilon}_{j}=\delta_{ij}$), $m$
atomic masses, $\omega$ frequencies. The vector $\boldsymbol{r}$ is a lattice
vector, and $u$ and $p$ are the position and momentum operators. We use the
notation $\bar{\boldsymbol{q}}=-\boldsymbol{q}$. I introduce scaled
eigenvectors (to simplify notation) via

$\displaystyle\upsilon_{\boldsymbol{q}s}^{i\alpha}=$
$\displaystyle\sqrt{\frac{\hbar}{2m_{i}\omega_{\boldsymbol{q}s}}}\epsilon_{\boldsymbol{q}s}^{i\alpha}$
(S47) In this notation the matrix elements pertaining to anharmonicity become
$\displaystyle\Phi_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}$
$\displaystyle=\frac{1}{3!}\sum_{ijk\alpha\beta\gamma}\upsilon_{\boldsymbol{q}s}^{i\alpha}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{j\beta}\upsilon_{\boldsymbol{q}^{\prime\prime}s^{\prime\prime}}^{k\gamma}\Phi_{ijk}^{\alpha\beta\gamma}e^{-i(\boldsymbol{q}\cdot\boldsymbol{r}_{i}+\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}_{j}+\boldsymbol{q}^{\prime\prime}\cdot\boldsymbol{r}_{k})}$
(S48)
$\displaystyle\Phi_{\lambda\lambda^{\prime}\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}$
$\displaystyle=\frac{1}{4!}\sum_{ijkl\alpha\beta\gamma\delta}\upsilon_{\boldsymbol{q}s}^{i\alpha}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{j\beta}\upsilon_{\boldsymbol{q}^{\prime\prime}s^{\prime\prime}}^{k\gamma}\upsilon_{\boldsymbol{q}^{\prime\prime\prime}s^{\prime\prime\prime}}^{l\delta}\Phi_{ijkl}^{\alpha\beta\gamma\delta}e^{-i(\boldsymbol{q}\cdot\boldsymbol{r}_{i}+\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}_{j}+\boldsymbol{q}^{\prime\prime}\cdot\boldsymbol{r}_{k}+\boldsymbol{q}^{\prime\prime\prime}\cdot\boldsymbol{r}_{l})}$
(S49) so that the anharmonic part of the Hamiltonian becomes $\displaystyle
H_{A}$
$\displaystyle=\Phi_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}A_{\lambda}A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}}+\Phi_{\lambda\lambda^{\prime}\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}A_{\lambda}A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}}A_{\lambda^{\prime\prime\prime}}+\ldots$
(S50)

In an analogous way we define the coefficients of the expansion of the
polarizability as

$\displaystyle P^{\mu\nu}_{\lambda}$
$\displaystyle=\sum_{i\alpha}\upsilon_{\boldsymbol{q}s}^{i\alpha}P_{i}^{\mu\nu,\alpha}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}_{i}}\Delta_{\boldsymbol{q}}=\sum_{i\alpha}\upsilon_{\boldsymbol{\Gamma}s}^{i\alpha}P_{i}^{\mu\nu,\alpha}=P^{\mu\nu}(s)$
(S51) $\displaystyle P^{\mu\nu}_{\lambda\lambda^{\prime}}$
$\displaystyle=\frac{1}{2!}\sum_{ij\alpha\beta}\upsilon_{\boldsymbol{q}s}^{i\alpha}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{j\beta}P_{ij}^{\mu\nu,\alpha\beta}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}_{i}}e^{-i\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}_{j}}\Delta_{\boldsymbol{q}\boldsymbol{q}^{\prime}}=$
(S52)
$\displaystyle=\frac{1}{2!}\sum_{ij\alpha\beta}\upsilon_{\boldsymbol{q}s}^{i\alpha}(\upsilon_{\boldsymbol{q}s^{\prime}}^{j\beta})^{\dagger}P_{ij}^{\mu\nu,\alpha\beta}e^{i\boldsymbol{q}\cdot(\boldsymbol{r}_{j}-\boldsymbol{r}_{i})}$
(S53)
$\displaystyle=\frac{1}{2!}\sum_{ij\alpha\beta}(\upsilon_{\boldsymbol{q}s}^{i\alpha})^{\dagger}\upsilon_{\boldsymbol{q}s^{\prime}}^{j\beta}P_{ij}^{\mu\nu,\alpha\beta}e^{-i\boldsymbol{q}\cdot(\boldsymbol{r}_{j}-\boldsymbol{r}_{i})}=P^{\mu\nu}(\boldsymbol{q},s,s^{\prime})=P^{\mu\nu}(\bar{\boldsymbol{q}},s,s^{\prime})^{\dagger}=P^{\mu\nu}(\bar{\boldsymbol{q}},s^{\prime},s)$
(S54) $\displaystyle
P^{\mu\nu}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}$
$\displaystyle=\frac{1}{3!}\sum_{ijk\alpha\beta\gamma}\upsilon_{\boldsymbol{q}s}^{i\alpha}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{j\beta}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{k\gamma}P_{ijk}^{\mu\nu,\alpha\beta\gamma}e^{-i(\boldsymbol{q}\cdot\boldsymbol{r}_{i}+\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}_{j}+\boldsymbol{q}^{\prime\prime}\cdot\boldsymbol{r}_{k})}$
(S55) and the dipole moment matrix elements become $\displaystyle
M^{\mu}_{\lambda}$
$\displaystyle=\sum_{i\alpha}\upsilon_{\boldsymbol{\Gamma}s}^{i\alpha}M_{i}^{\mu,\alpha}=M^{\mu}(s)$
(S56) $\displaystyle M^{\mu\nu}_{\lambda\lambda^{\prime}}$
$\displaystyle=\frac{1}{2!}\sum_{ij\alpha\beta}(\upsilon_{\boldsymbol{q}s}^{i\alpha})^{\dagger}\upsilon_{\boldsymbol{q}s^{\prime}}^{j\beta}M_{ij}^{\mu,\alpha\beta}e^{-i\boldsymbol{q}\cdot(\boldsymbol{r}_{j}-\boldsymbol{r}_{i})}=M^{\mu}(\boldsymbol{q},s,s^{\prime})=M^{\mu}(\bar{\boldsymbol{q}},s,s^{\prime})^{\dagger}=M^{\mu}(\bar{\boldsymbol{q}},s^{\prime},s)$
(S57) $\displaystyle M^{\mu}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}$
$\displaystyle=\frac{1}{3!}\sum_{ijk\alpha\beta\gamma}\upsilon_{\boldsymbol{q}s}^{i\alpha}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{j\beta}\upsilon_{\boldsymbol{q}^{\prime}s^{\prime}}^{k\gamma}M_{ijk}^{\mu,\alpha\beta\gamma}e^{-i(\boldsymbol{q}\cdot\boldsymbol{r}_{i}+\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}_{j}+\boldsymbol{q}^{\prime\prime}\cdot\boldsymbol{r}_{k})}$
(S58) Such that $\displaystyle P^{\mu\nu}$
$\displaystyle=P^{\mu\nu}_{0}+P^{\mu\nu}_{\lambda}A_{\lambda}+P^{\mu\nu}_{\lambda\lambda^{\prime}}A_{\lambda}A_{\lambda^{\prime}}+P^{\mu\nu}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}A_{\lambda}A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}}+\ldots$
(S59) $\displaystyle M^{\mu}$
$\displaystyle=M^{\mu}_{0}+M^{\mu}_{\lambda}A_{\lambda}+M^{\mu}_{\lambda\lambda^{\prime}}A_{\lambda}A_{\lambda^{\prime}}+M^{\mu}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}A_{\lambda}A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}}+\ldots$
(S60)

###  SVII.1 Phonon Green’s functions

We briefly sketch the derivation of the phonon self-energy in the presence of
third order anharmonicity since intermediate results will be used to express
the Raman and infrared spectra. Starting from the retarded Green’s function

$G^{X,Y}=-i\theta(t)\left\langle[X(t),Y^{\dagger}(0)]\right\rangle$ (S61)

we express the equations of motion of the phonon-phonon Green’s function via

$\displaystyle\partial_{t}A(t)$ $\displaystyle=-i[H,A]$ (S62a)
$\displaystyle\frac{d}{dt}G^{AA}_{\lambda\lambda^{\prime}}$
$\displaystyle=-i\omega_{\lambda}G^{BA}_{\lambda\lambda^{\prime}}$ (S62b)
$\displaystyle\frac{d}{dt}G^{BA}_{\lambda\lambda^{\prime}}$
$\displaystyle=i2\delta_{\lambda\lambda^{\prime}}\delta(t)-i\omega_{\lambda}G^{AA}_{\lambda\lambda^{\prime}}-i6\sum_{\mu_{1}\mu_{2}}\Phi_{\lambda\mu_{1}\mu_{2}}G^{AA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
(S62c) $\displaystyle i\frac{d}{dt}G^{AA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
$\displaystyle=\omega_{\mu_{1}}G^{BA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}+\omega_{\mu_{2}}G^{AB,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
(S62d) $\displaystyle i\frac{d}{dt}G^{AB,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
$\displaystyle=\omega_{\mu_{1}}G^{BB,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}+\omega_{\mu_{2}}G^{AA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}+\sum_{\nu_{1}\nu_{2}}6\Phi_{\mu_{2}\nu_{1}\nu_{2}}G^{AAA,A}_{\mu_{1}\nu_{1}\nu_{2},\lambda^{\prime}}$
(S62e) $\displaystyle i\frac{d}{dt}G^{BA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
$\displaystyle=\omega_{\mu_{1}}G^{AA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}+\omega_{\mu_{2}}G^{BB,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}+\sum_{\nu_{1}\nu_{2}}6\Phi_{\nu_{1}\nu_{2}\mu_{1}}G^{AAA,A}_{\nu_{1}\nu_{2}\mu_{2},\lambda^{\prime}}$
(S62f) $\displaystyle i\frac{d}{dt}G^{BB,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
$\displaystyle=\omega_{\mu_{1}}G^{AB,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}+\omega_{\mu_{2}}G^{BA,A}_{\mu_{1}\mu_{2},\lambda^{\prime}}$
(S62g) the four-point Green’s functions are decoupled [11, 12] via
$\left\langle abcd\right\rangle\approx\left\langle ab\right\rangle\left\langle
cd\right\rangle+\left\langle ac\right\rangle\left\langle
bd\right\rangle+\left\langle ad\right\rangle\left\langle bc\right\rangle$ to
yield
$\displaystyle\sum_{\nu_{1}\nu_{2}}6\Phi_{\mu_{2}\nu_{1}\nu_{2}}G^{AAA,A}_{\mu_{1}\nu_{1}\nu_{2},\lambda^{\prime}}$
$\displaystyle=12(2n_{\mu_{1}}+1)\sum_{\nu_{1}}\Phi_{\mu_{1}\mu_{2}\nu_{1}}G^{AA}_{\nu_{1}\lambda^{\prime}}$
(S62h)
$\displaystyle\sum_{\nu_{1}\nu_{2}}6\Phi_{\nu_{1}\nu_{2}\mu_{1}}G^{AAA,A}_{\nu_{1}\nu_{2}\mu_{2},\lambda^{\prime}}$
$\displaystyle=12(2n_{\mu_{2}}+1)\sum_{\nu_{1}}\Phi_{\mu_{1}\mu_{2}\nu_{1}}G^{AA}_{\nu_{1}\lambda^{\prime}}$
(S62i)

After Fourier transforming equations (S62a)–(S62i) to the frequency domain we
get a solvable set of equations that recover the usual definition of three-
phonon anharmonicity [6, 13, 14, 11]:

$G_{\lambda\lambda^{\prime}}(\Omega)^{-1}=G^{0}_{\lambda\lambda^{\prime}}(\Omega)^{-1}+\Sigma_{\lambda\lambda^{\prime}}$
(S63)

where

$\Sigma_{\lambda\lambda^{\prime}}(Z)=-18\sum_{\boldsymbol{q}_{1}\boldsymbol{q}_{2}s_{1}s_{2}}\Phi^{\lambda
s_{1}s_{2}}_{\boldsymbol{q}\bar{\boldsymbol{q}}_{1}\bar{\boldsymbol{q}}_{2}}\Phi^{\lambda^{\prime}s_{1}s_{2}}_{\bar{\boldsymbol{q}}\boldsymbol{q}_{1}\boldsymbol{q}_{2}}S(s_{1},s_{2},Z)$
(S64)

and

$S(s_{a},s_{b},Z)=(n_{a}+n_{b}+1)\left[\frac{1}{(\omega_{a}+\omega_{b}-Z)_{p}}-\frac{1}{(\omega_{a}+\omega_{b}+Z)_{p}}\right]+(n_{a}-n_{b})\left[\frac{1}{(\omega_{b}-\omega_{a}+Z)_{p}}-\frac{1}{(\omega_{b}-\omega_{a}-Z)_{p}}\right]$
(S65)

Naturally if one considers the four-phonon interactions you recover the
standard results, but the information we need right now is the intermediate
result for the three-point Green’s function expressed in terms of two-point
Green’s functions and three-phonon matrix elements:

$G_{\lambda\lambda^{\prime},\lambda^{\prime\prime}}=6\sum_{s}\Phi_{\bar{\lambda}\bar{\lambda}^{\prime}s}S(\lambda,\lambda^{\prime})G_{A_{s}A_{\lambda^{\prime\prime}}}$
(S66)

###  SVII.2 Raman spectra

If we start from the polarizability-polarizability correlation function and
insert the expansion of the polarizability in terms of phonon coordinates, Eq.
(S59), we get

$\begin{split}\left\langle
P(t)P(0)\right\rangle=i\overset{>}{G}(P,P)=&\boldsymbol{P}_{0}\otimes\boldsymbol{P}_{0}+\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}}i\overset{>}{G}(A_{\lambda},A_{\lambda^{\prime}})+\\\
+&\boldsymbol{P}_{\lambda\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}i\overset{>}{G}(A_{\lambda}A_{\lambda^{\prime}},A_{\lambda^{\prime\prime}}A_{\lambda^{\prime\prime\prime}})+\\\
+&\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}\lambda^{\prime\prime}}i\overset{>}{G}(A_{\lambda},A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}})+\boldsymbol{P}_{\lambda\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime}}i\overset{>}{G}(A_{\lambda}A_{\lambda^{\prime}},A_{\lambda^{\prime\prime}})+\\\
+&\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}i\overset{>}{G}(A_{\lambda},A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}}A_{\lambda^{\prime\prime\prime}})+\boldsymbol{P}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime\prime}}i\overset{>}{G}(A_{\lambda}A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}},A_{\lambda^{\prime\prime\prime}})\end{split}$
(S67)

where we have omitted terms of $A^{4}$ or higher. We will deal with the terms
in order. The constant term does not contribute to the Raman spectrum. The
first term that gives a contribution is

$I^{\textrm{I}}=\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}}\int
i\overset{>}{G}(A_{\lambda}(t),A_{\lambda^{\prime}}(0))e^{-i\Omega
t}dt=\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}}(n(\Omega)+1)J_{\lambda\lambda^{\prime}}(\Omega)$
(S68)

where $J$ is the phonon spectral function

$J(\Omega)=-\frac{1}{\pi}\Im\left\\{G(\Omega)\right\\}=\frac{i}{(n(\Omega,T)+1)}\overset{>}{G}(\Omega)$
(S69)

We identify Eq. (S68) as the first order Raman spectra. The next contribution
comes from

$I^{\textrm{II}}=\boldsymbol{P}_{\lambda\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}\int\overset{>}{G}(A_{\lambda}(t)A_{\lambda^{\prime}}(t),A_{\lambda^{\prime\prime}}(0)A_{\lambda^{\prime\prime\prime}}(0))e^{-i\Omega
t}dt$ (S70)

The first step is to decouple the four-point Green’s function [12]:

$\begin{split}\overset{>}{G}(A_{\lambda}(t)A_{\lambda^{\prime}}(t),A_{\lambda^{\prime\prime}}(0)A_{\lambda^{\prime\prime\prime}}(0))&\approx\overset{>}{G}(A_{\lambda}(t),A_{\lambda^{\prime\prime\prime}}(0))\overset{>}{G}(A_{\lambda^{\prime}}(t),A_{\lambda^{\prime\prime}}(0))+\\\
&+\overset{>}{G}(A_{\lambda}(t),A_{\lambda^{\prime\prime}}(0))\overset{>}{G}(A_{\lambda^{\prime}}(t),A_{\lambda^{\prime\prime\prime}}(0))+\\\
&+\underbrace{\overset{>}{G}(A_{\lambda}(t),A_{\lambda^{\prime}}(t))\overset{>}{G}(A_{\lambda^{\prime\prime}}(0),A_{\lambda^{\prime\prime\prime}}(0))}_{=\textrm{constant}}\end{split}$
(S71)

and note that the Fourier transform of a product becomes a convolution:

$\begin{split}I^{\textrm{II}}&=\boldsymbol{P}_{\lambda\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}\frac{1}{\pi^{2}}\bigg{[}\int(n(\Omega^{\prime})+1)J_{\lambda\lambda^{\prime\prime}}(\Omega^{\prime})(n(\Omega-\Omega^{\prime})+1)J_{\lambda^{\prime}\lambda^{\prime\prime\prime}}(\Omega-\Omega^{\prime})d\Omega^{\prime}+\\\
&+\int(n(\Omega^{\prime})+1)J_{\lambda\lambda^{\prime\prime\prime}}(\Omega^{\prime})(n(\Omega-\Omega^{\prime})+1)J_{\lambda^{\prime}\lambda^{\prime\prime}}(\Omega-\Omega^{\prime})d\Omega^{\prime}\bigg{]}=\\\
&=2\Re\left\\{\boldsymbol{P}_{\lambda\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}\right\\}\int(n(\Omega^{\prime})+1)J_{\lambda\lambda^{\prime\prime}}(\Omega^{\prime})(n(\Omega-\Omega^{\prime})+1)J_{\lambda^{\prime}\lambda^{\prime\prime\prime}}(\Omega-\Omega^{\prime})d\Omega^{\prime}\end{split}$
(S72)

where we relabelled the summation indices in the second term. This term is the
second-order Raman spectra, slightly more general than it is usually
presented. If we assume a diagonal self-energy (and consequently a diagonal
spectral function) and insert the non-interacting spectral function
$J^{0}_{\lambda}=\delta(\Omega-\omega_{\lambda})-\delta(\Omega+\omega_{\lambda})$
the integral becomes

$\begin{split}&\int
n(-\Omega^{\prime})\left[\delta(\Omega^{\prime}-\omega_{a})-\delta(\Omega^{\prime}+\omega_{a})\right]n(\Omega^{\prime}-\Omega)\left[\delta(\Omega-\Omega^{\prime}-\omega_{b})-\delta(\Omega-\Omega^{\prime}+\omega_{b})\right]d\Omega^{\prime}=\\\
=&n(-\omega_{a})n(\omega_{a}-\Omega)\left[\delta(\Omega-\omega_{a}-\omega_{b})-\delta(\Omega-\omega_{a}+\omega_{b})\right]-\\\
&n(\omega_{a})n(-\omega_{a}-\Omega)\left[\delta(\Omega+\omega_{a}-\omega_{b})-\delta(\Omega+\omega_{a}+\omega_{b})\right]=\\\
=&\left[n(\Omega)+1\right]\left[n(\omega_{a})+n(\omega_{b})+1\right]\left[\delta(\Omega+\omega_{a}+\omega_{b})-\delta(\Omega-\omega_{a}-\omega_{b})\right]+\\\
&+\left[n(\Omega)+1\right]\left[n(\omega_{a})-n(\omega_{b})\right]\left[\delta(\Omega-\omega_{a}+\omega_{b})-\delta(\Omega+\omega_{a}-\omega_{b})\right]=\\\
=&-\frac{n(\Omega)+1}{\pi}\Im\left\\{S(\omega_{a},\omega_{b},\Omega)\right\\}\end{split}$
(S73)

where we made use of $n(a+b)=(1+n(a)+n(b))/n(a)n(b)$ as well as
$n(-a)=-n(a)-1$. With this we recover the expression given by Cowley [9] as
the limit of weak anharmonicity. In this study, however, it was found
necessary to keep the interacting spectral function when calculating the Raman
spectra. It is intuitive: the second order Raman is proportional to the two-
phonon DOS, and if the phonons are broadened then the two-phonon DOS should
also be broadened.

Next we have a term given by

$I^{\textrm{III}}=\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}\lambda^{\prime\prime}}G^{>}(A_{\lambda},A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}})+\boldsymbol{P}_{\lambda\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime}}G^{>}(A_{\lambda}A_{\lambda^{\prime}},A_{\lambda^{\prime\prime}})$
(S74)

where we make use of the intermediate result for the three-point Green’s
functions in Eq. (S66) and arrive at

$I^{\textrm{III}}=-6(n(\Omega)+1)\sum_{\mu}\left(\boldsymbol{P}_{\lambda^{\prime\prime}\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda}\Phi_{\bar{\lambda}^{\prime\prime}\bar{\lambda}^{\prime}\mu}+\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}\lambda^{\prime\prime}}\Phi_{\bar{\lambda}^{\prime}\bar{\lambda}^{\prime\prime}\mu}\right)\Im\left\\{S_{\lambda^{\prime\prime}\lambda^{\prime}}\right\\}J_{\mu\lambda}(\Omega)$
(S75)

where, if we assume a diagonal self-energy, we recover previous results. [9]
It is worth noting is that – as far as we have discovered – the last two terms
are in practice several orders of magnitude smaller than the usual first and
second order terms. In term III we have a product between the two-phonon DOS
and the one-phonon spectral function, so to have any sizeable contribution
these need to peak simultaneously. The final term

$I^{\textrm{IV}}=\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}\lambda^{\prime\prime}\lambda^{\prime\prime\prime}}G(A_{\lambda},A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}}A_{\lambda^{\prime\prime\prime}})+\boldsymbol{P}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}\otimes\boldsymbol{P}_{\lambda^{\prime\prime\prime}}G(A_{\lambda}A_{\lambda^{\prime}}A_{\lambda^{\prime\prime}},A_{\lambda^{\prime\prime\prime}})$
(S76)

we approach with the same decoupling procedure as before:

$\begin{split}G^{>}(A_{\lambda}(t),A_{\lambda^{\prime}}(0),A_{\lambda^{\prime\prime}}(0)A_{\lambda^{\prime\prime\prime}}(0))&\approx
G^{>}(A_{\lambda}(t),A_{\lambda^{\prime\prime\prime}}(0))G^{>}(A_{\lambda^{\prime}}(0),A_{\lambda^{\prime\prime}}(0))+\\\
&+G^{>}(A_{\lambda}(t),A_{\lambda^{\prime\prime}}(0))G^{>}(A_{\lambda^{\prime}}(0),A_{\lambda^{\prime\prime\prime}}(0))+\\\
&+G^{>}(A_{\lambda}(t),A_{\lambda^{\prime}}(0))G^{>}(A_{\lambda^{\prime\prime}}(0),A_{\lambda^{\prime\prime\prime}}(0))=\\\
&=3G^{>}(A_{\lambda}(t),A_{\lambda^{\prime\prime\prime}}(0))\delta_{\lambda^{\prime}\lambda^{\prime\prime}}(2n_{\lambda^{\prime}}+1)\end{split}$
(S77)

where the factor 3 comes from relabelling indices. We also used the equal-time
thermal average of the phonon operators. This gives

$I^{\textrm{IV}}=3(n(\Omega)+1)\left(\boldsymbol{P}_{\lambda}\otimes\boldsymbol{P}_{\lambda^{\prime}\lambda^{\prime}\lambda^{\prime\prime}}+\boldsymbol{P}_{\lambda^{\prime\prime}\lambda^{\prime}\lambda^{\prime}}\otimes\boldsymbol{P}_{\lambda}\right)(2n_{\lambda^{\prime}}+1)J_{\lambda\lambda^{\prime\prime}}$
(S78)

What we derived in the section above holds for any operator that is expressed
as a constant multiplied with the phonon operator, i.e. a quantity that only
depends on the position of the atoms (it’s frequency dependence can be
ignored, or at least the cross-terms between the frequency dependence and the
position of atoms). The derivation of the infrared absorption spectrum is
identical, one only needs to replace the matrix elements. If the real part of
the susceptibility is of interest it is conveniently calculated via a Kramers-
Kronig transformation of the imaginary part.

In the main text we have defined compound matrix elements for the
polarizability to simplify the notation, they are given by:

$\displaystyle P^{(\textrm{I})}_{\mu\nu,\xi\rho}(s_{1})$
$\displaystyle=P^{\mu\nu}(s_{1})P^{\xi\rho}(s_{1})^{\dagger}$ (S79)
$\displaystyle P^{(\textrm{II})}_{\mu\nu,\xi\rho}(\boldsymbol{q},s_{1},s_{2})$
$\displaystyle=P^{\mu\nu}(\boldsymbol{q},s_{1},s_{2})P^{\xi\rho}(\boldsymbol{q},s_{1},s_{2})^{\dagger}$
(S80) $\displaystyle
P^{(\textrm{III})}_{\mu\nu,\xi\rho}(\boldsymbol{q},s_{1},s_{2},s_{3})$
$\displaystyle=\Phi_{\boldsymbol{q}\bar{\boldsymbol{q}}\boldsymbol{\Gamma}}^{s_{1}s_{2}s_{3}}P^{\mu\nu}(\boldsymbol{q},s_{1},s_{2})P^{\xi\rho}(s_{3})+\Phi_{\boldsymbol{q}\bar{\boldsymbol{q}}\boldsymbol{\Gamma}}^{s_{1}s_{2}s_{3}}P^{\xi\rho}(\boldsymbol{q},s_{1},s_{2})P^{\mu\nu}(s_{3})$
(S81) $\displaystyle
P^{(\textrm{IV})}_{\mu\nu,\xi\rho}(\boldsymbol{q},s_{1},s_{2})$
$\displaystyle=P^{\mu\nu}(\boldsymbol{\Gamma},\bar{\boldsymbol{q}},\boldsymbol{q},s_{1},s_{2},s_{2})P^{\xi\rho}(s_{1})+P^{\xi\rho}(\boldsymbol{\Gamma},\bar{\boldsymbol{q}},\boldsymbol{q},s_{1},s_{2},s_{2})P^{\mu\nu}(s_{1})$
(S82) In exactly the same way I can define the following dipole matrix
elements: $\displaystyle M^{(I)}_{\alpha\beta}(s_{1})$
$\displaystyle=M^{\alpha}_{s_{1}}(M^{\beta}_{s_{1}})^{\dagger}$ (S83)
$\displaystyle M^{(II)}_{\alpha\beta}(s_{1},s_{2},\boldsymbol{q})$
$\displaystyle=M^{\alpha}(^{s_{1}s_{2}}_{\boldsymbol{q}})M^{\beta}(^{s_{1}s_{2}}_{\boldsymbol{q}})^{\dagger}$
(S84) $\displaystyle
M^{(III)}_{\alpha\beta}(s_{1},s_{2},s_{3},\boldsymbol{q})$
$\displaystyle=\Phi_{\boldsymbol{q}\bar{\boldsymbol{q}}\boldsymbol{0}}^{s_{1}s_{2}s_{3}}M^{\beta}(^{s_{1}s_{2}}_{\boldsymbol{q}})M^{\alpha}_{s_{3}}+\Phi_{\boldsymbol{q}\bar{\boldsymbol{q}}\boldsymbol{0}}^{s_{1}s_{2}s_{3}}M^{\alpha}(^{s_{1}s_{2}}_{\boldsymbol{q}})M^{\beta}_{s_{3}}$
(S85) $\displaystyle M^{(IV)}_{\alpha\beta}(s_{1},s_{2},\boldsymbol{q})$
$\displaystyle=M^{\alpha}(^{s_{1}s_{2}s_{2}}_{\boldsymbol{0}\bar{\boldsymbol{q}}\boldsymbol{q}})M^{\beta}_{s_{1}}+M^{\beta}(^{s_{1}s_{2}s_{2}}_{\boldsymbol{0}\bar{\boldsymbol{q}}\boldsymbol{q}})M^{\alpha}_{s_{1}}$
(S86)

Here I have explicitly used momentum conservation in the matrix elements.
There is only a single q-vector that remains, which means that the
computational cost of determining the spectras (given interaction tensors) are
modest, at most a single sum over the Brillouin zone.

##  SVIII Details of calculations

###  SVIII.1 Self-consistency procedure

The interatomic force constants are determined self-consistently as explained
in Shulumba  _et al._ [15], but it is worth explaining the practical
procedure. The basic idea is that a set of force constants gives us normal
modes that can be thermally populated, generating thermally excited
structures. Thermally excited structures can, with the help of equation (S23)
yield a new set of force constants. The procedure is terminated when the input
and output force constants show no significant difference.

The self-consistency procedure is implemented via a geometric series:

  * •

Create an initial guess for the second order force constants (see Shulumba
_et al._ [15]) and generate 1 thermally excited supercell.

  * •

Generate new force constants and produce 2 thermally excited supercells.

  * •

Use both the one old and two new supercells to generate a new set of force
constants, and 4 new supercells.

  * •

Use 2+4 supercells to generate forceconstants and 8 supercells.

  * •

Use 4+8 supercells to generate forceconstants and 16 supercells.

  * •

Repeat until convergence.

This procedure has the benefit of introducing some mixing (since you use
supercells generated from different force constants) smoothing convergence. In
addition, the prescribed procedure leaves only 25% of calculations as waste. A
straight implementation (generate N supercells each iteration, repeat until
convergence) can waste a much larger fraction of calculations.

For this study all calculated spectra converged after 7 iterations, using
128+64 configurations. The calculations were done on a grid of volumes and
temperatures, and a few points were tested with one more iteration where no
appreciable difference in calculated quantities could be seen.

Figure S7: Phonon spectral functions for NaCl at different temperatures
calculated using the TDEP method. The white circles in panel (a) are
experimental neutron measurements from Raunio  _et al._ [16]. We show an
overall good agreement with experiment bar a slight uniform underestimation of
frequencies.

###  SVIII.2 Free energy minimization

To determine the equilibrium volume at any given temperature we did
calculations on a grid of volumes and temperatures. The grid is specified in
Table S1.

Temperature (K) | Lattice constants (Å)
---|---
0 | 5.233848 | 5.324865 | 5.419088 | 5.510051 | 5.603468 | 5.717206
150 | 5.255144 | 5.346531 | 5.441138 | 5.532472 | 5.626268 | 5.740469
300 | 5.282018 | 5.373873 | 5.468963 | 5.560764 | 5.655040 | 5.769825
500 | 5.326527 | 5.419155 | 5.515047 | 5.607621 | 5.702692 | 5.818444
700 | 5.380951 | 5.474526 | 5.571398 | 5.664918 | 5.760960 | 5.877895
1000 | 5.481180 | 5.576498 | 5.675174 | 5.770436 | 5.868267 | 5.987380
Table S1:  List of temperatures and lattice parameters used in the
simulations.

Once simulations are converged the free energy was determined via

$F=U_{0}+F_{\textrm{ph}}+\Delta F^{3\textrm{ph}}+\Delta F^{4\textrm{ph}}$
(S87)

where $F_{\textrm{ph}}$ is the usual phonon free energy and $U_{0}$ the
renormalized baseline free energy given by

$U_{0}=\left\langle U-\frac{1}{2!}\sum_{\begin{subarray}{c}ij\\\
\alpha\beta\end{subarray}}\overset{\textrm{lr}}{\Phi}_{ij}^{\alpha\beta}u_{i}^{\alpha}u_{j}^{\beta}-\frac{1}{2!}\sum_{\begin{subarray}{c}ij\\\
\alpha\beta\end{subarray}}\Phi_{ij}^{\alpha\beta}u_{i}^{\alpha}u_{j}^{\beta}-\frac{1}{3!}\sum_{\begin{subarray}{c}ijk\\\
\alpha\beta\gamma\end{subarray}}\Phi_{ijk}^{\alpha\beta\gamma}u_{i}^{\alpha}u_{j}^{\beta}u_{k}^{\gamma}-\frac{1}{4!}\sum_{\begin{subarray}{c}ijkl\\\
\alpha\beta\gamma\delta\end{subarray}}\Phi_{ijkl}^{\alpha\beta\gamma\delta}u_{i}^{\alpha}u_{j}^{\beta}u_{k}^{\gamma}u_{l}^{\delta}\right\rangle$
(S88)

Here it is important to note that we have to subtract the long-ranged polar
interactions to avoid double-counting them. The explicit anharmonic
contributions are given via [6, 13, 14]

$\Delta
F^{3\textrm{ph}}=-6\sum_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}\left|\Phi_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}\right|^{2}\left(\frac{3n_{\lambda}n_{\lambda^{\prime}}+3n_{\lambda}+1}{(\omega_{\lambda}+\omega_{\lambda^{\prime}}+\omega_{\lambda^{\prime\prime}})_{p}}+\frac{6n_{\lambda}n_{\lambda^{\prime\prime}}-3n_{\lambda}n_{\lambda^{\prime}}+3n_{\lambda^{\prime\prime}}}{(\omega_{\lambda}+\omega_{\lambda^{\prime}}-\omega_{\lambda^{\prime\prime}})_{p}}\right)+9\Phi_{\lambda\bar{\lambda}\lambda^{\prime\prime}}\Phi_{\lambda^{\prime}\bar{\lambda}^{\prime}\bar{\lambda}^{\prime\prime}}\frac{4n_{\lambda}(n_{\lambda^{\prime}}+1)+1}{(\omega_{\lambda^{\prime\prime}})_{p}}\,,$
(S89)

and

$\Delta
F^{4\textrm{ph}}=3\sum_{\lambda\lambda^{\prime}}\Phi_{\lambda\bar{\lambda}\lambda^{\prime}\bar{\lambda}^{\prime}}(2n_{\lambda}+1)(2n_{\lambda^{\prime}}+1)$
(S90)

The free energy was fitted to a Birch-Murnaghan equation of state, and
pressure was determined via $P=-dF/dV$.

###  SVIII.3 Parameters in calculations

All density functional theory calculations were done using the projector
augmented wave (PAW) [17] method as implemented in VASP [18, 19, 20, 21]. We
treated exchange-correlation within the AM05 approximation [22, 23]. We used a
288 atom supercell, and the Brillouim zone integrations used a $2\times
2\times 2$ Monkhorst-Pack mesh, [24] and a Gaussian smearing of 0.2 eV was
applied. Phonon self-energies and Raman spectra were integrated on a $51\times
51\times 51$ q-point mesh. The plane-wave energy cutoff was set to 350 eV.

The statistics needed to determine the polarizability and Born charges on a
large set of supercells used a slightly smaller simulation cell to offset the
increased computational cost. Here we used a 64 atom supercell with a tighter
$6\times 6\times 6$ k-point mesh.

## References

  * Burstein  _et al._ [1965] E. Burstein, F. A. Johnson, and R. Loudon, Selection rules for second-order infrared and raman processes in the rocksalt structure and interpretation of the Raman spectra of NaCl, KBr, and NaI, Physical Review 139, 10.1103/PhysRev.139.A1239 (1965).
  * Birman [1962] J. L. Birman, Space Group Selection Rules: Diamond and Zinc Blende, Physical Review 127, 1093 (1962).
  * Elliott and Loudon [1960] R. Elliott and R. Loudon, Group theory of scattering processes in crystals, Journal of Physics and Chemistry of Solids 15, 146 (1960).
  * Yu and Cardona [2010] P. Y. Yu and M. Cardona,  _Fundamentals of Semiconductors: Physics and Materials Properties_, Graduate Texts in Physics, Vol. 1 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010) pp. 1829–1841.
  * Aroyo  _et al._ [2006] M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, A. Kirov, and H. Wondratschek, Bilbao Crystallographic Server: I. Databases and crystallographic computing programs, Zeitschrift für Kristallographie - Crystalline Materials 221, 15 (2006).
  * Leibfried and Ludwig [1961] G. Leibfried and W. Ludwig, Theory of Anharmonic Effects in Crystals, Solid State Physics - Advances in Research and Applications 12, 275 (1961).
  * Maradudin and Vosko [1968] A. A. Maradudin and S. Vosko, Symmetry Properties of the Normal Vibrations of a Crystal, Reviews of Modern Physics 40, 1 (1968).
  * Born and Huang [1964] M. Born and K. Huang,  _Dynamical theory of crystal lattices_ (Oxford University Press, Oxford, 1964).
  * Cowley [1964] R. A. Cowley, The theory of Raman scattering from crystals, Proceedings of the Physical Society 84, 281 (1964).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Physical Review B 55, 10355 (1997).
  * Semwal and Sharma [1972] B. S. Semwal and P. K. Sharma, Thermal Conductivity of an Anharmonic Crystal, Physical Review B 5, 3909 (1972).
  * Zubarev [1960] D. N. Zubarev, Double-time Green functions in statistical physics, Soviet Physics Uspekhi 3, 320 (1960).
  * Cowley [1963] R. Cowley, The lattice dynamics of an anharmonic crystal, Advances in Physics 12, 421 (1963).
  * Wallace [1998] D. C. Wallace,  _Thermodynamics of crystals_, Dover Books on Physics (Dover Publications, Incorporated, 1998).
  * Shulumba  _et al._ [2017] N. Shulumba, O. Hellman, and A. J. Minnich, Intrinsic localized mode and low thermal conductivity of PbSe, Physical Review B 95, 014302 (2017).
  * Raunio  _et al._ [1969] G. Raunio, L. Almqvist, and R. Stedman, Phonon Dispersion Relations in NaCl, Physical Review 178, 1496 (1969).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996a).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59, 1758 (1999).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B 54, 11169 (1996b).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for open-shell transition metals, Physical Review B 48, 13115 (1993).
  * Armiento and Mattsson [2005] R. Armiento and A. Mattsson, Functional designed to include surface effects in self-consistent density functional theory, Physical Review B 72, 085108 (2005).
  * Mattsson and Armiento [2009] A. Mattsson and R. Armiento, Implementing and testing the AM05 spin density functional, Physical Review B 79, 155101 (2009).
  * Monkhorst and Pack [1976] H. Monkhorst and J. Pack, Special points for Brillouin-zone integrations, Physical Review B 13, 5188 (1976).


# Signatures of a Quantum Griffiths Phase close to
an Electronic Nematic Quantum Phase Transition

Pascal Reiss<EMAIL_ADDRESS>Clarendon Laboratory, Department
of Physics, University of Oxford, Oxford, UK  David Graf  National High
Magnetic Field Laboratory, Florida State University, Tallahassee, FL, USA
Amir A. Haghighirad  Clarendon Laboratory, Department of Physics, University
of Oxford, Oxford, UK  Institute for Quantum Materials and Technologies,
Karlsruhe Institute of Technology, Karlsruhe, Germany  Thomas Vojta
Department of Physics, Missouri University of Science and Technology, Rolla,
MO, USA  Amalia I. Coldea<EMAIL_ADDRESS>Clarendon
Laboratory, Department of Physics, University of Oxford, Oxford, UK

###### Abstract

In the vicinity of a quantum critical point, quenched disorder can lead to a
quantum Griffiths phase, accompanied by an exotic power-law scaling with a
continuously varying dynamical exponent that diverges in the zero-temperature
limit. Here, we investigate a nematic quantum critical point in the iron-based
superconductor FeSe0.89S0.11 using applied hydrostatic pressure. We report an
unusual crossing of the magnetoresistivity isotherms in the non-
superconducting normal state which features a continuously varying dynamical
exponent over a large temperature range. We interpret our results in terms of
a quantum Griffiths phase caused by nematic islands that result from the local
distribution of Se and S atoms. At low temperatures, the Griffiths phase is
masked by the emergence of a Fermi liquid phase due to a strong nematoelastic
coupling and a Lifshitz transition that changes the topology of the Fermi
surface.

11footnotetext: See Supplemental Material at [URL will be inserted by
publisher] for further experimental data and analyses

## Introduction

A central characteristic of finite- and zero-temperature phase transitions is
how the spatial and temporal correlation lengths evolve as the transition is
approached. For clean and continuous phase transitions, scaling theory
predicts power-law divergences of both correlation lengths as a function of
control parameter, with the critical exponents reflecting the universality
class. Moreover, the spatial and temporal correlation lengths are closely
related by the dynamics of the system Hertz (1975, 1976); Millis (1993);
Moriya (1985). In the presence of quenched disorder, this relation may be
lost. Quenched disorder is perfectly correlated in time, but can harbor a
spatially varying order parameter. In this situation, a smeared phase
transition can occur, where ordered islands form within a disordered bulk Kuo
_et al._ (2016); Cui and Fernandes (2018). Moreover, when order parameter
fluctuations within the islands are non-negligible, a Griffiths phase can
emerge which leads to continuously varying critical exponents as a function of
temperature and control parameter, fundamentally different to clean systems
Griffiths (1969); Fisher (1992, 1995); Del Maestro  _et al._ (2008); Vojta
(2003, 2006, 2010).

Figure 1: (a) Pressure-temperature phase diagram and (d) crystal structure of
FeSe0.89S0.11. The relative position of samples B and C under a pressure of
$p=4.7$ kbar are indicated by vertical lines. (b,c) The spatial distribution
of small S atoms induces locally varying critical pressures $p_{c}(r)$ and
random local strains. Close to the nematic quantum phase transition (QPT),
this leads to the formation of nematic islands. The scale shows the
experimental mean free path length $\lambda$ Reiss  _et al._ (2020).

Experimentally, quantum Griffiths phases have been identified in ferromagnetic
nickel vanadium alloys Ubaid-Kassis  _et al._ (2010); Wang  _et al._ (2017) as
well as a range of superconducting thin films Seidler  _et al._ (1992);
Schneider  _et al._ (2012); Shi  _et al._ (2014); Xing  _et al._ (2015); Saito
_et al._ (2018); Lewellyn  _et al._ (2019); Liu  _et al._ (2019); Lin  _et
al._ (2015). In the latter systems, sharp crossings of the magnetoresistivity
isotherms emerged as a distinctive experimental signature. A scaling analysis
revealed a temperature-dependent critical exponent $z\nu$ that diverges in the
low-temperature limit. This is the hallmark of a quantum Griffiths phase
($\nu$ is the critical correlation length exponent, and $z$ the dynamical
exponent) Fisher (1992, 1995); Del Maestro  _et al._ (2008); Vojta (2003,
2006, 2010).

In this Letter, we report the magnetoresistivity of the quasi-2D bulk
superconductor FeSe0.89S0.11 when tuned to the vicinity of its nematic quantum
critical point (QCP) using a hydrostatic pressure of $4.7$ kbar (Fig. 1(a))
Reiss  _et al._ (2020). Here, the magnetoresistivity isotherms show a
remarkably sharp crossing at about $30$ T over nearly two decades in
temperature up to $30$ K. Scaling of the magnetoresistivity yields a critical
exponent $z\nu$ that diverges at low temperatures, in agreement with the
quantum Griffiths scenario. We argue that the Griffiths phase is induced by
the local distribution of isoelectronic Se and S atoms that promote the
formation of nematic islands in the vicinity of the nematic QCP, as shown in
Fig. 1(b) and (c). Below a crossover temperature $T\approx 10$ K, the quantum
Griffiths phase and the associated QCP appear to be masked by an emergent non-
zero energy scale which coincides with the re-entrance of Fermi liquid
behavior attributed to a strong nematoelastic coupling, as well as a
topological Lifshitz transition of the Fermi surface.

## Methods

Single crystal of FeSe1-xSx with $x=0.11$ sulfur substitution were grown using
the KCl/AlCl3 chemical vapour transport method as described elsewhere Böhmer
_et al._ (2016). High-pressure, high-field measurements for samples B and C
were carried out in the 45 T hybrid DC facility in Tallahassee. We used Daphne
Oil 7575 as pressure medium which ensures hydrostatic conditions for much
higher pressures than reported here, and we used the Ruby fluorescence shifts
below 4 K to determine the pressure. Low-field measurements up to 13.5 T were
carried out on sample A in a QuantumDesign PPMS in Oxford. Here, Daphne Oil
7373 was used, and the pressure was determined by the superconducting
transition of tin. Samples were aligned with the magnetic field parallel to
the crystallographic $c$ axis to avoid breaking an in-plane symmetry.
Transport measurements were performed using a standard 4 or 5 contact setup,
using the AC LockIn technique with a low frequency $f\approx 20$ Hz, and a low
excitation current $I_{p}=1$ mA within the $(ab)$ plane.

## Results

Figure 2:  (a,b) The isothermal magnetoresistivity of samples B and C at a
pressure of $p=4.7$ kbar cross at $\mu_{0}H^{\star}\approx 28$ T,
$\rho^{\star}\approx 33\,\mu\Omega$ cm. Up and down sweeps show no hysteresis.
The pressure difference in brackets represents the distance to the critical
pressure (see Fig. 1(a) and in the SM Note (1)). (c) The same data as in panel
(a), but as a function of temperature in fixed field. (d) The actual and
extrapolated zero-field resistivities, $\rho(T)$ (solid line) and $\rho_{H\to
0}$ (points, see the SM Note (1)). Error bars are smaller than the symbol
size. The red dashed line is a fit to Fermi liquid behavior ($\rho_{0}\approx
4.9\,\mu\Omega$ cm, $A\approx 0.036\,\mu\Omega\,\text{cm}/\text{K}^{2}$). (e)
The first derivative reveals the onset of superconductivity, indicated by
arrows. Large quantum oscillations (QO) can be seen for $H>H^{\star}$. (f) The
extrapolation of the superconducting (SC) onset coincides with $H^{\star}$
only at $T=0$. All reported data are measured at $p=4.7$ kbar.

Figure 2(a) and (b) show the temperature dependence of the magnetoresistivity
of two different single crystals B and C of FeSe0.89S0.11 under a hydrostatic
pressure of $p=4.7$ kbar, which is in the immediate vicinity of their nematic
QCPs ($p_{c}=4.8$ kbar for sample B and $5.5$ kbar for sample C, respectively
as shown in Fig. 1(a) and in the Supplemental Material (SM) Note (1); Reiss
_et al._ (2020); Xiang  _et al._ (2017)). All magnetoresistivity isotherms
cross around a similar magnetic field, $\mu_{0}H^{\star}\approx 28.6$ T for
sample B and $28.0$ T for sample C, with similar resistivities
$\rho^{\star}\approx 32\,\mu\Omega$ cm and $34\,\mu\Omega$ cm, respectively.
This crossing occurs over nearly two decades in temperature
$0.3\,\text{K}\lesssim T\lesssim 30\,\text{K}$ and its significance becomes
evident in the resistivity plots as a function of temperature in constant
field, shown in Fig. 2(c). For $H<H^{\star}$, the resistivity follows a
metallic-like behavior with $\partial\rho/\partial T>0$ before the sample
becomes superconducting below $T_{c}^{\text{on}}\approx 10$ K (Fig. 2(d)).
Equivalently, the onset magnetic field, $H_{c2}^{\text{on}}$, between the
superconducting and normal phases can be identified in magnetic fields smaller
than $H^{\star}$ in Fig. 2(e) whose zero-temperature extrapolation coincides
with $H^{\star}$ (Fig. 2(f)). Thus, the magnetoresistivity crossing occurs
strictly within the non-superconducting normal phase for all finite
temperatures ($H^{\star}>H_{c2}^{\text{on}}(T)$), implying that this behavior
describes the normal phase in the vicinity of the nematic QCP.

In high magnetic fields above $H^{\star}$, the resistivity shows insulating-
like behavior ($\partial\rho/\partial T<0$), before it saturates below
$T\approx 2$ K (Fig. 2(c)), similar to previous reports Bristow  _et al._
(2020a). Despite this insulating-like behavior, the large magnetoresistivity
is a feature of the metallic, compensated multi-band system FeSe1-xSx Watson
_et al._ (2015); Coldea  _et al._ (2019); Bristow  _et al._ (2020a). Quantum
oscillations are visible for temperatures below $\approx 5$ K (Fig. 2(e)),
demonstrating the existence of a Fermi surface and highlighting the high
quality of the samples Reiss  _et al._ (2020); Shoenberg (1984). A two-band
analysis of the magnetoresistivity allows us to extrapolate the zero-field
resistivity from high magnetic fields Note (1), which indicate Fermi liquid
behavior below a crossover temperature $T_{\text{FL}}\approx 10$ K, shown in
Fig. 2(d) Reiss  _et al._ (2020); Bristow  _et al._ (2020a). The orbitally
averaged effective masses from quantum oscillations show non-divergent
electronic correlations in the vicinity of the nematic QCP, as discussed in
detail in Ref. 24, likely due to a coupling between the nematic order
parameter and the lattice Paul and Garst (2017); Wang and Berg (2019); de
Carvalho and Fernandes (2019); Vieira  _et al._ (2019).

Next, we use a prototypical power-law scaling ansatz to describe the
magnetoresistivity of FeSe0.89S0.11, previously applied in thin-film
materials, including dirty films of FeSe Seidler  _et al._ (1992); Schneider
_et al._ (2012); Shi  _et al._ (2014); Xing  _et al._ (2015); Saito  _et al._
(2018); Lewellyn  _et al._ (2019); Liu  _et al._ (2019); Lin  _et al._ (2015).
In $d$ dimensions, the scaling is given by

$\displaystyle\rho(H,T)/\rho^{\star}=T^{(2-d)/z}f\left(\mu_{0}|H-H^{\star}|/T^{1/z\nu}\right)$
(1)

with $f(0)=1$ and the critical exponent $z\nu$ Fisher (1990). Clearly, a
crossing of the magnetoresistivity isotherms at a finite $\rho^{\star}$ is
only possible for a two-dimensional system. Indeed, FeSe1-xSx have strongly
two-dimensional electronic and superconducting properties Terashima  _et al._
(2014); Reiss  _et al._ (2020); Coldea  _et al._ (2019); Bristow  _et al._
(2020b); Farrar  _et al._ (2020); Coldea (2020).

For a clean QCP, $z\nu$ is a constant given by the appropriate universality
class, and as described in the SM Note (1), this should lead to a straight
line in Fig. 3(a). This is evidently not the case here where $z\nu$ varies as
a function of temperature, shown in Fig. 3(b). Using this extracted $z\nu(T)$,
all magnetoresistivity data collapse onto a single curve for both samples,
reflecting the form of the scaling function $f$, shown in Fig. 3(c).
Deviations for this scaling only occur for the superconducting transition at
lowest fields and temperatures, and at the highest temperatures and fields.
These deviations indicate the limits of the scaling relation, as shown in the
SM Note (1).

Figure 3: (a) Log-log plot of the first derivative at the crossing field
$H^{\star}$ (refer to Fig. 2(e)). Error bars are smaller than the symbol size.
The slope of the dashed line corresponds to $1/z\nu(T)$. (b) Temperature
dependence of $z\nu$ extracted from panel A (crosses), and from the piece-wise
extraction shown in the SM (dots) Note (1). Error bars indicate a 1$\sigma$
confidence interval. (c) Scaled magnetotransport data using $z\nu\sim
T^{-1.5}$ (sample B) and $z\nu\sim T^{-0.95}$ (sample C). The superconducting
transitions (SC) deviate from this scaling form. (d) The low-field, low-
temperature resistivity in the mixed state follows a power-law form. All
reported data are measured at $p=4.7$ kbar. Figure 4: Evolution of
superconductivity for sample A. (a)-(c) Magnetoresistivity and (d)-(f) the
corresponding first derivatives showing the development of the superconducting
transition in magnetic fields for pressures across the nematic quantum phase
transition. (g) Three-dimensional $H$-$p$-$T$ superconducting phase diagram.
The inset shows the extrapolated critical fields, $H_{c2}^{\text{off}}$,
$H_{c2}^{\text{on}}$ and $H_{c2}^{\text{G}}$, defined in panels (d)-(f), in
the zero-temperature limit. The relative positions of samples B and C under a
pressure of 4.7 kbar in the phase diagram is indicated by dashed and the
dotted lines, respectively (see also Figs. 2 and 3). Error bars indicate a
1$\sigma$ confidence interval.

This scaling analysis reveals an interesting and unexpected feature. While the
zero-temperature divergence of the effective critical exponent $z\nu(T)$ is a
key signature of quantum Griffiths phases, the observed power-law form
$z\nu(T)\sim T^{\alpha}$ (with non-universal exponents $\alpha\approx-1.5$ for
sample B and $\approx-1.0$ for sample C) is much stronger than the logarithmic
(‘activated’) dependence predicted within the infinite-randomness criticality
scenario Vojta (2010); Lewellyn  _et al._ (2019); Note (1), as shown in Fig.
3(b). In fact, a power-law divergence of $z\nu(T)$ is incompatible with the
notion of a normal QCP because the temperature term $T^{1/z\nu(T)}$ in Eq. 1
remains finite for $T\to 0$, which implies the persistence of a non-zero
energy scale at lowest temperatures. Interestingly, we find that $z\nu$
deviates from the activated behavior dependence below $T\approx 5$–$10$ K
which coincides with a re-entrance of Fermi liquid behavior, Fig. 2(b). This
suggests a suppression of Griffiths fluctuations and/or dimensional changes,
e.g. topological changes of the Fermi surface Reiss  _et al._ (2020); Bristow
_et al._ (2020c); Coldea  _et al._ (2019); Coldea (2020).

We now focus on the nature of the underlying phases separated by $H^{\star}$.
Figure 3(d) shows that for fields smaller than $H^{\star}$, the resistivity
within the superconducting mixed state follows a power-law form $\rho\propto
T^{\alpha(H)}$ over almost one decade in temperature. We attribute this power-
law form to a disordered vortex-liquid phase that freezes into a vortex-glass
in the zero-temperature limit, as found in underdoped cuprates Shi  _et al._
(2014). Crossing over into the high-field regime above $H^{\star}$ where
quantum oscillations are present, the resistivity reflects the behavior of a
metallic phase with a partial charge-carrier localization, as discussed in the
SM Note (1).

To elucidate the origin and extend of the low-field disordered vortex phase,
we investigate the pressure dependence of the superconducting to normal
transitions in magnetic fields on sample A ($p_{c}\approx 5.8$ kbar Reiss  _et
al._ (2020)). Figure 4(a)-(f) shows the magnetoresistivity and its derivative
up to $13.5$ T inside the nematic phase ($0.9$ kbar), close to the nematic
quantum phase transition ($5.5$ kbar) and within the tetragonal phase ($8.8$
kbar). In the nematic and tetragonal phases, the normal-to-superconducting
transition widths are nearly temperature and field independent. In contrast, a
visible broadening of the transition is found close to $p_{c}$, but only for
high fields/low temperatures, thus coinciding with the vortex-liquid phase in
sample B. To quantify this additional broadening, we extract the
superconducting offset and onset critical fields, $H_{c2}^{\text{off}}$ and
$H_{c2}^{\text{on}}$, as shown in Fig. 4(d)-(f). Furthermore, we define a
critical magnetic field $H_{c2}^{\text{G}}$, where the magnetoresistivity
derivative has an additional shoulder before it returns to its high-
temperature normal state background, which is observable only in the vicinity
of $p_{c}$ (Fig. 4(e)). Figure 4(g) summarizes all extracted critical fields
and their zero-temperature extrapolations, see also the SM Note (1).
Interestingly, the zero-temperature superconducting transition width peaks at
the nematic quantum phase transition, doubling the extent of the
superconducting mixed state. The width of the $H_{c2}^{\text{G}}$(0) peak in
pressure is estimated to be $\sigma_{p}\approx 0.7(2)$ kbar, which agrees well
with an estimate for the pressure range of the quantum Griffiths phase, as
discussed below. Figure 4(g) also shows that the zero-field superconducting
transition does not display any similar broadening. This demonstrates that the
peak in $H_{c2}^{\text{G}}$ is a low-temperature and high-field effect, ruling
out effects of possible pressure inhomogeneities, as discussed in the SM Note
(1).

## Discussion

Quantum Griffiths phases were previously detected in inhomogeneous
superconductor-to-insulator transitions in thin films, including FeSe Seidler
_et al._ (1992); Schneider  _et al._ (2012); Shi  _et al._ (2014); Xing  _et
al._ (2015); Saito  _et al._ (2018); Lewellyn  _et al._ (2019); Liu  _et al._
(2019); Lin  _et al._ (2015). Here, in bulk FeSe0.89S0.11, the situation is
very different. The scaling relation only describes the normal state
resistivity and holds for magnetic fields up to $45$ T and temperatures up to
$30$ K, vastly exceeding the bulk superconducting phase. We therefore propose
that the quantum Griffiths phase in FeSe0.89S0.11 emerges from the suppression
of the nematic phase with pressure Reiss  _et al._ (2020); Xiang  _et al._
(2017); Matsuura  _et al._ (2017) and the formation of rare nematic islands in
a tetragonal matrix due to the random distribution of sulfur atoms (Fig. 1),
as suggested before (see the SM to Ref. 5). To demonstrate how this can lead
to a quantum Griffiths phase, we sample a random distribution of $11$ % S
atoms over a square lattice, and average the effective sulfur density $x(r)$
over the experimental quasiparticle mean-free path length $\lambda\approx 350$
Å Reiss  _et al._ (2020), as shown in Fig. 1(b) and (c). The intrinsic local
variation $\Delta x(r)\approx 0.4\%$ (std. dev.) gives rise to regions with
higher (lower) S content which have a locally lower (higher) critical pressure
$p_{c}(r)$. This is the prototypical case of random-mass disorder that smears
the quantum phase transition over a region $\Delta p_{c}(r)$. By comparing the
reported nematic transition temperatures from pressure and isoelectronic
substitution studies, we estimate $\Delta x(r)\propto\Delta p_{c}(r)\approx
0.4$ kbar Coldea  _et al._ (2019); Bristow  _et al._ (2020a); Reiss  _et al._
(2017); Xiang  _et al._ (2017). This estimate is similar to the observed
pressure range of a broadened superconducting transition, $\sigma_{p}=0.7(2)$
kbar. This suggests that the peak in $H_{c2}^{\text{G}}(0)$ occurs either due
to enhanced superconducting fluctuations within the nematic islands, and/or
superconducting nematic islands below the percolation threshold, which get
suppressed at $H^{\star}$. These effects could also provide a favorable
environment for the observed inhomogeneous superconducting vortex phase in the
vicinity of the nematic QCP. Finally, we note that the spatial arrangement of
the S atoms locally breaks the $C_{4}$ symmetry of the lattice and thus
introduces random-field effects. In the two-dimensional regime, they may limit
the size of the nematic domains, but for weak disorder, the corresponding
breakup length is exponentially large Seppälä and Alava (2001).

It is rather surprising that in such a clean system signatures of a quantum
Griffiths phase are detected in the vicinity of the electronic nematic quantum
phase transition. It is clear that there are a further ingredients at play
here as the nematoelastic coupling sets in and quenches the two-dimensional
quantum critical nematic fluctuations, below a cross-over temperature
$T_{\text{FL}}\approx 10$ K. As a result, Fermi liquid behavior with finite
electronic correlations is restored Reiss  _et al._ (2020); Bristow  _et al._
(2020c); Coldea  _et al._ (2019); Coldea (2020), and the quantum Griffiths
phase is cut off, leading to the overly strong divergence of $z\nu$.
Additionally, a band with likely 3D character is formed due to a Lifshitz
transition of the Fermi surface in the proximity of the nematic QCP Coldea
_et al._ (2019); Reiss  _et al._ (2020) which may change the effective
dimensionality of the system at low temperatures.

The observation of a quantum Griffiths phase in an iron-based superconductor
has a number of important implications and provides new insights into the
nature of nematic quantum phase transitions. Our study provides evidence that
the quantum Griffiths phase significantly affects the mixed state of the
superconducting phase. Moreover, it may be present in other families of iron-
based superconductors in which nematic and tetragonal phases may form over
limited compositional ranges around QCPs Böhmer  _et al._ (2015); Hosoi  _et
al._ (2016); Bristow  _et al._ (2020c); Coldea  _et al._ (2019). Finally, we
uncover a number of new experimental signatures in electronic transport
measurements, most notably the reported power-law behavior of $z\nu(T)$, that
could provide new insights into the dynamics of (quenched) nematic order
parameter fluctuations. We hope that our results will guide further
theoretical and experimental research in understanding nematic quantum
Griffiths phases.

## Acknowledgements

We thank R. Fernandes for insightful discussions. This work was mainly
supported by the EPSRC (EP/I004475/1, EP/I017836/1). P.R. and A.A.H.
acknowledge financial support of the Oxford Quantum Materials Platform Grant
(EP/M020517/1). A portion of this work was performed at the National High
Magnetic Field Laboratory, supported by NSF Cooperative Agreement DMR-1157490
and the State of Florida. Work in Missouri was supported by the NSF
(DMR-1828489). We acknowledge financial support of Oxford University John Fell
Fund. A.I.C. acknowledges an EPSRC Career Acceleration Fellowship
(EP/I004475/1) and Oxford Centre for Applied Superconductivity for financial
support.

## References

  * Hertz (1975) John A Hertz, “Zero- and low-temperature phase transitions and the renormalization group,” AIP Conference Proceedings 24, 298–299 (1975).
  * Hertz (1976) John A Hertz, “Quantum critical phenomena,” Physical Review B 14, 1165–1184 (1976).
  * Millis (1993) A J Millis, “Effect of a nonzero temperature on quantum critical points in itinerant fermion systems,” Physical Review B 48, 7183–7196 (1993).
  * Moriya (1985) Tôru Moriya,  _Spin Fluctuations in Itinerant Electron Magnetism_ , Springer Series in Solid-State Sciences 56 (Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1985).
  * Kuo  _et al._ (2016) H.-H. Kuo, J.-H. Chu, J. C. Palmstrom, S. A. Kivelson, and I. R. Fisher, “Ubiquitous signatures of nematic quantum criticality in optimally doped Fe-based superconductors,” Science 352, 958–962 (2016).
  * Cui and Fernandes (2018) Tianbai Cui and Rafael M. Fernandes, “Smeared nematic quantum phase transitions due to rare-region effects in inhomogeneous systems,” Physical Review B 98, 085117 (2018).
  * Griffiths (1969) Robert B. Griffiths, “Nonanalytic Behavior Above the Critical Point in a Random Ising Ferromagnet,” Physical Review Letters 23, 17–19 (1969).
  * Fisher (1992) Daniel S. Fisher, “Random transverse field Ising spin chains,” Physical Review Letters 69, 534–537 (1992).
  * Fisher (1995) Daniel S. Fisher, “Critical behavior of random transverse-field Ising spin chains,” Physical Review B 51, 6411–6461 (1995).
  * Del Maestro  _et al._ (2008) Adrian Del Maestro, Bernd Rosenow, Markus Müller, and Subir Sachdev, “Infinite Randomness Fixed Point of the Superconductor-Metal Quantum Phase Transition,” Physical Review Letters 101, 035701 (2008).
  * Vojta (2003) Thomas Vojta, “Disorder-Induced Rounding of Certain Quantum Phase Transitions,” Physical Review Letters 90, 107202 (2003).
  * Vojta (2006) Thomas Vojta, “Rare region effects at classical, quantum and nonequilibrium phase transitions,” Journal of Physics A: Mathematical and General 39, R143–R205 (2006), 0602312 [cond-mat] .
  * Vojta (2010) Thomas Vojta, “Quantum Griffiths Effects and Smeared Phase Transitions in Metals: Theory and Experiment,” Journal of Low Temperature Physics 161, 299–323 (2010).
  * Ubaid-Kassis  _et al._ (2010) Sara Ubaid-Kassis, Thomas Vojta, and Almut Schroeder, “Quantum Griffiths Phase in the Weak Itinerant Ferromagnetic Alloy Ni1-xVx,” Physical Review Letters 104, 066402 (2010).
  * Wang  _et al._ (2017) Ruizhe Wang, Adane Gebretsadik, Sara Ubaid-Kassis, Almut Schroeder, Thomas Vojta, Peter J. Baker, Francis L. Pratt, Stephen J. Blundell, Tom Lancaster, Isabel Franke, Johannes S. Möller, and Katharine Page, “Quantum Griffiths Phase Inside the Ferromagnetic Phase of Ni1-xVx,” Physical Review Letters 118, 267202 (2017).
  * Seidler  _et al._ (1992) G. T. Seidler, T. F. Rosenbaum, and B. W. Veal, “Two-dimensional superconductor-insulator transition in bulk single-crystal YBa2Cu3O6.38,” Physical Review B 45, 10162–10164 (1992).
  * Schneider  _et al._ (2012) R. Schneider, A. G. Zaitsev, D. Fuchs, and H. v. Löhneysen, “Superconductor-Insulator Quantum Phase Transition in Disordered FeSe Thin Films,” Physical Review Letters 108, 257003 (2012).
  * Shi  _et al._ (2014) Xiaoyan Shi, Ping V. Lin, T. Sasagawa, V. Dobrosavljević, and Dragana Popović, “Two-stage magnetic-field-tuned superconductor-insulator transition in underdoped La2-xSrxCuO4,” Nature Physics 10, 437–443 (2014).
  * Xing  _et al._ (2015) Ying Xing, H.-M. Zhang, H.-L. Fu, Haiwen Liu, Yi Sun, J.-P. Peng, F. Wang, Xi Lin, X.-C. Ma, Q.-K. Xue, Jian Wang, and X. C. Xie, “Quantum Griffiths singularity of superconductor-metal transition in Ga thin films,” Science 350, 542–545 (2015).
  * Saito  _et al._ (2018) Yu Saito, Tsutomu Nojima, and Yoshihiro Iwasa, “Quantum phase transitions in highly crystalline two-dimensional superconductors,” Nature Communications 9, 778 (2018).
  * Lewellyn  _et al._ (2019) Nicholas A. Lewellyn, Ilana M. Percher, Jj Nelson, Javier Garcia-Barriocanal, Irina Volotsenko, Aviad Frydman, Thomas Vojta, and Allen M. Goldman, “Infinite-randomness fixed point of the quantum superconductor-metal transitions in amorphous thin films,” Physical Review B 99, 054515 (2019).
  * Liu  _et al._ (2019) Yi Liu, Ziqiao Wang, Pujia Shan, Yue Tang, Chaofei Liu, Cheng Chen, Ying Xing, Qingyan Wang, Haiwen Liu, Xi Lin, X. C. Xie, and Jian Wang, “Anomalous quantum Griffiths singularity in ultrathin crystalline lead films,” Nature Communications 10, 3633 (2019).
  * Lin  _et al._ (2015) Yen Hsiang Lin, J. Nelson, and A. M. Goldman, “Superconductivity of very thin films: The superconductor-insulator transition,” Physica C: Superconductivity and its Applications 514, 130–141 (2015).
  * Reiss  _et al._ (2020) Pascal Reiss, David E. Graf, Amir A. Haghighirad, William Knafo, Loïc Drigo, Matthew Bristow, Andrew J Schofield, and Amalia I. Coldea, “Quenched nematic criticality and two superconducting domes in an iron-based superconductor,” Nature Physics 16, 89–94 (2020).
  * Böhmer  _et al._ (2016) A. E. Böhmer, V. Taufour, W. E. Straszheim, Thomas Wolf, and P. C. Canfield, “Variation of transition temperatures and residual resistivity ratio in vapor-grown FeSe,” Physical Review B 94, 024526 (2016).
  * Note (1) See Supplemental Material at [URL will be inserted by publisher] for further experimental data and analyses.
  * Xiang  _et al._ (2017) Li Xiang, Udhara S Kaluarachchi, Anna E Böhmer, Valentin Taufour, Makariy A Tanatar, Ruslan Prozorov, Sergey L Bud’ko, and Paul C Canfield, “Dome of magnetic order inside the nematic phase of sulfur-substituted FeSe under pressure,” Physical Review B 96, 024511 (2017).
  * Bristow  _et al._ (2020a) Matthew Bristow, Pascal Reiss, Amir A Haghighirad, Zachary Zajicek, Shiv J Singh, Thomas Wolf, David E Graf, William Knafo, Alix McCollam, and Amalia I Coldea, “Anomalous high-magnetic field electronic state of the nematic superconductors FeSe1-xSx,” Physical Review Research 2, 013309 (2020a), 1904.02522 .
  * Watson  _et al._ (2015) M. D. Watson, Takuya Yamashita, Shigeru Kasahara, William Knafo, M Nardone, J Béard, Frédéric Hardy, Alix McCollam, A Narayanan, S. F. Blake, Thomas Wolf, A. A. Haghighirad, Christoph Meingast, A. J. Schofield, H. v. Löhneysen, Yuji Matsuda, A. I. Coldea, and Takasada Shibauchi, “Dichotomy between the Hole and Electron Behavior in Multiband Superconductor FeSe Probed by Ultrahigh Magnetic Fields,” Physical Review Letters 115, 027006 (2015).
  * Coldea  _et al._ (2019) Amalia I. Coldea, Samuel F. Blake, Shigeru Kasahara, Amir A. Haghighirad, Matthew D. Watson, William Knafo, Eun Sang Choi, Alix McCollam, Pascal Reiss, Takuya Yamashita, Mara Bruma, Susannah C. Speller, Yuji Matsuda, Thomas Wolf, Takasada Shibauchi, and Andrew John Schofield, “Evolution of the low-temperature Fermi surface of superconducting FeSe1-xSx across a nematic phase transition,” npj Quantum Materials 4, 2 (2019).
  * Shoenberg (1984) D Shoenberg,  _Magnetic oscillations in metals_ (Cambridge University Press, Cambridge, 1984).
  * Paul and Garst (2017) I. Paul and M. Garst, “Lattice Effects on Nematic Quantum Criticality in Metals,” Physical Review Letters 118, 227601 (2017).
  * Wang and Berg (2019) Xiaoyu Wang and Erez Berg, “Scattering mechanisms and electrical transport near an Ising nematic quantum critical point,” Physical Review B 99, 235136 (2019).
  * de Carvalho and Fernandes (2019) V. S. de Carvalho and Rafael M Fernandes, “Resistivity near a nematic quantum critical point: Impact of acoustic phonons,” Physical Review B 100, 115103 (2019).
  * Vieira  _et al._ (2019) Lucas E. Vieira, Vanuildo S. de Carvalho, and Hermann Freire, “DC resistivity near a nematic quantum critical point: Effects of weak disorder and acoustic phonons,” Annals of Physics 419, 168230 (2019).
  * Fisher (1990) Matthew P. A. Fisher, “Quantum phase transitions in disordered two-dimensional superconductors,” Physical Review Letters 65, 923–926 (1990).
  * Terashima  _et al._ (2014) Taichi Terashima, Naoki Kikugawa, Andhika Kiswandhi, Eun-Sang Choi, James S. Brooks, Shigeru Kasahara, Tatsuya Watashige, Hiroaki Ikeda, Takasada Shibauchi, Yuji Matsuda, Thomas Wolf, Anna E. Böhmer, Frédéric Hardy, Christoph Meingast, Hilbert v. Löhneysen, Michi-To Suzuki, Ryotaro Arita, and Shinya Uji, “Anomalous Fermi surface in FeSe seen by Shubnikov-de Haas oscillation measurements,” Phys. Rev. B 90, 144517 (2014).
  * Bristow  _et al._ (2020b) M Bristow, A. A. Haghighirad, M. D. Watson, P. Reiss, Z. Zajicek, J. Prentice, S. J. Blundell, A. McCollam, and A. I. Coldea, “Multi-band effects and dominant inter-band pairing responsible for the upper critical fields of bulk FeSe,” in preparation (2021b).
  * Farrar  _et al._ (2020) Liam S Farrar, Matthew Bristow, Amir A Haghighirad, Alix McCollam, Simon J Bending, and Amalia I Coldea, “Suppression of superconductivity and enhanced critical field anisotropy in thin flakes of FeSe,” npj Quantum Materials 5, 29 (2020).
  * Coldea (2020) Amalia I. Coldea, “Electronic nematic states tuned by isoelectronic substitution in bulk FeSe1-xSx,”  arXiv:2009.05523 (2020).
  * Bristow  _et al._ (2020c) M. Bristow, P. Reiss, A. A. Haghighirad, Z. Zajicek, S. J. Singh, T. Wolf, D. Graf, W. Knafo, A. McCollam, and A. I. Coldea, “Anomalous high-magnetic field electronic state of the nematic superconductors ${\mathrm{FeSe}}_{1-x}{\mathrm{S}}_{x}$,” Phys. Rev. Research 2, 013309 (2020c).
  * Matsuura  _et al._ (2017) Kohei Matsuura, Yuta Mizukami, Y. Arai, Y. Sugimura, N. Maejima, A. Machida, T. Watanuki, T. Fukuda, T. Yajima, Z. Hiroi, K. Y. Yip, Y. C. Chan, Q. Niu, Suguru Hosoi, Kousuke Ishida, K. Mukasa, Shigeru Kasahara, J.-G. Cheng, Swee K Goh, Yuji Matsuda, Yoshiya Uwatoko, and Takasada Shibauchi, “Maximizing $T_{c}$ by tuning nematicity and magnetism in FeSe1-xSx superconductors,” Nature Communications 8, 1143 (2017).
  * Reiss  _et al._ (2017) Pascal Reiss, Matthew D Watson, Timur K Kim, Amir A Haghighirad, D. N. Woodruff, Mara Bruma, S. J. Clarke, and Amalia I Coldea, “Suppression of electronic correlations by chemical pressure from FeSe to FeS,” Physical Review B 96, 121103(R) (2017).
  * Seppälä and Alava (2001) E. T. Seppälä and M. J. Alava, “Susceptibility and percolation in two-dimensional random field Ising magnets,” Physical Review E 63, 066109 (2001).
  * Böhmer  _et al._ (2015) A E Böhmer, F Hardy, L Wang, T Wolf, P Schweiss, and C Meingast, “Superconductivity-induced re-entrance of the orthorhombic distortion in Ba1-xKxFe2As2,” Nature Communications 6, 7911 (2015).
  * Hosoi  _et al._ (2016) Suguru Hosoi, Kohei Matsuura, Kousuke Ishida, Hao Wang, Yuta Mizukami, Tatsuya Watashige, Shigeru Kasahara, Yuji Matsuda, and Takasada Shibauchi, “Nematic quantum critical point without magnetism in FeSe1-xSx superconductors,” Proceedings of the National Academy of Sciences 113, 8139–8143 (2016), 1604.00184 .