Tuning Einstein Oscillator Frequencies of Cation Rattlers: A Molecular Dynamics Study of the Lattice Thermal Conductivity of CsPbBr$_3$

The pure CsPbBr$_3$ perovskite is an archetypal example of a strongly anharmonic crystal that poses a major challenge for computational methods to describe its thermodynamic properties. Its lattice dynamics exhibits characteristics of a phonon liquid: mode coupling, low lifetimes and `rattlers'. To study the thermal conduction in this crystal, including the effect of dynamic disorder introduced by the Cs rattlers, we apply large-scale molecular dynamics (MD) simulations combined with machine-learning interatomic potentials. We simulate its ultra-low lattice thermal conductivity in the cubic phase and obtain phonon spectra by measuring velocity autocorrelation functions. The thermal conductivity at 500 K is computed to be $0.53\pm0.04\frac{\rm W}{\rm m K}$, which is similar to that of demineralized water under normal indoor conditions. MD based insight into the heat transport mechanism of halide perovskites is presented. In the analysis the Cs cations are interpreted as damped Einstein oscillators. The phonon bandstructure of a system with artificially raised Cs masses demonstrates an increased interference of the Cs rattling with the acoustic phonon modes. We show that the thermal conductivity of the CsPbBr$_3$ perovskite can still be slightly decreased by tuning the cation rattling frequency into the range of the low lying acoustic mode


Introduction
Around room temperature, condensed matter is often represented as a 'wiggling and jiggling' of atoms connected by springs.When this matter displays periodicity, as in a crystal lattice, energy can be stored in collective lattice oscillations.The phonon gas model is a 0 K approximation that assumes that these springs are harmonic.The phonons in the lattice possess infinite lifetimes and can transport kinetic energy unhindered throughout the crystal.In reality, there is always a degree of anharmonicity of the interatomic potential and disorder in the atomic positions.This causes phonon modes to scatter frequently with one another, thereby losing their phase coherence, resulting in finite lifetimes.Even though present day simulations go beyond the harmonic approximation, much is still unknown about thermal transport in crystals that can be classified as strongly anharmonic 1 .The CsPbBr 3 perovskite is a clear example of a strongly anharmonic crystal 2,3 , and forms a challenge for firstprinciples simulation methods 1,[4][5][6][7] .This type of perovskite is technologically interesting because of its opto-electronic effects [8][9][10] and its ultra-low thermal conductivity (κ < 1 W mK ) [10][11][12][13][14] .The material is cheap and easy to produce, because it can be crystallized from its liquid solutions, but it is also toxic and not very stable 10,15 .Its perovskite lattice contains a Cs + 'rattler' 2,16 in the center of the lead-bromide cage, which surrounds the cation as shown in Figure 1.Neighbouring Cs cations show a close to negligible level of dynamic correlation.The effect of rattling has been linked to the low thermal conductivity of several (types of) halide perovskites 2,[16][17][18][19][20] , possibly inspired by resemblances with rattlers observed in hostguest complexes: filled skutterudites 21,22 or clathrates 23 .Interest in these materials results from efforts to find a phonon-glass electroncrystal 24 , which would make an ideal thermoelectric.The effect of the guest rattlers is often described by localised, independent Einstein oscillators 21 .In this scenario these loosely bound, often heavy, atoms produce strong incoherent scattering that impedes thermal transport as illustrated in Ref. 25 and by molecular dynamics (MD) simulations in Ref. 26 .However, this explanation has been challenged by experimental work indicating that there is a strong coupling between guest and host modes [27][28][29] .Furthermore, simulations have indicated that it is the phonon lifetimes (over a broad spectral range) that are affected by the rattlers rather than their group velocities 30,31 .These differences arise because of uncertainties in the coupling strength between the host's phonons and the rattler phonons.These uncertainties are also present in the halide perovskites, which raises the question how rattling affects CsPbBr 3 .In this work we simulate the lattice thermal transport of this perovskite in its cubic phase using an on-the-fly trained machine learning (ML) potential 6,32 .We study the influence of rattling by a computational experiment in which we tune the rattling frequency and the level of disorder of the Cs atoms.
The challenge in modeling thermal transport in strongly anharmonic crystals stems from the complexity of their potential energy surface (PES).As a result a highly interacting phonon gas, or maybe even phonon liquid, is formed in the cubic phase of CsPbBr 3 32 .Most of the phonon quasi-particles are ill-behaved; with non-Lorentzian multipeak spectra and/or lifetimes smaller than the inverse of their frequency (Ioffe-Regel limit).In CsPbBr 3 these two criteria are not satisfied, since the phonon spectra show mode coupling 1,32 .This implies that lattice thermal transport may not be accessible with models (κ L ) relying on the Peierls-Boltzmann (PB) transport equation in the relaxation time approximation, PB relies on decoupled and well separated phonon branches with frequencies ω α (q) and group velocity v α (q) = ∇ q ω α (q).Each phonon mode α has a relaxation time τ α , a population n and contributes C α (q) = (hωα(q)) 2 k B T 2 n(ω α (q)) (n(ω α (q)) + 1) to the lattice heat capacity.
The phonon dispersion ω α (q) can be readily obtained from standard ab-initio calculations, by computing the second order (harmonic) force constants.However, higher order potential derivatives are required to obtain fi-  nite lifetimes [33][34][35][36] .These calculations are exceedingly expensive and can quickly become computationally intractable for complex materials, but have been done for CsPbBr 3 including effects of up to 3rd and/or 4th order derivatives 1,37,38 .
The work of Simoncelli et al. 1 goes beyond PB transport and should be better at describing thermal transport in anharmonic lattices.To summarize, they approximate the lattice thermal conductivity as the sum of two contributions, This includes a Peierls (populations) term corresponding to intraband transport (κ P in Eq. ( 1)) and a coherences term corresponding to inter band transport (κ C ).These contributions scale differently with the lifetime; with decreasing τ , κ P reduces, whereas κ C increases.The latter describes diffusive heat transport carried by diffusons as identified by Allen and Feldman 39 .Applying the thermal conductivity equation of Simoncelli et al. to the case of low temperature (50-300 K), orthorhombic CsPbBr 3 shows a temperature dependent thermal conductivity κ L (T ) in agreement with experiment 1 .Recently, Tadano et al. 37 have applied the same formalism, but used renormalized frequencies {Ω α (q)} originating from a GW -type approach that describes the broadening of the quasi-particle spectrum by phonon-phonon interactions.They find, in the high temperature (400-800 K) cubic phase, a small temperature dependence with a maximum κ L (∼ 450K) = 0.51 W mK , and κ L (500K) = 0.50 W mK .Supporting findings have been very recently reported by Wang et al. 38 .However, depending on the level of self-consistent phonon theory to calculate {Ω α (q), τ α (q)}, differences in κ L of up to 50% are reported.
We make use of recent developments in ML interatomic potentials [40][41][42][43] to simulate heat transport using the full PES and (in principle) including all higher order scattering processes.These developments have opened up a computationally efficient description of the interatomic potential.Using MD and this potential we simulate large thermodynamic ensembles of cubic CsPbBr 3 at 500 K with near-first principles accuracy.Properties related to ion dynamics, such as thermal conductivity or phonon dispersion are extracted by a 'measurement' of the temperature gradient or the projected velocity auto-correlation functions.This is a purely classical approach, valid above the Debye temperature, that naturally includes the effect of dynamical disorder that is absent or approximated in the before mentioned lattice thermal conductivity models.Disorder is expected to drive the system further away from the harmonic regime.The heat transport coefficient (κ) is calculated by thermal conductivity measurements in a non-equilibrium MD (NEMD) setup [44][45][46][47][48] .As with the lattice thermal conductivity models, this method has its own convergence conditions.Using the trained ML potential, we here observe a relatively fast convergence of κ with MD simulation size/length.However, we note that the required simulation size is too large for a fully ab-initio MD to be computationally feasible.In order to increase our phenomenological understanding of heat transport in this system, we compare to the predictions made by the lattice models 1,37,38 .we probe phonon frequencies, lifetimes and populations under full phonon-phonon scattering conditions in the microcanonical ensemble.Even though the interband contribution to the thermal transport may be substantial, we show that the remaining Peierls (intraband) contribution is affected by Cs rattling.Therefore there is room for changing the thermal conductivity by tuning the effective rattling frequency.

Computational experiment
In Figure 2 the harmonic phonon dispersion 49,50 ω α (q) of cubic CsPbBr 3 is shown.By highlighting (in red color) the Cs contribution in the respective eigenmode |e Cs,q (ω)| 2 |eq(ω)| 2 , we identify nearly dispersionless 'rattling bands'.In the middle panel we see that the Cs cations mainly contribute to the phonon bandstructure in an energy window around ∼0.7 THz.In the left/right panel we have artificially lowered/raised the masses of only the Cs atoms, which results in a blue/red shift of the central rattling frequency, respectively.For the light Cs cations the central frequency shifts up to ∼ 2.2 THz, thereby almost removing the overlap with the linear acoustic bands.Conversely, the rattling band of the heavy Cs cations now overlaps at ∼ 0.2 THz with the acoustic bands.Here, we will explore the effects of this rattling band on the lattice thermal conductivity of CsPbBr 3 .In the analysis we view the Cs rattlers as Einstein oscillators with a characteristic frequency (ω E ) and a Lorentzian distribution in the vibrational density of states (VDOS) as plotted in Figure 1.We perform MD simulations for different Cs masses (M ) to study the influence of rattlers on the phonon spectrum and the thermal conductivity.Altering the mass does not change the PES as described by Density Functional Theory (DFT), but modifies the Einstein frequency as ω E = k E M , with k E the spring constant.Furthermore, it affects the level of disorder exhibited by the Cs atoms.The same perturbation in the force has a larger effect on the displacement of a light compared to a heavy Cs atom.Only by performing the MD simulation can we find out how these two effects affect κ.Therefore, for the systems with a low, normal and high Cs mass, we apply two types of simulations: i) a passive measurement of the phonon spectrum using VACFs, and ii ) an active measurement of the thermal conductivity.

Computational Method
We start by briefly summarizing the applied computational approach as a step-by-step plan: • Train MLFFs for CsPbBr 3 with light, normal and heavy Cs masses.
Then, for each of the three different mass systems we: • Simulate N V E ensembles in a 10×10×10 cubic supercells of CsPbBr 3 at 500 K.
• Calculate q-VACFs in post-processing of the N V E-trajectories for all phonon modes.Calculate PVACFs for the first acoustic mode and fit Lorentzians to the power spectrum.
• Simulate heat flow in a NEMD setup of 4 × 4 × 20 supercells with leads set at 500±25 K. Five independent trajectories starting from snapshots of the N P T ensemble are calculated.
• The ensemble average temperature gradient and heat flux are computed in postprocessing and result in the thermal conductivity.
We now expand on some of the most important computational details of these steps.More details have added to the supplementary information (SI).
We simulate the heat flow through the lattice using large-scale MD and the Machine-Learning Force-Fields (MLFF) method 6,51,52 as implemented in the DFT code vasp 53,54 version 6.3.The MLFF training procedure and its accuracy benchmarking are described in the SI.Here, we want to mention that the MLFFs were trained with respect to the PES as described by DFT using the SCAN 55 exchange-correlation functional, which has been shown to be accurate for this type of perovskite 56,57 .
The q-resolved VACFs 32,58-60 are computed from microcanonical ensembles.The velocity field of atoms type s in reciprocal space is given by Here v s (n, t) is the velocity at position r s (n, t) at time t inside the unit cell denoted by n with mass m s .We compute the self-correlation and perform a temporal Fourier transform, i.e.
where ω denotes the temporal frequency axis and t = t ′ − t ′′ .By taking the power spectrum of this quantity the q-resolved form of the phonon density of states (DOS) is obtained.The projected VACF is obtained by projecting v s (n, t) along the harmonic phonon eigenvectors e s,α (q), summing over the species index s, and computing the transform as in Eq. ( 4).The resulting PVACF H α (q, ω) is then decomposed into 15 phonon branches α of the CsPbBr 3 cubic phase.The (P)VACFs were computed with the open-source code dsleap 32,61 .
To measure the thermal conductivity in NEMD, we apply the approach of a constant temperature gradient ( ∆T L ) [46][47][48] , imposed by the Langevin thermostat 62 with a low friction coefficient of 1 ps −1 .We observe that this approach converges the thermal conductivity value within acceptable simulation times, see SI.The geometry is split into three regions: a hot region (Γ H ), a cold region (Γ C ) and a scattering region of length L where no thermostating is applied.The desired temperatures for Γ H and Γ C are pinned to T H and T C by the thermostat.As a result, an energy flux in opposite direction of the temperature gradient is induced.Here, E hot and E cold are the amounts of kinetic energy the thermostats insert into and extract from the simulation cell in one time step ∆t, respectively, and A is the area perpendicular to the heat flow.To compute the thermal conductivity Fourier's law is used; where ⟨.⟩ denote ensemble averages.The factor two arises due to the applied periodic boundary conditions, resulting in two scattering regions in each simulation cell (see Fig. 5).Therefore, mirror symmetry is applied in determining the average temperature profile.

Results
We start with a characterization of the anharmonic vibrations in the three CsPbBr 3 model systems as they appear in MD. q-resolved VACFs were computed from 10 microcanonical MD runs at 500 K in 10 × 10 × 10 cubic supercells.The averaged power spectrum of Eq. ( 4) is shown in Figure 3 and shows the Cs and Pb contribution to the phonon spectrum.
As in the harmonic approximation, decreasing the mass in the MD simulation shifts the Cs rattling band up, as can be seen in Fig. 3(a).The shift is somewhat larger than in Fig. 2. In Fig. 3(c) we see that the rattling frequency is shifted down in the case of heavy Cs atoms.We are interested in their impact on Pb vibrations in the broad rattling window indicated by the red rectangles in Figs.3(d-f).We will get back to this point later on.The rattling Cs atoms introduce dynamical disorder into the system, which we measure by the displacement d Cs (t) of the Cs atom from the geometric center of its enclosing PbBr host lattice.We extract effective rattling frequencies ω E and correlation times from the selfcorrelation functions ⟨C Cs (t)⟩ = ⟨d Cs (t)d Cs (0)⟩ as shown in Figure 4(a).The normal and heavy mass (blue and green) behave as damped Einstein oscillators, the light mass shows overdamped dynamics.All have very short selfcorrelation times, ranging from sub-picoseconds (low mass) to several picoseconds (heavy mass).The power spectra of the Cs motion shown in Fig. 4(b) show clear peaks in the case of normal and heavy Cs masses that are fitted with Lorentzians.These peaks are located at ω E = 0.8 THz and 0.3 THz, respectively.The lighter systems show broader Lorentzians, indicative of more dynamic disorder.The spectrum of the light mass system is so broad, that it more closely resembles white noise than a Lorentzian.In this case the rattling picture does not apply anymore.
The lattice thermal conductivity of CsPbBr 3 in its cubic phase at 500 K is determined by com-   , where L is the distance between the leads, and T H = 525 and T C = 475 K are the target temperatures of the leads.These temperatures keep all parts of the simulation cell in the cubic crystal phase.By applying weakly coupled Langevin thermostats to the leads, a temperature profile without interface effects [46][47][48] or interface thermal resistances at the lead|scattering-region interface is obtained, see SI.
The computational efficiency of the MLFF method allows the simulation of atomic scale heat flow in large 4 × 4 × 20 supercells of the cubic unit cell.To converge to the bulk conductivity value L should be (much) larger than the phonon mean free path (Λ).In this work, we assume that the length of the scattering region (L = 8a lat ≈ 46 Å) is sufficient, considering the already very low Λ values reported for CsPbBr 3 at 350 K in Ref. 64 .See the SI for a comparison to L = 3a lat simulations, which result in comparable κ values.
The obtained thermal conductivity values for the different independent simulations are tabulated in Table 1.We determine the value of the thermal conductivity κ for CsPbBr 3 by the mean of the five independent simulations to be 0.53 W mK , with a 95% uncertainty interval of ±0.044.This is in agreement with experimental values 10,16,65 and slightly above the modeled thermal conductivity coefficient of Simoncelli et al. 1 (who reported a value of ∼ 0.4 W mK ) and within the spread of the three approaches by Tadano et.al. 37.Our simulations indicate that decreasing the mass of the Cs atoms (blue shifting ω E ) by an order of magnitude only changes the thermal conductivity slightly, resulting in 0.58 ± 0.049 W mK .This is a small but statistically significant increase in κ.Our simulations also show that increasing the mass of the Cs atoms (red shifting ω E ) by an order of magnitude changes the thermal conductivity more substantially, resulting in ∼ 0.36 W mK .These observations combined show that in this crystal that exhibits ultra-low κ, merely tuning the Einstein rattling frequency still has a noticeable effect on the thermal conductivity.
We will now analyze whether the high/low mass system enhances/decreases scattering with low frequency acoustic modes, respectively.
From a Peierls-Boltzmann (PB)-transport point of view this would either decrease or enhance κ.We will first establish whether the acoustic modes are affected, by comparing the acoustic branches of the Pb VACF spectra in Figs.3(d-f) close to Γ, we observe a reduced intensity of the acoustic bands in Fig. 3(f).This perturbation of the acoustic modes occurs slightly off the Γ point precisely around ω E =∼ 0.3 THz, whereas the higher frequencies are largely unchanged.A similar effect is observed in the Br spectra, see SI.The coupling of the heavy mass rattling band with the acoustic mode can be more clearly identified in Figure 6(a).In the power spectra of the first acoustic phonon branch α = 1 (black lines), we observe Lorentzian line-shapes for both the light and the normal Cs mass along the whole Γ-X line in q-space.This indicates that α = 1 is a well-defined phonon-quasi-particle, albeit increasingly short lived as indicated by their lifetimes shown in Fig. 6(b).For the heavy mass system we observe multi-peak spectra at q = [ 1 5 00] and [ 3 10 00].Around ∼ 0.7 THz a second peak appears that is much wider.This indicates that the acoustic mode is coupled with low lying optical modes, which makes the Lorentzian lineshape approximation not applicable here and it cannot be uniquely fitted.Spectra of the other α = 2, 3 acoustic branches are in the SI and show similar behavior.Having established this, we fit Lorentzians (colored lines) to the highest peak in each of the 15 spectra.We extract all required parameters for the PB conductivity of Eq. ( 1): the central frequency ω of the peak, the occupancy n from the peak height, and the phonon lifetime τ from the full width at half maximum (Γ) as τ = 1 2Γ .The resulting frequencies and lifetimes are shown in Fig. 6(b) and the relative peak heights can be read off Fig. 6(a), taking into account the indicated magnification factors.The dispersion of the acoustic mode as observed in the normal and light mass system are identical and lie below the rattling frequency ω E (blue dashed line) of the normal mass system.However, ω E of the heavy mass system (green dashed line) intersects with the acoustic band.As a result, the dispersion splits, showing a red shift close to Γ, which turns into a blue shift upon approach- ing X.The lifetimes of the acoustic phonons rapidly decrease with the distance away from q = Γ.
These results are combined in the calculated PB conductivity of this mode κ α=1 P (q) as shown in Figure 7. Details of the computation have been added to the SI.The PB conductivity at low q is significantly reduced for the heavy mass compared to the normal mass system, whereas it is of the same order when approaching the X-point.This is mainly due to the mode's low occupancy and further reduced by the lower group velocity.The κ α=1 P of the low mass system shows comparable behavior to the normal mass system.

Discussion
The ML MD approach applied here presents an ultra-low lattice thermal conductivity of CsPbBr 3 at 500 K that is slightly higher than the available experimental values at room temperature.This is expected, since the simulation covers a single crystalline domain without impurities, defects and grain boundaries.The presented lifetimes and lattice thermal conductivities are of the same order as those in recent lattice thermal conductivity model works of Refs. 37,38,64.The model of Simoncelli et.al.
shows that only ∼ 30% of the thermal conductivity of orthorhombic CsPbBr 3 at 300 K originates from κ P , and most of it is linked to the acoustic modes 1 .The remaining ∼ 70% originates from κ C .The model of Wang et.al.supports this assessment, stating that in cubic CsPbBr 3 optical phonons, ranging from 3.1 to 3.9 THz, dominate the heat transport 38 .On the contrary, Tadano et.al.report that ∼ 90% of the thermal conductivity of cubic CsPbBr 3 at 500 K originates from κ P , dominated by modes below ∼ 1.5 THz 37 .It should be noted that these works apply different levels of theory to obtain the phonon basis.
The coherences' contribution κ C was not calculated, since the strongly interacting phonons observed in the MD make a decomposition into projected phonon frequencies and lifetimes far from unique.The fingerprints of these strong interactions are evidenced by the non-Lorentzian lineshapes that are found (apart for some low frequency acoustic modes) across most of the spectrum 32 .Therefore, with the current NEMD study we can not make a quantitative statement about the relative importance of κ C vs. κ P .However, the observation made here that the thermal conductivity can still be significantly tuned, and that this can be explained by combining the rattling picture with a PB-view on thermal transport, seems to advocate for a κ P contribution of importance.We note that we cannot rule out a back-coupling effect, where the decrease of κ P affects the value of κ C .We would therefore highly encourage a study with lattice models making the same alterations to the Cs masses as in the present study.In future work, a modal decomposition of the thermal conductance 66 or conductivity 67 as computed by MD could help in resolving this question.

Conclusion
A near first-principles accuracy ML interatomic potential has been used to simulate large thermodynamic ensembles of the cubic phase of CsPbBr 3 .The Cs cations introduce dynamic disorder into the crystal, and, depending on their mass, rattle at different effective frequencies.We have shown that the ultra-low lattice thermal conductivity can be tuned by creating a substantial change in the Einstein oscillator (rattling) frequency.The lattice thermal conductivity was computed to be 0.53 ± 0.04 W mK for normal Cs mass and reduced to 0.36 ± 0.02 W mK for a Cs mass ten times heavier.By lowering the mass to a tenth of normal Cs we shift the rattling band far up in the spectrum.However, only a slight increase of the thermal conductivity to 0.58 ± 0.05 W mK is observed, which indicates that the main cause for the ultra-low thermal conductivity is not incoherent scattering with guest Einstein modes.Rather, the high degree of mode coupling across the spectrum results in a phonon liquid in which diffusive heat transport dominates.These results provide new insights applicable to the wider class of halide perovskites.For practical thermoelectric engineering purposes we have shown that reducing κ by a suitably chosen cation based on its rattling frequency in the Pb-Br framework is possible, but has only a limited effect.However, this does not exclude chemical possibilities, such as cation replacement Cs→ X, whereby X is non-spherical or has different volume/internal structure 11,68,69 , introducing multiple species X ′ , and mixing the metal and halide compositions.

Data Availability Statement
Research data for this paper has been made available through a data set in the 4TU.ResearchData repository; see Ref. 70 .The following data is stored: (i) The electronic structure databases, including crystal structures and corresponding DFT energies, forces, and stress tensors, used to train the MLFFs (ML AB files).(ii) A high level analysis of the electronic structure databases presented by pdf fact sheets generated with open-source fpdataviewer software 71 .(iii) The generated ML force fields (ML FF files) applied in the NEMD simulations in this study.the q-VACF spectra, results of test calculation to determine the computational setup of the active measurements, additional NEMD simulation results for a small 4×4×10 supercell, explanation of the applied Lorentzian fitting procedure of the PVACF spectra, additional PVACF spectra for the three acoustic modes, and an explanation of the Peierls' conductivity calculation.
structure is heated during an on-the fly training run starting at 50 K to 500 K in the NPT ensemble.A time step of 5 fs was used and 40,000 steps were simulated in total, resulting in a temperature gradient of 2.25 K/ps.The external stress was set to 0.001 kbar.The cutoff radii for the radial and angular descriptors were set to 5 Å.The number of spherical Bessel functions for the radial part and the angular part of the descriptor were set to 8. The maximum angular momentum number for the expansion of the angular part was adjusted to l max = 4.After the training during heating was done, another training round is done at 500 K.For the high temperature training we took the last structure of the heating run and propagated this structure for 100 ps.During training we were using Bayes linear regression to determine the regression coefficients.Before the production runs we recomputed the regression coefficients of our force fields with singular value decomposition (SVD).This procedure was reported to give lower root mean square errors (RMSE) in the energy forces and stress tensor 7 .The number of used reference configurations for the different atoms can be found in table S1

Error analysis
To quantify the error in the energy, the forces and the stress tensor of the MLFF a set of independent test structures was created.To create a set of test structures the CsPbr 3 force field of Refs. 1,8was used.With this force field test structures were created during a heating run, starting from 50 K and going up to 500 K.The temperature gradient was adjusted as in the training procedure.From this heating run 200 equidistant structure files were extracted, giving a POSCAR file after every 1 ps or every 2.25 K.For this 200 structures the energy, forces and stress tensor were computed with the SCAN functional and the force-fields.The average energy error, was computed according to an absolute energy difference divided by the number of atoms.
The error in the forces is computed as a RMSE, x,y,z α=1 where the index α is running over the Cartesian components of the force vector and the index i runs over all the atoms in the supercell.The force vector is plotted against the average magnitude of the DFT force of the considered structure.The RMSE in the stress tensor T is computed over the 6 independent components of the tensor, The errors for the different Cs mass force-fields are summarized in Fig. S1 for the light mass, in Fig. S2 for the normal mass and the error of the heavy mass is shown in

q-VACF
The dynamic properties of the crystal are accessed with the q-resolved velocity-autocorrelationfunction (q-VACF).The q-point is one of the commensurate reciprocal grid points.In this manner we are able to qualitatively determine how a change of the Cs mass affect the phonon band structure.We compute the q-VACF by Fourier transforming the mass weighed velocity field to q-space ṽs (q, t) Then Eq. ( 4) is self-correlated and the temporal Fourier transform is computed gs (q, ω) = ṽs (q, t ′ )ṽ s (−q, t ′′ )e −iωt dt, with t = t ′ − t ′′ and dt the corresponding differential.The power spectrum of the so obtained function is a q-resolved form of the phonon DOS 8 .How the power-spectrum and the ensemble average of Eq. ( 5) are computed is described in detail in Ref  To compute the thermal conductivity we used the method of active measurement with a constant temperature gradient [9][10][11] .A finite size analysis was carried out by simulating an NVT ensemble.The resulting temperature profiles after averaging the temperature per layer over time are shown in Figure S7.Artefacts appear in the temperature profile near the thermostated leads.Whereas the jump up in temperature at the cold lead (position left) could still be interpreted as an interface thermal resistance, the jump down at the hot lead (position right) makes no physical sense.In an attempt to resolve this artefact we stay closer to the linear response regime, in which the leads cause a small perturbation out of equilibrium.Therefore, the temperature difference between the leads was lowered to 500 ± 25 K, followed by lowering the Langevin friction coefficient to 1 ps −1 .This removes the unphysical temperature jumps between the leads and the scattering region.Knowing this, we attribute the artefact to the MLFF (as described in section SI 1.1) being inferred for local atomic environments that lie too far outside the phase space on which it was trained.It was trained at 500 K with a friction coefficient of 10 ps −1 , however the hot lead was set to 550 K, and the same friction was applied.By lowering both ∆T and the Langevin friction, the MLFF remains able to assign a correct local energy to the hot From this point on, and in all simulations presented in the main text, a Langevin friction coefficient of 1 ps −1 was used.We have observed that the variance in the κ values from independent simulations in the same dimension supercell is more important, than the effect of a changing its size (Fig. S7).Since the expected mean free path of almost all phonon modes at 500 K is below 10 Å (see discussion in main text), a simulation box with a scattering region length of L = 3a =∼ 18 Å (region in which no thermostat is applied) can start to approach the bulk thermal conductivity.For this small supercell the temperature profile that develops on average over 1 ns of MD simulation is shown in Figure S8(left).In the energy balance of Fig. S8(right) we see that after ∼ 200 ps the NEMD simulation reaches a steady state.On average, the hot lead pumps in as much kinetic energy as the cold lead extracts from the system.In Table S2 we have compared the κ values obtained from 5 independent simulations of a 4 × 4 × 10 supercell (L = 3a), to the values presented in the main text for a 4 × 4 × 20 supercell (L = 8a).What we see is that the mean value is practically the same, but the confidence interval is higher in the smaller supercell.
In this way we approximate the spectrum by the best-fitting single broadened oscillator. ) for the first acoustic mode (α = 1).The solid/dashed lines correspond to optimizing w.r.t.L 1 /L 2norm, respectively.
We have tested whether L 2 -norm affects the fitted {Γ, ω 0 }-values.Figure S10 shows that the L 2 -norm results in the same frequency (ω 0 ) and, the same lifetimes ( 1 2Γ ).Only at q = [1/10 0 0] the coarseness of the frequency grid combined with the narrow peak width S12 makes that Γ depends on the applied norm.In case of the heavy and the normal mass system, the L 2 -norm results in a more narrow Lorentzian and therefore a higher lifetime (heavy mass: 1000 ps, normal mass: 303 ps).
For comparison with Figure 6 in the manuscript, we include in Figure S11 the PVACF spectra projected onto the eigenvector of the first, second and third acoustic mode (α = 1, 2, 3).It shows that a similar plot as Figure 6(a) in the manuscript can be made for the α = 2 mode.The multi-peak spectra observed for the α = 3 mode become very broad.

Calculation of Peierls' conductivity
The extracted phonon quasi-particle properties for the first acoustic mode (α = 1) based on Lorentzian fitting as explained in the previous section are shown in Table S3.The group velocity v α (q) = ∇ q ω α (q) is computed by numerical differentiation.The central differences scheme is applied on the interior points on the line Γ − X, and the Euler forward/backward differentiation is applied at the start/end points, respectively.The occupancy of the mode n is determined by the height of the Lorentzian and is expressed as a relative number.The same holds for the group velocity, ie.prefactors are not taken into account.The absolute value of the resulting Peierls' conductivity (Eq.(1) in the main text) has therefore an arbitrary prefactor, but ratio's between the 15 numbers in Table S3 are correct.To indicate this, the κ P values in Figure 7 of the main text have been expressed in percentages.

Figure 1 :
Figure 1: Cs rattlers as Einstein oscillators in CsPbBr 3 .(a) The vibrational density of states (VDOS) of CsPbBr 3 at 500 K decomposed in a Pb+Br and a Cs contribution.A Lorenzian distribution (dashed lines) centered at frequency ω E approximates the Cs contribution of the VDOS.Changing the Cs mass can blue or red shift ω E .(b) Opposite to Pb and Br, Cs atoms are approximated as not interconnected, damped oscillators in the PbBr 3 framework.

Figure 2 :
Figure 2: Shifting the 'rattling band'.Harmonic phonon bandstructure of cubic CsPbBr 3 with light ( m Cs 10 ), normal (m Cs ) and heavy (10m Cs ) effective mass of the Cs cation.The colorscale from black to red indicates the Cs contribution to the respective eigenmode.

Figure 3 :
Figure 3: q-VACF of Cs and Pb atoms for light, medium and heavy Cs masses.The frequency domain of the Cs rattling band is indicated by the red rectangles in (d-f).

10 M
Figure 4: (a) Averaged Cs displacement selfcorrelation function and its (b) power spectrum for light, normal and heavy Cs mass.The displacement vector connects the Cs atom and the center of the surrounding PbBr cage.

Figure 5 :
Figure 5: Temperature profile obtained from five independent NEMD simulations of CsPbBr 3 (normal mass).Thermostats keep the (red) Γ H -layers at high and the (blue) Γ C -layers at low constant temperature.A stable temperature profile under periodic boundary conditions develops by time averaging the layer-resolved temperatures.The orange line is the fit determining ⟨∇T ⟩.

Figure 6 :
Figure 6: Projection of the first acoustic mode (α = 1) on the MD trajectory.a) q-PVACF spectra for q-points on the Γ−X line.The spectra (black lines) for the low, normal and high Cs mass are shown from top to bottom.In some of the subplots the intensity of the spectra have been multiplied by the indicated factors.A Lorentzian function (color lines) is fitted to each spectrum.b) The frequencies and lifetimes of the fitted Lorentzians.

Figure 7 :
Figure 7: Peierls-Boltzmann thermal conductivity based on the Lorentzian fits of PVACFs of the first acoustic mode (α = 1).

FigFigure S1 :Figure S2 :Figure S3 :
Figure S1: Error analysis of MLFF trained with m Cs = 13.29054.a: Error in total energy per atom.b: RMSE in of the force.c: RMSE of the stress tensor.The x-axis is the temperature at which the test structure was taken.

8 .Figure S4 :
Figure S4: Pb contribution to the power-spectrum of the q-VACF for varying Cs mass.The mass increases from left to right.

Figure S5 :Figure S6 :
Figure S5: Br contribution to the power-spectrum of the q-VACF for varying Cs mass.The mass increases from left to right

1 ns long single trajectories in 2 × 2 × 2 K
FigureS7: Convergence of the heat gradient for various system sizes in the active measurement approach using a too large Langevin friction (10 ps −1 ).The system sizes and the calculated κ are shown in the plots.

Figure S8 :
Figure S8: CsPbBr 3 (normal mass) NEMD simulations of a short 4 × 4 × 10 supercell under periodic boundary conditions.Left: Temperature profile, thermostats keep the (red) Γ Hlayers at high and the (blue) Γ C -layers at low constant temperature.A stable temperature gradient develops in the time average of the layer temperatures, the orange line is the fit determining ⟨∇T ⟩.Right: Energy balance during the NEMD trajectory visualised by a moving average with window size of 50 ps.In both the left and right figure the lines represent data averaged over five independent trajectories.

FigureFigure S9 :Figure 6 6 )
Figure S9 for the individual trajectories.The ensemble averages in Fourier's law, κ = −⟨∆Q⟩ 2⟨∇T ⟩ , are evaluated for an increasing part of the in total near 1 ns trajectories.In the left figure we convergence of κ in the L = 3a (orange lines) is compared to L = 8a (blue lines).The values in Table S2 are the κ values at the ends of the curves.In the right figure the mass of the Cs cations was varied.Here, all trajectories are for 4 × 4 × 20 supercells.What we see is that κ values obtained from simulations starting from snapshots of the NPT ensemble converge within a uncertainty interval.This interval is narrow enough such that we can conclude that on average κ light > κ normal > κ heavy .

Table 1 :
Lattice thermal conductivity (κ) in W mK for CsPbBr 3 at 500 K with different Cs masses.Confidence intervals at 95% uncertainty are computed with the Student's t test 63 .

Table S1 :
Number of reference configurations in the various force fields

Table S2 :
Lattice thermal conductivity (κ) in W mK for CsPbBr 3 at 500 K with normal Cs masses.Results form a 4 × 4 × 10 and a 4 × 4 × 20 supercell, with length of the scattering region of 3a and 8a, respectively.