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

We combine ab initio simulations and Raman scattering measurements to demonstrate explicit anharmonic effects in the temperature dependent dielectric response of a NaCl single crystal. We measure the temperature evolution of its Raman spectrum and compare it to both a quasi-harmonic and anharmonic model. Results demonstrate the necessity of including anharmonic lattice dynamics to explain the dielectric response of NaCl, as it is manifested in Raman scattering. Our model fully captures the linear dielectric response of a crystal at finite temperatures and may therefore be used to calculate the temperature dependence of other material properties governed by it.


I. INTRODUCTION
The microscopic manifestation of temperature in lattice dynamics is generally understood through the harmonic approximation, where atomic motion is mapped onto a set of vibrational normal modes (i.e., phonons) [1].Modern studies of the lattice dynamics of solids significantly expand on this view to account for anharmonic phenomena such as phase transitions [2,3], thermal expansion [4], and thermal conductance [5][6][7].
In contrast, studies that focus on the electronic properties of solids rarely go beyond a quasiharmonic (QH) treatment, where the normal modes are renormalized for each given temperature to account for thermal expansion [8].The effects of lattice dynamics are usually incorporated through the prism of electron-phonon interaction.This can manifest in modified carrier lifetimes [9,10], corrections to the electronic band structure [11], mobility estimations [12][13][14], and calculations for vibrational spectroscopy [15][16][17][18].Formally, almost any ab initio calculation of electron-phonon interactions relies on the harmonic approximation [19].A different manifestation of lattice dynamics in electronic properties is the temperature evolution of the material's dielectric response.The need to extend anharmonic modeling to dielectric properties has been demonstrated perhaps most strikingly in halide perovskite semiconductors [20][21][22][23][24][25][26].Importantly, by relying on molecular dynamics (MD) simulations, calculations of the dielectric response of a crystal are not limited to a perturbative treatment of anharmonicity.
Because inelastic light scattering originates in the vibrational modulation of a crystal's polarizability autocorrelation function [27], it is the ideal probe to study the effect of anharmonicity on dielectric properties.A main drawback of optical inelastic light scattering is that momentum conservation limits vibrational contributions to the , or zero crystal momentum point.This limitation does not apply in secondorder Raman scattering, which is driven by contributions from the entire Brillouin zone (BZ) [28].Simply put, second-order Raman describes the inelastic scattering of a single photon with two phonons [29].The scattering intensity inside a given frequency interval I ( )d emerges out of contributions from all possible phonon combinations with the appropriate energy and momentum conservation.It is therefore sensitive to any anharmonic effect in either the phonon dispersion or the dielectric response.
In this study, we introduce an ab initio approach incorporating finite-temperature effects into the dielectric response of a crystal.The approach is based on the effective interatomic force constants used in the temperature-dependent effective potential (TDEP) method [30][31][32].The computational method is benchmarked with experimental measurements of dielectric response in terms of second-order Raman scattering.As a model system, we choose NaCl.It is an ideal showcase since first-order Raman scattering is forbidden by symmetry, necessitating a computational methodology that goes beyond the lowest order.Our treatment is not compromised by the natural isotopic distribution in NaCl, as discussed in Appendix E. Furthermore, the central role of anharmonicity in NaBr was recently demonstrated [33], motivating further examination of the rock-salt structure.
The comparison between simulated and experimental results shows that the harmonic treatment is inadequate when considering dielectric response at finite temperatures.In contrast, the method developed here not only reproduces experimental results more faithfully than conventional approaches, it also shows how a higher-order dielectric response is linked to other manifestations of anharmonicity.More specifically, it accounts for phonon broadening and thermal transport, and demonstrates how a measurement of the second-order Raman spectrum can be used to conclude the dominant anharmonic expressions in a given material.Importantly, due to its ab initio nature, our method provides a full description of the linear dielectric response of the material and its dependence on atomic displacements.

A. Experimental and calculated spectra
Our main results are presented in Fig. 1, where we compare the measured, temperature-dependent, unpolarized Raman spectra of NaCl and corresponding simulations.The technical details for the construction of an unpolarized spectrum are provided in Sec.SI in the Supplemental Material (SM) [34].Details of the corresponding computations are given in Appendix C. Our experimental spectra are in agreement with past measurements [35,36].We observe a continuous spectrum with some sharp non-Lorentzian features centered around 280 cm −1 .
Unlike early analytical treatments [37][38][39], our computation is able to predict the entire shape of the spectrum.It is immediately evident that even at 300 K, the QH calculation for the spectrum is imperfect.It overestimates the high-frequency (≈300 cm −1 ) to low-frequency (≈100 cm −1 ) scattering intensity ratio.It also fails to capture the broad decaying intensity above 350 cm −1 .The discrepancies between experiment and QH calculation at 300 K are remarkable, as they demonstrate the significance of anharmonic effects in the dielectric response of an archetypal rock-salt crystal, even at room temperature.These discrepancies get even larger at higher temperatures.Above 500 K, the main spectral feature around 300 cm −1 broadens and flattens, but the QH calculation predicts an opposite trend.These findings constitute direct evidence for the failure of a QH approach to explain the dielectric response of NaCl, due to finite-temperature anharmonic effects.
Contrary to the QH approach, the TDEP-based computation shows much better agreement with experiment.At 300 K, the line shape of the TDEP-based Raman spectrum (blue) more accurately follows the experimental spectrum (black).The agreement between experimental and TDEP spectra is even better above 500 K.The main discrepancy remains the decaying intensity above 350 cm −1 which is not captured by either QH or TDEP.We suspect this is the result of higherorder (>2) scattering terms since, as discussed below, our calculation does not include them.Having established that our computational method is both required and successful in predicting the Raman spectrum, we turn to discuss the TDEP-based calculation in more detail.

B. Dielectric response of a crystal from first principles
We start by briefly reiterating the main idea behind TDEP to describe the ionic motion.The starting point is the lattice dynamical expansion for the crystal energy [1], where u α i denote the displacement of atom i in direction α.The interatomic force constants are determined by minimizing the difference between the model system and ab initio calculated forces f , = arg min f model − f ab initio . ( The ab initio calculations are performed in the canonical ensemble, yielding an explicit temperature-dependent model Hamiltonian-either via MD or a self-consistent stochastic sampling [40].In this study, we used the latter.Since the sampling is done at finite temperature, an explicit temperature dependence is built into the effective interaction parameters.The TDEP Hamiltonian is constructed to be the best possible fit at a given temperature. To assess the accuracy of our model Hamiltonian for the lattice dynamics, we calculate the phonon spectral function and compare it with a neutron-scattering experiment [41].The calculated spectral functions agree well with experiments and previous nonharmonic calculations [42] (see Sec. SIV B in the SM [34]).
The need for a nonperturbative treatment when describing anharmonic lattice dynamics motivates an analogous treatment for the dielectric response.To extend the formalism to dielectric properties, we similarly expand the polarizability P and dipole moment M in terms of atomic displacements u [1], Here the indices μν denote Cartesian components of electric field derivatives.In principle, this expansion is done for dynamic quantities-including the electronic frequency dependence-but in the present case for NaCl, we are interested in a frequency range far from any resonances and can safely ignore the electronic frequency dependence.
In analogy with the process described for the force constants, the terms in the dipole moment and polarizability expansions are determined by minimizing the difference between the model values and those obtained from simulation, The end result is a set of interaction parameters that explicitly depends on temperature, incorporating all orders of nonharmonic effects in the linear dielectric response of a material.We note that the inclusion of frequency dependence adds no conceptual difficulty, only the practical challenge of accurately determining the frequency-dependent electronic dielectric response for a simulation with hundreds of atoms.The algorithm for determining the interaction parameters is described in detail in Appendix B, and the specific computational details are given in Secs.SIII [43] and SIV [40,[44][45][46][47][48][49][50][51][52][53][54][55] in the SM [34].
With the expansion of the dipole moment and polarizability at hand, a host of physical quantities may potentially be calculated.The dielectric response of a material governs numerous material properties such absorption, reflection, beam propagation [56], capacitance [57], and ferroelectricity [58], all of which are determined by the dielectric function.Here, we are interested in the inelastic light scattering character of the material and therefore turn to compute the Raman scattering cross section from the dielectric response obtained above.

C. Ab initio calculation of the Raman spectrum of a crystal
Theoretical interpretations of second-order Raman spectra usually limit themselves to high-symmetry areas of reciprocal space where most density of states (DOS) singularities are expected [59][60][61][62][63][64].This approach adequately accounts for pronounced peaks, but severely limits the possibility of incorporating anharmonic effects.Also, because of the symmetry breaking of moving away from the point and generalization into multiple phonon states, selection rules in second-order Raman scattering usually become too relaxed to place meaningful restrictions on the measured Raman tensor (see Appendix D for more details).Existing ab initio approaches to second-order Raman scattering all operated strictly within the harmonic approximation [15,16].A shell-model-based calculation by Bruce [65] did incorporate phonon-phonon interactions, but his partial treatment led him to conclude these are negligible in the Raman spectrum of NaCl.
The construction of the full second-order spectrum from first principles requires a few computational steps.The starting point is to use the TDEP method to construct an anharmonic Hamiltonian that describes the dynamics of the ions.This anharmonic Hamiltonian is then coupled to an expansion of the polarizability in terms of atomic displacements.The coupled ionic-polarizability system can then be solved with many-body methods, providing the finite-temperature spectrum.The steps involved are discussed in detail below.
Following Cowley [45], the Raman scattering cross section is given by where E in and E out are the electric field vectors of the incoming and outgoing light, with Greek letters standing for Cartesian components.The tensor I depends on the probing frequency and constitutes the material property governing Raman scattering.It is given by which is the Fourier-transformed thermal average of the polarizability-polarizability autocorrelation function.
The full tensorial formalism for second-order Raman scattering has already been rigorously worked out, but is rarely used and is repeated and extended here [66,67].Far from resonance, to lowest order, the tensor I relating incident and outgoing light has four terms: with the different terms given by where the matrix elements P (I) − P (IV) are defined via the Fourier components of the terms in Eq. ( 3) (see Appendix C for the explicit expressions and derivations of the above terms).The phonon spectral function for mode s, J s = −Im{G s ( )}/π , is obtained from where qs (Z ) = −18 and The first term, Eq. ( 10), is the first-order Raman scattering which comes down to the one-phonon spectral function weighted by the first-order Raman matrix elements.The spectral function contains any broadening and shifts due to anharmonicity, and also any deviation from a Lorentzian line shape.This term is what is commonly referred to as Raman scattering.
The second term, Eq. ( 11), is the second-order Raman scattering.The matrix elements contain the momentum conservation, and the convolution term contains the energy conservation.If one were to set the matrix elements to unity, it would yield a spectrum that follows the two-phonon DOS.
The third and fourth terms are more subtle, as they contain contributions from the whole BZ but are multiplied with the one-phonon line shape.They will thus decay rapidly away from the first-order Raman peaks, with the net effect of slightly shifting and altering the shapes of the first-order peaks.These terms are one reason that neutron-scattering and Raman scattering measurements might not coincide exactly.
In NaCl, first-order Raman scattering is forbidden by symmetry, which means P (I) , P (III) , and P (IV) are zero.This statement remains valid for isotopic impure crystals of the rock-salt structure (such as NaCl), as discussed further in Appendix E. We stress that the spectrum is calculated from the interacting phonons, i.e., a convolution of one-phonon spectral functions, and includes explicit anharmonic effects.
Given an ab initio calculation, it is possible to track the contributions of each phonon-pair component in Eq. ( 11) and identify the specific ones that are more significant to the spectrum.Our analysis of NaCl showed that all phonon branches and BZ parts contributed appreciably, so no phonon pair proved noteworthy beyond the expected DOS dependency, demonstrating the necessity for summation over the entire BZ to produce accurate predictions.
Equation ( 8) says that the complete description for Raman scattering is given by a tensor detailing the vectorial relationship between incident and scattered fields.It is therefore desirable to inspect polarization dependence when comparing theory and measurement.

D. Polarization dependence of the Raman spectrum of NaCl
To realize a detailed comparison between our ab initio calculation and experiment, we performed a polarizationorientation (PO) Raman measurement of a NaCl single crystal.In a PO measurement, the polarization of the incident laser beam is rotated in a plane parallel to the NaCl (001) crystal face.We separately probe the parallel and perpendicular scattered polarizations with respect to the incident beam orientation.A diagram depicting the different polarizations involved in the measurement is presented in Fig. 2.
The temperature-dependent spectra presented in Fig. 1 are actually calculated sums of the scattered intensities for a full rotation of incident polarization.Details of the experimental setup are given in Appendix A.
Our theoretical calculation also fully accounts for the tensor nature of Raman scattering, making a detailed comparison possible.Figure 3 shows the full PO color map for measurement and simulation in 400 K (PO maps for all temperatures are given in Sec.SI in the SM [34]).The theoretical PO maps in Fig. 3 are calculated from Eq. (11).Dashed colored lines follow corresponding spectral features in simulation and experiment.The polarization angle (vertical axis) is offset by a random angle controlled by the crystal orientation in the laboratory frame.We can actually use the observed PO pattern to determine this orientation (up to a π/4 rotation).Specifically for O h structures in the (001) backscattering geometry, the fourfold period peaks must bisect or align with the crystallographic axes.
The calculated model closely follows the PO dependence throughout the spectrum.Our calculation correctly predicts the periodicity of each feature, as well as the relative phase between them, thereby faithfully reproducing the tensorial nature of the modulated dielectric response.The colored FIG. 3. NaCl Raman PO color map for measurement and TDEP simulation in 300 K. Unpolarized spectra are on top, measured spectra in the middle (Expt.), and calculated spectra on the bottom (Calc.).Parallel polarization (//) maps are to the left; perpendicular (⊥) maps are to the right.The frequency of the spectral features in simulation slightly misses experimental values, so matching dashed color lines help guide the comparison.
lines locations indicate some absolute shift in frequencies between experiment and simulation.Although our calculated dispersion curves match neutron experiments (Fig. S6 of Sec.SIV B in the SM [34,41]), any inaccuracy is doubled by the two-phonon scattering mechanism.Accordingly, the shift grows larger for higher frequencies.Another discrepancy is the low-frequency (<50 cm −1 ) checkered pattern apparent in all perpendicular configuration measurements up to 600 K (Sec.SI in the SM [34]), not reproduced by the calculated spectra.This might stem from flat regions in the effective surface potential, which demand high-order terms absent in our polynomial expansion of Eq. (1).We already speculated these to have observable contributions in the high-frequency (>350 cm −1 ) range (Fig. 1), but they may also affect low frequencies through difference combinations.

E. Anharmonic effects in second-order Raman scattering
It is worth dwelling on the temperature effects in secondorder Raman scattering, and their connection to anharmonic behavior.The scattering is generated by phonons throughout the BZ.Since the TDEP method implicitly contains all orders of anharmonicity, this already imparts a nontrivial temperature dependence, but the effects of anharmonicity go beyond the temperature dependence of the bare phonon dispersion curve.Inspecting Eq. (11), one sees that it is a convolution of two one-phonon spectral functions.If we replace the interacting spectral functions with the noninteracting ones, we will recover the traditional result [66] [see Eq. (C29)].In the present study, we found it crucial to include the interacting spectral functions when calculating the second-order Raman spectra; without them, agreement with experiment was markedly worse.Intuitively, if the one-phonon DOS is broadened, the two-phonon DOS must also broaden-and since the two-phonon DOS is a convolution, non-Lorentzian line shapes have a larger impact on two-phonon than on one-phonon spectra.
Besides its own nontrivial temperature dependence, second-order Raman scattering emerged as a uniquely sensitive probe for anharmonic lattice dynamics.Three-phonon scattering is ostensibly a straightforward perturbative anharmonic effect.Only a small fraction of all possible three-phonon combinations satisfies energy and momentum conversion-this is what we refer to as the scattering phase space.Mathematically, this is defined via the imaginary part of Eq. ( 16).The temperature dependence of available phonon scattering phase space is, however, a nonperturbative anharmonic effect and can lead to a host of surprising trends with temperature.Thermal conductivity that does not follow T −1 [6], phonon line widths that do not grow linearly with temperature [68], or, in our case, a Raman spectrum that does not follow the quasiharmonic approximation are all consequences of nonperturbative anharmonicity.
The second-order Raman spectrum follows the same set of momentum and energy conservation rules as three-phonon scattering.In fact, the spectra are directly proportional to the available scattering phase space.This makes Raman scattering invaluable as a direct probe of anharmonicity: we can see how (a part) of the scattering phase space evolves with temperature, providing insight into the underlying mechanisms governing strongly anharmonic materials.In this specific case, we can verify prior theoretical studies.Ravichandran and Broido [42] found it necessary to include temperature-dependent phonons to accurately describe thermal transport in NaCl-this can be corroborated experimentally by the change in the Raman spectra with temperature in Fig. 1.Moreover, four-phonon scattering was predicted to be important.The presence of higherorder scattering is evident in the tail at large wave numbers (second-order Raman scattering can only contribute up to wave numbers equal to 2ω max , and third-order Raman, which FIG. 4. Temperature dependence of the Born effective charge and ∞ in NaCl treated in three different ways.The most realistic, at constant pressure (i.e., volume changing with temperature) with full temperature dependence of all parameters, is denoted NPT.The same treatment without thermal expansion is denoted NVT.For reference, we also include the quasiharmonic results, denoted QH, where all interaction parameters are set to their values at 0 K, and temperature is only included via thermal expansion.
would involve four-phonon scattering, can contribute up to 3ω max ).

III. TEMPERATURE DEPENDENCE OF ADDITIONAL DIELECTRIC PROPERTIES
The methodology described above is general with regard to the interplay between anharmonicity and dielectric response.Experimentally, we have the Raman spectra to verify our results with.To illustrate the general nature of the method, we calculate a few other quantities accessible by it.We present the infrared absorption spectra, as well as the temperature dependence of the dielectric tensor and Born charges.
In Fig. 4, we show the temperature dependence of the Born charges and static dielectric tensor.The temperature dependence is modest as one would expect in a wide-bandgap material.Even so, the quasiharmonic treatment of Born charges overestimates the temperature dependence by a factor of 2 when compared with the TDEP calculation.It is easy to imagine that in a material with a much smaller band gap, this discrepancy could be significantly larger.Even more striking is the difference between the NPT (fully anharmonic at constant pressure) and NPC (fully anharmonic at constant volume) curves.This demonstrates that thermal expansion alone cannot give a full account of temperature evolution, with anharmonic effects giving rise to an opposite competing trend.
In Fig. 5, we show the IR spectra of NaCl for a few temperatures.The procedure to calculate these is identical to the Raman spectra, the only difference being the matrix elements.Past measurements of the IR spectrum of NaCl are consistent with our results (see Appendix E for a discussion regarding the isotopic effects absent from our calculation), but only covered a limited range of frequencies [69,70], so we leave the full spectra as predictions.

IV. CONCLUSIONS
In conclusion, we have demonstrated the crucial role of anharmonic lattice dynamics effects in determining the phonon-modulated linear dielectric response of a crystal.This was achieved by introducing an ab initio method to calculate the dielectric response of a crystal at finite temperatures and comparing it to the first continuous PO spectrum of an exclusively second-order Raman structure.Our generalized TDEP method incorporates all orders of nonharmonic effects through a sampling of temperature-dependent lattice configurations, expressed as effective interatomic force constants.We introduce analogous effective interaction tensors that govern the dielectric response, i.e., we expand the dipole moment and polarizability in terms of atomic displacements.Using the new generalized TDEP method, other dielectric material properties such as ferroelectricity and capacitance may be calculated, along with their temperature dependence.By isolating and evaluating the normal modes comprising the emergent temperature dependence of a material property, new design rules may be inferred for better functional materials operating at finite temperatures.spontaneous emission by 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. 6) 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. 6).The polarized beam then enters a microscope and is 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 Sec. SII in the SM [34]) is placed inside a (TS1000 Hi-Temp) Linkam stage under argon purge.Purging with a monoatomic gas proved crucial since under vacuum local laser heating would create defects, while nitrogen would dominate the Raman signal with its rotational excitations.A FIG. 6.Schematic of experimental setup.Laser light enters from top.Scattered signal goes right for detection.monoatomic 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 backscattered light is collected by the same objective and passes back through HW1 and through the BS, where 90% of the Rayleigh is eliminated and the Raman signal is reflected down the beam path (to the right in Fig. 6).Note that the beam's polarization at this point is not known since we cannot assume the scattered light, especially the Raman signal, scattered back into the same incident polarization.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 also be polarized perpendicular to the incident polarization.After the BS, the scattered beam (orange line in Fig. 6) passes through an achromatic half-wave plate (HW2 in Fig. 6) which either leaves or rotates the signal polarization by 90 • .Another identical polarizer (designated "analyzer" in Fig. 6) filters the polarization of the signal.We thus separate the two scattered polarization components.With HW2 at 0 • , only the parallel signal will make it past the analyzer; with HW2 at 45 • , only the perpendicular signal is transmitted.The system achieved an extinction ratio (I ⊥ /I ) 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 ultranarrowband) notch filters and the signal is directed into a one meter (Horiba FHR 1000) spectrometer and detected by a (Horiba Synapse) CCD.

APPENDIX B: 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 polarizability interactions are determined and how interatomic force constants are determined.The starting point is a set of displacements u and forces 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 × 1 vectors u and f (for clarity, we will make note of the dimensions of the matrices in each step): This can equivalently be written as with the Kronecker product ⊗.To express the forces from N c supercells, the matrices are stacked on top of each other, ⎛ ⎝ f 1 . . .
Without any loss of generality, the second-order force constants can be expressed as a linear combination of N x irreducible components x II , where N x (3N a ) 2 , via a matrix The set of symmetry relations and invariances determines C and N x , [1,44,71] and in general N x (3N a ) 2 which considerably simplifies the problem of numerically determining x.In practice, only the matrix A is determined and stored.The nominally very large matrices that go into the construction of 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 where the contracted matrix A is many orders of magnitude smaller than the matrices it is built from.The interatomic force constants are determined in succession, where f denotes the reference forces to be reproduced.Equations (B6)-(B8) 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 ∂ 2 U/∂E 2 as well as ∂ 2 U/∂E ∂R, which are 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 The terms in the polarizability expansion are determined slightly differently than the interatomic force constants.For a set of N c supercells, we build the set of differences in polarizability between supercells i and j, and solve for x PI from where the set of all pairs of supercell i j is used.The fluctuations from the first order are removed to construct and so on for all orders considered.The baseline polarizability is then determined via where P th is the contribution to the polarizability from the second-order terms in the polarizability expansion [66], with the matrix elements 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, where the terms are determined in succession:

APPENDIX C: DETERMINING THE RAMAN SPECTRA
In this section, we outline the procedure to determine the 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 q and polarization s (when there is no loss of clarity, we use the compound index λ to denote qs), Here, are phonon eigenvectors ( i • j = δ i j ), m atomic masses, and ω frequencies.The vector r is a lattice vector, and u and p are the position and momentum operators.We use the notation q = −q.We introduce scaled eigenvectors (to simplify notation) via In this notation, the matrix elements pertaining to anharmonicity become so that the anharmonic part of the Hamiltonian becomes In an analogous way, we define the coefficients of the expansion of the polarizability as = P μν (q, s, s ) = P μν (q, s, s ) † = P μν (q, s , s), (C10) and the dipole moment matrix elements become such that

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, we express the equations of motion of the phonon-phonon Green's function via and the four-point Green's functions are decoupled [72,73] via abcd ≈ ab cd + ac bd + ad bc to yield After Fourier transforming equations (C18a)-(C18i) to the frequency domain, we get a solvable set of equations that recover the usual definition of three-phonon anharmonicity [44][45][46]72], where λλ (Z ) = −18 and 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,

Raman spectra
If we start from the polarizability-polarizability correlation function and insert the expansion of the polarizability in terms of phonon coordinates, given by Eq. (C15), we get 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 where J is the phonon spectral function, We identify Eq. (C24) as the first-order Raman spectra.The next contribution comes from The first step is to decouple the four-point Green's function [73], and note that the Fourier transform of a product becomes a convolution, where we relabeled 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 noninteracting spectral function With this, we recover the expression given by Cowley [66] 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: second-order Raman scattering 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 where we make use of the intermediate result for the three-point Green's functions in Eq. (C22) and arrive at where, if we assume a diagonal self-energy, we recover previous results [66].The final term, we approach with the same decoupling procedure as before: where the factor 3 comes from relabeling indices.We also used the equal-time thermal average of the phonon operators.This gives It is worth noting that as far as we have discovered, the last two terms I III , I IV 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 onephonon spectral function, so to have any sizable contribution these need to peak simultaneously.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 (its 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 polarizability is of interest, it is conveniently calculated via a Kramers-Kronig transformation of the imaginary part.

APPENDIX D: FACTOR GROUP ANALYSIS FOR SECOND-ORDER RAMAN SCATTERING
The full factor group analysis for the rock-salt structure has already been worked out by Burnstein [59], following methods established by Birman [37] and Elliot [74].These ostensibly provide a complete toolset to interpret second-order Raman spectra by assigning each spectral feature to a known combination of vibrational modes (besides having a Raman active irrep component, the phonon pair must also conserve momentum q 1 = −q 2 ).As mentioned above, this treatment may prove useful for mode assignment of 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 second-order Raman scattering are too relaxed to offer complete predictive power, even with reliable dispersion relations at hand.
To demonstrate this, we perform the factor group analysis of one high-symmetry point in the Brillouin zone; specifically, for the X 5 (TO), X 5 (TA) phonon pair in the reciprocal space X-point of a rock-salt structure.
The selection rule governing Raman scattering is [27] f ⊆ 15 ⊗ 15 , (D1) with f the representation of the excited phonon state, and 15 the representation for the momentum operator in BSW notation (Raman scattering can be viewed as the product of two dipole transitions).The crux of the matter lies in the fact that in second-order Raman scattering, this final phonon state is generally composed of two off -point phonons, with its representation given by the Kronecker product of the two single-phonon states.For example, the quantum state |X 2phonon = |X 5 (TO), X 5 (TA) has the following corresponding representation: with all the irreducible representations in the direct sum being Raman active according to Eq. (D1).This leads to the symmetry-allowed Raman tensor form [75], 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 ( )d from a single-phonon combination.In practice, Eq. ( 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 exacerbates for lower-symmetry structures, where more phonon combinations will be Raman active and have more general Raman tensors.

APPENDIX E: EFFECTS OF ISOTOPIC IMPURITY IN THE OPTICAL SPECTRA OF ROCK-SALT CRYSTALS
Generally, isotopic impurities in crystals can have a significant implication on its optical spectra.Specifically, they contribute to the anharmonic character of its structural dynamics by modifying the self-energy of vibrational modes [76].From a theoretical point of view, isotopic effects may modify the predicted spectra for a crystal through the three main mechanisms described below.
First, like a mass on a spring, a change of mass modifies the vibration frequency (ω −2 ∝ M ).In the case of a crystal, isotopes will change the phononic dispersion curve.This can introduce a shift to the observed optical spectra and has been observed, for example, in silicon [77].Since our simulation used the natural isotopic distribution of NaCl (Cl 3 5 75.53%,Cl 3 7 24.47%[78]), our spectral functions include this kind of shift.
Stated more generally, the distribution of isotopes modifies the dynamical matrix for the system.This suggests its eigenvector can change as well as its eigenvalue, and with it the mode's symmetry and irrep assignment.A potential second mechanism for spectra modification is therefore the "turning on" of previously forbidden transitions.However, in the case of the rock-salt structure, its first-order Raman inactivity may be understood through its vanishing real-space polarizability derivatives.The isotopic distribution in NaCl does not alter the equilibrium structure potential surface, so all point-group symmetry operations remain valid, and the symmetry constraints placed on tensors such as (C9)-(C11) are left unaffected.
Strictly speaking, isotopic substitutions reduce the symmetry of the system from the crystal space group to the point group of a lattice site [79].This breakdown in translational invariance leads to a partial lifting of momentum conservation [76] and offers a third mechanism for spectra augmentation.Accordingly, first-order Raman spectra may acquire contributions from the entire BZ and display a component reminiscent of the phonon DOS [80].How this kind of relaxation will affect second-order spectra, where the entire BZ is already incorporated, is not obvious.Measurements of isotopically controlled diamond [69,81] and silicon [77] crystals suggest it to be negligible, with the isotopic effect in their second-order Raman spectra adequately explained by the mass renormalization described above.Note that this line of argumentation does not hold for the IR spectrum of NaCl.Harmonic selection rules dictate that the rock-salt structure should be first-order IR active, and relaxation of momentum conservation has indeed been reported to induce small DOS components in the absorption spectrum of NaCl [69,70].
Isotopic effects are not expected to have a strong temperature dependence compared to non-Hooke-like interactions in the second-order Raman spectrum of NaCl.Our calculation accounted for the isotopic role in the self-energy, as detailed in Sec.SIV C in the SM [34,47,78], and found it to be small compared to temperature-activated effects.This was also corroborated by our measurements, which demonstrated a dramatic evolution in the spectrum with temperature.