Anharmonic effects on the reflectivity of CaS and MgS: a first-principles based study

We employ systematic calculations based on density functional theory to model the reflectivity of CaS and MgS in the infrared region. We show that in addition to the modeling using the harmonic approximation, an accurate spectral description requires the inclusion of anharmonic effects. Due to their conceptual simplicity, CaS and MgS are excellent systems for the explicit consideration of the anharmonicity, which we include here using a perturbative approach up to three-phonon scattering processes, and the consideration of isotopic disorder. All physical quantities, such as Born effective charges and dielectric constant, necessary for the calculation of the reflectivity within the Lorentz model are extracted from our first-principles computations. To validate our predicted optical and transversal modes, and reflectivity spectra, we compare them to available experimental results. We find that the overall agreement is good, which supports the importance of the inclusion of anharmonic terms in the modeling of optical properties in the infrared region.


I. INTRODUCTION
The alkaline-earth metal sulfides XS (X = Be, Mg, Ca, Sr, Ba) have recently attracted interest in science and technology because of their remarkable physical properties and wide applications, ranging from catalysis to microelectronics 1 .In particular, CaS-and MgS-based compounds are known for their extensive use as highly efficient photoluminescents, cathodoluminescents, and Xray phosphors [2][3][4][5] .Additionally, the alkaline-earth and transition metal sulfides have been investigated for a potential use as electrode materials for electrochemical energy storage [6][7][8][9][10] .Also, the optical response of MgS and CaS in the radio and infrared range is of high interest in the study of Mercury's surface composition as previous research and surveys of the planet speculate that it is rich in volatiles 11,12 .Lastly, these monosulfides can also be found in meteorites [13][14][15][16] , therefore their thorough understanding could help to better describe the evolution of our solar system.
There have been a few experimental and theoretical studies of the optical properties in the infrared region of CaS and MgS.From experimental investigations, data obtained from absorption spectra are used to determine various physical quantities such as reflectivity, emissivity 17,18 , etc. On the other hand, theoretical studies based on density functional theory (DFT) within the local density approximation (LDA) 19 have only modeled the atomic dynamics of these systems in the limit of the harmonic approximation 20,21 .Unfortunately, the lack of neutron scattering studies on these sulfides makes impossible a direct comparison of the predicted phonon dispersions with their experimental counterparts.The only a) Electronic mail: chmeruk@gfz-potsdam.de b) Electronic mail: mari nv@gfz-potsdam.dephysical quantity directly comparable from experimental and modeling results is the splitting between longitudinal and transverse optical modes, i.e., the LO/TO splitting, which originates from the degeneracy elimination between the LO and TO phonons at the Brillouin-zone center 22 .In this regard, there is an overall reasonable agreement among reported measurements and predicted data (Table I).However, the harmonic approximation predicts a reflectivity with a step-like behavior and sharp edges at both ends of the spectrum 23 .This predicted reflectivity behavior is in drastic disagreement with experimental spectra, which appear with smeared edges, especially at the high-wavenumber tail.This inconsistency is a direct consequence of the simplifications made in the harmonic approximation and can be remedied by inclusion of anharmonic effects.
Thus, it has long been desired to include anharmonic effects in the first-principles atomic and molecular dynamics (MD) simulations.However, due to computational limitations, the explicit treatment of the phononphonon interactions has been limited, even though the analytical base has been available for a long time 24,25 .But in the last decade, the significant increase in computational power and efficiency has made possible to include anharmonic effects explicitly and on a rigorous basis.One way to consider directly anharmonic effects in a MD simulation is by coupling the system to a thermostat of a given temperature 26,27 .The consequent normal modes of motion obtained then contain the frequencies renormalized by the temperature, i.e., anharmonic interactions.Such an approach has been shown to successfully remedy the negative frequencies in certain elemental systems, persistent within the harmonic approximation 27 .The vibrational eigenfrequencies calculated in this way are "dressed" by anharmonic phononphonon interactions up to, formally speaking, all orders.Therefore, in order to estimate individual contributions, a mapping scheme onto an effective model, containing all the contributions separately needs to be constructed 28 .Alternatively, one might wish to introduce anharmonicity in a progressive and more controlled way.This implementation is achieved in a perturbative manner based on the Cowley's expansion 24 .At the lowest order one obtains two three-phonon and one four-phonon contributions to the total energy as shown in Fig. 1(a)-Fig.1(c).The main computational challenge to evaluate these Cowley's diagrams is the calculation of their force constants (FC) of various orders.Nowadays, it is possible to compute them completely ab intio using density functional theory (DFT) 29,30 by following one of two main approaches.The first method employs density functional perturbation theory (DFPT) 31 , and the FC are expressed directly in terms of the electron density and its derivatives 32 .In the second one, the finite displacement method (FDM) 33,34 , the FC are extracted by performing a series of self-consistent calculations for different atomic displacement patterns and using the Hellman-Feynman theorem.However, the main issue in this approach is that, the computations grow rapidly since a large number of configurations is needed for higher order FC.The calculation of the FC of the m-th order requires (6N n b ) m configurations, where n b is the number of atoms in the lattice basis, and N is the number of primitive cells in the supercell 35 .Although, this number can be greatly reduced by identifying equivalent displacement patterns using crystal symmetries 35,36 , it is still quite large, especially for the fourth-order FC.In addition to Cowley's contributions to the anharmonicity, the effects derived from isotope-disorder scatterings have been shown to play an important role at low temperature 38 .The formulation to deal with this isotopic scattering can also be described in the second-order perturbation theory 37 .In this framework, the unperturbed Hamiltonian, H 0 , for the crystal is defined in the harmonic approximation with the atomic mass replaced by an averaged mass that depends on the number of unit cells in the crystal and the fraction of a given atomic isotope.Then, the first perturbation, H I , which depends on the deviation of the isotope mass from the average mass, has its main contribution from the second-order term shown in Fig. 1(d).
The first-order and other higher-order diagrams derived from H I have either zero contribution (due to the average mass definition) or are not significant in comparison to the second-order term 37 .Therefore, in order to study the most substantial anharmonic effects in MgS and CaS, we use the FDM including only three-phonon processes and the isotopic scattering shown in Fig. 1.
The rest of this article is structured as follows.In Section II we summarize the necessary elements to model the infrared response as well as the means to calculate them.We present our results on optical properties in Section III and compare them to available theoretical and experimental findings within the harmonic limit and highlight the anharmonic contributions.Finally, in Section IV we draw our conclusions and argue for the applicability of the approach used in this work to model infrared optical Phonons are indexed by their wavevector q and branch j.Incoming and outgoing phonons are represented by qj and qj , and internal phonon lines are given by qiji respectively 25 .Diagrams (a) and (b) represent three-phonon scattering processes known as "bubble" and "tadpole" contributions.Diagram (c) shows a four-phonon scattering process known as "loop".In our study, the only process entering in our model, as explained in the main text, is the bubble contribution.Diagram (d) represents isotope-disorder scattering 37 , the dashed lines represent phonon scatterings by the isotope shown here as the black dot.This process is also considered in our anharmonic modeling and it is obtained from second-order perturbation theory, using mass deviation as a perturbation.
properties of other materials.

II. METHODOLOGY AND COMPUTATIONAL DETAILS
Our main goal is to calculate optical response of MgS and CaS in the infrared region and to this end, we use the generalized Lorentz model.In this scheme the dielectric function as a function of frequency ω is given by: where ∞ is the high-frequency dielectric constant, s j and ω j,0 are the strength and frequency of the j-th oscillator, respectively.The oscillator strength is related to the Born effective charges Z, as S j = 4πZ 2 /(vµ j ω 2 j,0 ), where v is the unit-cell volume and µ is the reduced mass of the oscillator 38 .Lastly, Π j (ω) is the self-energy of the j-th vibrational mode.
The self-energy can be improved gradually according to Cowley's perturbative approach 24 .In our work, we restrict its expansion to the following form: where Π 3ph is the anharmonic contribution to the total self-energy due to three-phonon scattering processes, and Π isot is the isotope-disorder induced scattering (Figs.1(a)-(b) and 1(d)).The three-phonon part of the self-energy consists of two terms known as the "bubble" and "tadpole" diagrams.The bubble term has both real and imaginary parts whereas the tadpole has only the real part.Thus, the finite lifetime of the phonons due to the phonon-phonon interaction is solely determined by the imaginary part of the bubble (B) diagram, Im Π B 3ph (ω) , which is the only term considered in our calculation.We note that in general, the bubble diagram is not diagonal in branch indices jj , however, for the purpose of calculating optical properties, it is sufficient to consider only the diagonal terms 25 .The second term in Eq. ( 2) also contributes to the broadening of the harmonic frequencies.Π isot (ω) is mostly visible at low frequencies and temperatures and, as we show in the next section for our sulfide systems, it becomes smeared out by the phonon-phonon interaction at high temperatures 37,38 .
In this study, we mainly seek to model the reflectivity, which is characterized by a frequency band of high reflectivity, i.e., the reststrahlen band 22 .The top and bottom of this band are defined by the transverse (ω T O ) and longitudinal (ω LO ) optical modes, respectively.These frequencies can be calculated from the harmonic approximation alone by employing the non-analytical correction (NAC) at almost zero q-vector 39 .It also means that the summation in Eq. ( 1) is reduced to only those modes that fall within the reststrahlen band, i.e., the transverse optical (TO) mode.
All our calculations are performed DFT and the projector-augmented plane wave method (PAW) 40 as implemented in the VASP code (version 5.4.4) 41,42 .The valence configurations are 3p 6 4s 2 for Ca, 2p 6 3s 2 for Mg, and 3s 2 3p 4 for S. The exchange-correlation (XC) term in the effective Kohn-Sham potential is approximated according to the Perdew-Burke-Ernzerhof parameterization for solids (PBEsol) 43 of the generalized gradient approximation (GGA) 44 .Both CaS and MgS crystallize into rock-salt (RS) cubic structure.We use their conventional unit-cells with 8 atoms.Integration in the Brillouin zone is done on a Γ-centered grid of uniformly distributed kpoints with a spacing of 2π × 0.3 Å−1 .The selected plane-wave kinetic energy cutoff is 500 eV and convergence of our structural optimizations is assumed when the total energy changes are less than 10 −8 eV and the forces on each atom smaller than 10 −3 eV/ Å.We use fully relaxed crystal structures (in unit-cell shape and ionic degrees of freedom) as the underlying unit-cells for the generation of atomic displacements.
The displacement configurations are generated in accordance with the crystal symmetry as implemented in the Phonopy 45 code.We use 2 × 2 × 2 supercells, which contain 64 atoms, and to obtain second-order FC, two displacement configurations are sufficient.Then, we calculate the harmonic eigenfrequencies with NAC to estimate the size of the reststrahlen band.To compute the third-order FC necessary to calculate the self-energy, i.e., Π B 3ph (ω) 25 , we also use the FDM method as implemented in the Phono3py 46 code.Unlike for the second-order FC, for the third-order case every pair of atoms in the supercell needs to be displaced, but in this instance, the crystal symmetry once again helps reducing the overall number of required configurations.The phonon interaction distance (i.e., the distance between the displaced atoms in the configurations) is chosen to be 7 Å and 8 Å for CaS and MgS, respectively.These parameters allow us to achieve the largest number of displacement patterns possible, that is, 146 unique displacement configurations for both systems.This condition of maximal displacements is necessary to converge the self-energy.Finally, we compute the self-energy on a uniform 12 × 12 × 12 qpoint grid, resulting in the overall number of 720 q-points in the first Brillouin zone.

III. RESULTS: OPTICAL PROPERTIES A. Harmonic approximation limit
Firstly, we analyze the dynamical properties of CaS and MgS within the harmonic approximation.As we discuss above, it is imperative to establish correctly the bounds of the reststrahlen band needed subsequently to calculate the reflectivity.In Table I, we summarize our predicted values for the lattice parameter a, the transverse ω T O and longitudinal ω LO optical modes, the Born effective charge Z, and the high frequency dielectric constant ∞ for CaS and MgS.Their computed phonon band spectra are shown in Fig. 2(a) and Fig. 2(b), respectively.
Our calculated phonon band structure for CaS is in good qualitative agreement with a previous LDA/Pseudopotential/FDM study 20 using the SIESTA implementation 53 as shown in Fig. 2(a).The main quantitative differences are due to lower frequencies predicted in our work for all phonon bands.This incosistency in frequency magnitudes could be partially attributed to our slightly smaller lattice parameter a CaS = 5.633 Å with respect to a CaS = 5.670 Å of Ref. [20], which is slightly in better agreement with single crystal measurements 17 of a CaS = 5.697 Å.However, our phonon band structure seems to be more consistent with experimental data 17,18 , as our predicted LO/TO splitting of 112.86 cm −1 agrees more with reported values of 113 cm −1 by Ref. [17].From Table I, we notice that the two experimental studies find nearly the same value for ω T O , but they differ in the magnitude of ω LO .An apparent reason for this discrepancy could be traced back to the method employed to extract ω T O and ω LO from experiments.Both studies use a single oscillator model in their dispersion analysis, and adopt an expression virtually identical to Eq. ( 2), but they treat the self-energy as frequency independent, giving it the role of a damping constant.In EXP1 17 (in combination with ∞ taken from Ref. [50], the real and imaginary parts of the dielectric function, Re( ) and Im( ), are constructed from reflectivity data, TABLE I. Lattice parameter (a), transverse and longitudinal optical modes (ωT O and ωLO), Born effective charge (Z) and high frequency dielectric constant ( ∞) calculated in this work [*].We compare to other theoretical and experimental data available in the literature: †-Ref.[17], -Ref.[18], ‡-Ref.[47], -Ref.[20], -Ref.[21], •-Ref.[48], -Ref.[49], -Ref.[50].Note that in the Ref. [17], the so-called Szigeti effective charge, ZS, is given instead of Z.The Born effective charge we report here can be obtained using the relationship 51,52 Z = ( ∞ +2) 3 ZS.The high-frequency dielectric constant ∞ used in Ref. [17] and Ref. [20]  was obtained from a semiempirical (SE) model in Ref. [49]  then ω T O and ω LO are taken as the maxima of Re( ) and Inv[Im( )], respectively.On the other hand, in EXP2 18 the TO and LO frequencies are treated as adjustable parameters to parametrize the reflectivity and fit the data.Thus, although the resolution of ω T O seems to be independent of the method, the latter seems to be more in agreement with our computation.Nevertheless, the calculated LO/TO splitting of 70 cm −1 by Ref. [20] is underestimated by at least 38% with respect to both measurements.This underestimation could be linked to their Born charge of Z Ca = 1.802e and high-frequency dielectric constant ∞ = 4.150 taken from Ref. [48] and Ref. [49], respectively.In contrast, our Born charge of Z Ca = 2.350e, calculated directly from our relaxed CaS structure, is much closer to the experimental value of 2.111e.Here, we observe clearly that a self-consistency among the required elements entering in the calculation of ω T O and ω LO is highly desirable, in particular, the LO/TO splitting is rather sensitive to the Born effective charge, as it enters the NAC as Z 2 .An accurate determination of the TO mode is rather significant for the modeling of the reflectivity as it controls the edge of the reststrahlen band 22 .
In the case of MgS, our predicted lattice parameter a MgS = 5.180 Å (Table I) and phonon band spectrum, (Fig. 2(b)) are in excellent agreement with an earlier LDA/Pseudopotential/DFPT study 21 as implemented in the quantum espresso code 54 .Our Born effective charge Z Mg = 2.31e and high-frequency dielectric constant ∞ = 4.150 are also almost identical to the values found in Ref. [21] of Z Mg = 2.314e and ∞ = 4.150.These findings further support the idea that the main reason for the significant difference between our results and Ref. [20] in the case of CaS stemmed from the sizeable difference in the values of Z Ca and ∞ .However, our LO/TO splitting of 152.4 cm −1 , and the other ab initio value of 156 cm −1 are both, about 20% smaller than absorption spectra data 18 reporting a LO/TO splitting of 195 cm −1 .But, as discussed previously for the CaS system, the method used by Ref. [18] to determine ω T O and ω LO seems to give an overestimated value for ω LO (Table I), which could be the main source of divergence between the calculated and the experimental values.

B. Anharmonic effects
We now turn our attention to the anharmonic effects and their role in the optical properties of CaS and MgS.As explained in the previous section, the only contributions to the self-energy that we consider are the threephonon Π 3ph scattering processes and the isotopic disorder Π isot .
We notice that the Π isot contribution is only prominent at very low temperature, and its effect can be seen clearly in the imaginary part of the dielectric function, Im( ), as shown in Fig. 3. Three-phonon processes are essentially absent at low temperature and low wavenumbers, thus the isotopic disorder is the leading term in this region.But as temperature increases, the three-phonon processes become dominant and smear out the isotopedisorder part.In addition, one can also observe the Im( ) peak slightly shifting towards lower wavenumbers.Higher-order scattering processes (e.g., four-phonon contributions) are expected not to have a drastic effect on the computed Im( ) of MgS and CaS.This assumption is in line with modeling results for MgO 38 , as it was demonstrated that the four-phonon processes start being substantial at relatively high wavenumbers when the three-phonon contributions vanish.
Our calculated reflectivity at three different temperatures for CaS and MgS is shown in Fig. 4. Here, we can immediately see the importance of the precise prediction of the ω T O value, as it defines the reflectivity maximum.The width of the plateau at low temperature is mostly determined by the high-frequency dielectric constant ∞ and the Born effective charge Z.As temperature increases, both edges of the plateau smear out, although the change is much more obvious at the higher-wavenumber end.This effect is due to the anharmonic interactions, which become more notorious with increasing temperature.For MgS, Fig. 4(b), we observe characteristic shoulders between ∼325 cm −1 and ∼425 cm −1 , these features appear in MgO 22,38 as well.We also notice that for MgS, these shoulders show a somewhat spiky structure.For CaS, these shoulders are more prominent as temperature increases, that is, for ≥ 300 K.The only available experimental report is on the reflectivity of CaS 17 as shown in Fig. 4(a).If we compare it to our reflectivity computed at 300 K, we observe a good qualitative agreement with some systematic quantitative deviation.Our predicted position of the TO peak is only 0.14% lower than the measured one at 229 cm −1 .In general, we notice that the experimental reflectivity values are smaller than our simulation results in the region considered.It is highly unlikely that the inclusion of additional anharmonic terms in the Cowley's expansion of Eq. ( 2) would reduce the overall magni- tude of our predicted reflectivity.Namely, we think that the main cause for the overestimation of the reflectivity has its roots in the calculated force constants, and consequently, the estimated strength of the anharmonic effects.Perhaps, larger supercells could improve the overall FC values because they would allow for larger interaction distances between phonons, however, the computational cost would increase significantly.The other peculiar feature of our predicted reflectivity curves for CaS and MgS is the presence of several peaks in the high-wavenumber shoulders.Such structures are not observed in the experimental study on CaS 17 , where all possible anharmonic effects are, of course, present.It is conceivable that the inclusion of higher-order phonon processes could remove these spiky structures in that region as it is shown in the case of MgO 38 , in which both three-and four-phonon processes were included and a smooth computed reflectivity is achieved.Therefore, we speculate that at the level of our simulations, both CaS and MgS seem to be more anharmonic in the sense that, more anharmonic terms are needed to be included in order to reach a better agreement with the experimentally observed reflectivity, but such studies are outside of our current resources and are left for future explorations.
At this point, it is worth mentioning that the experimental results used for comparison in our previous discussion are from single crystal measurements.However, there are experiments in which, powdered (pressed pellets) samples are studied also with the infrared spectroscopy technique.This experimental approach is particularly relevant in planetary investigations, as dust is likely to be studied directly on the planet's surface thanks to spectrometers in missions or in simulated conditions in the laboratory 55 .However, the use of pressed pellets leads to scattering, and diffused reflection or refraction, resulting in an overall weaker reflectivity measurement 56 .It has been also observed that, sometimes, reflectivity spectra from powdered samples retain qualitative features of the reflectivity of single crystals, that is, the measured peaks do not shift in wavenumber but only decrease in magnitude. 57.In such scenario, we expect the harmonic modeling using the generalized Lorentz model could still be applicable, but adjusted to measurements from sulfide powders, by modifying the number of oscillators and the damping coefficient.Thus, the self-energy (which plays the role of the damping factor) in Eq. ( 2), would be then dominated not by anharmonic contributions, but by scatterings, grain sizes, etc., i.e., processes specifically differentiating between single crystal and pressed pellets.
In recent years, with the objective of building a spectral library of sulfide minerals to support investigations of Mercury's surface chemistry, a study 55 reported reflectance (R) measurements of synthetic CaS and MgS powders in the range from far infrared (FIR) until visible (VIS), at day Mercury's surface temperature (∼773 K during the day).Reflectance is related to reflectivity in that the latter is the reflectance of a semi-infinite slab 38 .Thus a direct and straight comparison between their reflectance magnitudes is not possible, nonetheless, the position of their maxima, which should coincide, can still be investigated.The measured R maximum in the FIR region for MgS (sample denoted "MgS-2" in Ref. [55] occurs at the wavelength of 40 µm, equivalent to a wavenumber of 250 cm −1 .In comparison, we predict the maximum of the reflectivity at 240 cm −1 and 300 K. Finally, for CaS, the reflectance measurements 55 position its peak at a wavelength of about 8 µm, that is, a wavenumber of about 555 cm −1 , this value is in drastic disagreement with our calculated position for the reflectivity peak at 228.69 cm −1 and the single crystal measurements 17,18 of 229 and 232 cm −1 .This sharp difference in the position of the maxima from reflectance and reflectivity could be due to the use of CaS powder in the case of the reflectance experiment 55 , however, further experiments would be necessary to clarify the discrepancy.

IV. CONCLUSIONS
In summary, we have conducted a full first-principles study of the reflectivity of CaS and MgS in the infrared region.We have considered both harmonic and anharmonic contributions.Within the harmonic limit, we have demonstrated that it is crucially important to employ the non-analytic correction to obtain the correct LO/TO splitting.An accurate determination of the TO and LO peaks is highly desirable as they provide the boundaries of the reststrahlen band, i.e., the low-and highfrequency edges of the maximum reflectivity.We have also shown the effects of the anharmonicity by considering three-phonon scatterings and isotope-disorder processes at the lowest perturbation level.The anharmonic terms' main influence occurred in the smear of the reflectivity spectra's edges, being more prominent in the higher-wavenumber region.
Although we have not included higher-order anharmonic terms, e.g., four-phonon scatterings, we think that these processes would be only noticeable in the highwavenumber region of our study, outside of the maximum reflectivity peak, as it has been shown to be the case for MgO 38 .However, these higher-order anharmonic terms could eliminate the spiky structures in our predicted re-flectivities.Finally, it might be worthwhile to try to incorporate the anharmonic terms into the modeling of the polycrystalline and powder samples.Although, in this case, we expect them to play a secondary role since the self-energy, i.e, the damping constant would be dominated by the processes that are characteristic of powder pellets such as diffused reflection and refraction.

FIG. 1 .
FIG.1.Diagrammatic anharmonic contributions considered in this work.Phonons are indexed by their wavevector q and branch j.Incoming and outgoing phonons are represented by qj and qj , and internal phonon lines are given by qiji respectively25 .Diagrams (a) and (b) represent three-phonon scattering processes known as "bubble" and "tadpole" contributions.Diagram (c) shows a four-phonon scattering process known as "loop".In our study, the only process entering in our model, as explained in the main text, is the bubble contribution.Diagram (d) represents isotope-disorder scattering37 , the dashed lines represent phonon scatterings by the isotope shown here as the black dot.This process is also considered in our anharmonic modeling and it is obtained from second-order perturbation theory, using mass deviation as a perturbation.

FIG. 2 .
FIG. 2. Our simulated phonon bands within the harmonic approximation are shown in solid lines for (a) CaS and (b) MgS.Different colors correspond to individual normal modes of motion: three acoustic bands (red, orange, and blue), and three optical bands (violet, brown, and green are the two TO and one LO mode, respectively).Hollow squares 20 (CaS) and hollow circles (MgS) 21 represent data taken from previous modeling studies.Experimental values shown in solid symbols for the TO and LO modes are denoted by EXP1 18 and EXP2 17 .

FIG. 3 .
FIG. 3. Calculated imaginary part of the dielectric function, Im[ ], for (a) CaS and (b) MgS at three different temperatures.The isotope-disorder scattering effects are responsible for the oscillatory behavior in the region 150-220 cm −1 and are only visible at sufficiently low temperature.The threephonon processes prevail over Πisot at higher temperatures.

FIG. 4 .
FIG. 4. Predicted reflectivity (solid lines) for (a) CaS and (b) MgS at three different temperatures.For CaS, we compare our results to an experimental trend 17 (dashed line) at 300 K.The dotted-dashed gray lines indicate the position of our predicted ωT O modes.