Temperature- and vacancy-concentration-dependence of heat transport in Li$_3$ClO from multi-method numerical simulations

Despite governing heat management in any realistic device, the microscopic mechanisms of heat transport in all-solid-state electrolytes are poorly known: existing calculations, all based on simplistic semi-empirical models, are unreliable for superionic conductors and largely overestimate their thermal conductivity. In this work, we deploy a combination of state-of-the-art methods to calculate the thermal conductivity of a prototypical Li-ion conductor, the Li$_3$ClO antiperovskite. By leveraging ab initio, machine-learning, and force-field descriptions of inter-atomic forces, we are able to reveal the massive role of anharmonic interactions and diffusive defects on the thermal conductivity and its temperature dependence, and to eventually embed their effects into a simple rationale which is likely applicable to a wide class of ionic conductors.


I. INTRODUCTION
The lithium-rich antiperovskite solid-state electrolyte (SSE) Li 3 ClO has emerged in the last decade as a promising candidate for all-solid-state lithium-metal batteries: it is superionic at room temperature with a large ionic conductivity; it is environmental friendly and made of light and cheap elements; it is not flammable and has demonstrated a good cyclability. 1 Furthermore, its wide electronic band gap leads to a very low electronic conductivity and a large electro-chemical stability window. 2 Finally, it is also chemically stable against Li-metal formation, which would negatively affect the battery performance via dendritic short-circuits.
A proper account of heat dissipation is key in the design and actual production of batteries: in fact, an excessively low thermal conductivity may lead to overheating, especially during fast charging cycles, which may itself prompt catastrophic incidents, such as melting or explosion. In light of this, thermal runaway can be rightly considered "the key scientific problem in battery safety research" (verbatim from Ref. 3). In addition, thermal dissipation governs energy saving and scavenging: a compromise must be established between minimising heat losses while maximising electric flow during the charging cycle.
A relatively small number of experimental studies on heat dissipation in fast ionic conductors have been reported, as compared to the extensive literature on electric transport. Measurement were performed on superionic materials such as α-LiIO 3 and Li 2 B 4 O 7 ; 4,5 sulfates * federico.grasselli@epfl.ch Li 2 SO 4 and Ag 2 SO 4 ; 6 the quasi-1D material LiCuVO 4 ; 7 lithium aluminium germanium phosphate glass-ceramics compounds of the form Li 1+x Al x Ge 2−x (PO 4 ) 3 ; 8 yttriumstabilised lithium zirconate phosphates, of the form Li 1+x+y Y x Zr 2−x (PO 4 ) 3 . 9 Of these, only the last two classes of materials can can be directly used as SSEs. From the theoretical standpoint, even fewer works are to be found. Even though some models that account for the contribution of the diffusing ions to thermal transport were introduced in the past decades, 10,11 a thorough study of heat dissipation in Li 3 ClO-and in SSE in general-is still missing. At the best of our knowledge, there exists only one calculation of the thermal conductivity, κ, of Li 3 ClO, reporting κ = 22.49 W m −1 K −1 at ambient temperature, i.e., more than one order of magnitude larger than the standard value found in ceramic SSEs. 8,9 Nonetheless, this seemingly promising result is obtained via a rather crude approximation to the Peierls-Boltzmann transport equation (BTE), namely the Slack model, known since its development to have a satisfactory agreement with the experiments (i.e., within ±20%) only for exceedingly simple materials, such as the rare-gas solids, while it is in general poorer for other systems. 12 Furthermore, the Slack approximation totally neglects the effects of defects/vacancies and, just like any BTEbased model, it cannot handle the spurious contributions to heat transport induced by the diffusion of Li ions. These limitations, stemming from a phonon(/normalmode)-based approach to heat transport, can be naturally and easily bypassed through the Green-Kubo (GK) theory of linear response, [13][14][15][16][17][18] which holds for solids with anharmonic interactions of any strength, as well as for diffusing systems, such as liquids and superionic solids. Recently, the GK theory of thermal transport has been combined with state-of-the art quantum simulation methods based on density-functional theory (DFT) 17,19,20 and successfully employed to evaluate κ from equilibrium molecular dynamics simulations of superionic systems, even at extreme pressure and temperature conditions. 21,22 In this paper we first address the structural and mechanical properties of Li 3 ClO, along with their temperature dependence, by means of state-of-the-art ab initio calculations based on DFT. Leveraging these results, we then explore the thermal transport of Li 3 ClO by calculating its thermal conductivity within i) the Slack model; ii) a state-of-the art implementation of the BTE; and iii) the GK theory of linear response and classical molecular dynamics (MD) simulations using both classical force fields and machine-learnt interatomic potentials, trained on DFT data. Our results show that the presence of LiCl divacancies in non-stoichiometric systems, though increasing the Li-ion diffusivity, strongly reduces the thermal conductivity-and thus heat dissipation within the electrolyte-with respect to stoichiometric conditions. We find that the dependence of κ on temperature is also reduced, thus making it potentially easier to engineer devices that can safely and efficiently operate in a wide range of temperatures.

II. RESULTS
Before addressing thermal transport (Sec. II E), which represents the main goal of our work, we investigate the structural (Sec. II A), electronic (Sec. II B), vibrational (Sec. II C), and mechanical (Sec. II D) properties of Li 3 ClO, along with their temperature dependence, and we extensively compare our ab initio results with the existing literature. These results are not only preliminary to the calculation of transport coefficients, but they also make it possible to draw some general conclusions on the deployment of Li 3 ClO for mass production of SSEbased batteries: ductility, stiffness, the magnitude of the mismatch in the lattice constants or in the thermal expansion coefficients between the SSE and the electrodes govern the extent and amplitude of local stresses, particularly at the SSE-electrode junction, and affect the overall performance of Li 3 ClO in a real device by either hindering or inducing cracks, mechanical instabilities, or spurious electric fields.

A. Structural properties
The crystal structure of Li 3 ClO has the Pm3m perovskite space group. The dependence of the lattice parameter on temperature is obtained in the quasiharmonic approximation (QHA), 23 using the the vibrational frequencies computed as explained in Sec. II C and the Murnaghan equation of state. 24 The resulting zero-temperature value, explicitly accounting for zero-point vibrational effects, is 3.91Å, while the full temperature dependence is reported in Fig. 1  tice parameter and thermal expansion coefficient computed at room temperature (RT = 300 K) are 3.93Å and 2.77 · 10 −5 K −1 , respectively, in good agreement with the values previously obtained in other works: Zhang et al., 25 Emly et al., 26 and Wu et al. 27 reported an optimised lattice parameter without the zero-point of 3.85Å, 3.90Å, and 3.91Å, respectively; Braga et al. 1 provided both experimental and theoretical data in accordance with one another, namely 3.91Å. 1,28 For the linear thermal expansion coefficient, Zhang et al. 25 reported a value of 2.11 × 10 −5 K −1 using molecular dynamics simulations in the NPT ensemble, while Wu et al. 27 reported α = 3.12 × 10 −5 K −1 within the quasi harmonic approximation (QHA), in good agreement with our calculation. Braga et al. 1 found the larger value 4.65×10 −5 K −1 from the slope of the lattice parameter as a function of temperature.

B. Electronic band structure
The electronic band structure of Li 3 ClO (Fig. 2) exhibits a direct band gap of 6.46 eV at the M in the BZ at the HSE06 level of theory, in good agreement with the values obtained by Wu et al. 27 (6.46 eV) and by Emly et at. 26 , thus making Li 3 ClO a wide band-gap insulator and preluding a good electrochemical stability.

C. Vibrational properties
The phonon dispersion along high-symmetry lines in the BZ are plotted in Fig. 3 together with the corresponding vibrational density of states (VDOS). Notice the presence of a large longitudinal-transverse (LO-TO) splitting in the infrared-active mode at the centre of the BZ (Γ point), due to the non-analytic behaviour of the dynamical matrices induced by the long-range (Coulomb) tails of the inter-atomic interactions. It has been pointed out that a proper account of these tails is essential not only for a correct qualitative description of the vibrational spectrum in the optic region, but also for an accurate evaluation of the thermal conductivity, which is more sensitive to the low-frequency portion of the spectrum. 29 The effect of a proper account of the LO-TO splitting on the value of the lattice thermal conductivity of Li 3 ClO is examined in Appendix B 6. The stability at zero temperature of this material has been debated, with different authors claiming the system to be either unstable 30 or stable 27 : methods are employed, the authors find soft modes (imaginary frequencies) at the M and R points in the BZ using a 6 × 6 × 6 supercell. Since the instability is larger when a smaller (3 × 3 × 3) supercell is used, the occurrence of lattice instabilities may be an artefact due to lack of convergence even when a supercell as large as 6 × 6 × 6 unit cells is used in finite-difference calculation. Using perturbation theory at the zone-center of a supercell, 1 the authors of Ref. 27 do not find any dynamical instabilities. Our calculations, performed within density functional perturbation theory (DFPT), confirm the dynamical stability of this system.

D. Temperature dependence of the mechanical properties
Mechanical properties can have a deep influence on the fabrication of any device and, in particular, of alkaliion solid-state batteries. An important aspect is that there must be a good contact between electrolyte and electrodes during the activity of a device. A good SSE candidate must be able to sustain large strains 31 to prevent the interfaces with cathode and anode to deteriorate in response to the deformation thereof. In this light, excessive stiffness is a feature to be avoided. Another element to consider is the problem of Li deposition at the interface with the cathode, and the subsequent dendrite formation, 32 that affects especially liquids. A solid material may partially overcome this obstacle by fine-tuning elastic properties such as the shear modulus and Poisson's ratio, 33 even if it has been reported that dendrite growth can still occur for other reasons. 34 The three elastic moduli, namely the bulk modulus B, the Young modulus E, and the shear modulus G, measure how a material responds to volumetric, tensile, and shear stress, respectively; Poisson's ratio, ν, measures the strain response in a direction perpendicular to an applied strain; Pugh's ratio, B/G, is related to how materials are ductile as opposed to brittle. The Chung-Buessem 35 elastic anisotropy index, A G , and the universal elastic anisotropy, 36 A U , provide a measure to quantify the extent of anisotropy in the elastic response of a crystal, and are obtained from the Voigt and Reuss estimates of B and G. 36 The Grüneisen parameter γ is the variation of pressure with thermal energy density at constant volume. 37 All these quantities are computed as a function of temperature within the QHA. The isothermal and adiabatic elastic constants (C ij ) as a function of temperature are computed within the quasi-static approximation (QSA).
The temperature dependence of these quantities is shown in Fig. 4. Their values at RT and at T = 0 are reported in Table I, together with results already available in the literature. 25,31 While our results at T = 0 are comparable with the ones found in other works, at the working temperature the quantities have lower values, in accordance to the fact that the material becomes softer when temperature is increased. In particular, the temperature dependence of B and G entails that Pugh's ratio is also decreased, hinting at a greater brittleness 38 100 300 500 Temperature (K) of Li 3 ClO than previously predicted. 27,31 Furthermore, the computed value of A G is higher than what found in the literature, suggesting that Li 3 ClO is slightly more anisotropic than previously thought; this is confirmed by the value of A U , that allows a broader comparison with other materials. Thermal conductivity (κ) plays a key role in the quest for promising SSE for battery production. A high value of the thermal conductivity is a favourable quality for a SSE candidate to have, since efficient heat dissipation allows the manufacture of safer batteries that do not overheat.
A quite high value κ = 22.49 Wm −1 K −1 was previously reported. 27 This result was obtained from Slack's model (SM), 12 which is based on a rather crude approximation of the BTE 39-41 . We briefly review SM in Appendix A. At the same level of theory, we obtain 16.55 W m −1 K −1 that is 26% lower and it is likely to be closer to a realistic value. Yet, the coarseness of the approximation makes it impossible to draw conclusions: as already mentioned in the Introduction, SM is known, since its formulation dating back to 1960s, to be fairly accurate (to within 20%) only for simple systems such as rare-gas crystals, while it grossly fails for more complex materials. The main issue is to be addressed to the cubic dependence of κ on the Debye temperature Θ D of the crystal, which itself depends dramatically on the parameters used to estimate it. For instance, if the lattice parameter is reduced by less than 2%, i.e., from the value we compute (3.91Å) to a value found in the literature 25 (3.85Å), Θ D increases by 4% (from 630 K to 653 K) and κ by 12% (from 16.6 W m −1 K −1 to 18.5 W m −1 K −1 ). This and further approximations of the model-it is for instance insensible to the crystal structure of the material-makes it impossible to produce reliable estimates for κ within SM. For a realistic calculation the BTE must be treated in a detailed fashion relying on the explicit calculation (and inversion) of the phononic scattering-rate matrix to obtain the out-of-equilibrium occupation numbers of the phonons. A quick review of the BTE-based approach we follow is given in Appendix B.
The lattice thermal conductivity in the Single Mode Relaxation Time Approximation (SMRTA) is computed ab initio via genuine third-order perturbation theory as implemented in the D3Q code 42 distributed with Quantum ESPRESSO, and is displayed in Fig. 5 as a function of the temperature. For a comparison, we also report our calculation within the SM, which agrees qualitatively to the ab initio SMRTA result better than the existing literature. 27 We test the reliability of the classical FF in probing thermal properties of Li 3 ClO by employing it to compute the lattice thermal conductivity in the SMRTA. In general, a full solution of the linearised BTE may change significantly the computationally cheaper SMRTA result. In practice, we explicitly verified that the full solution reduces the SMRTA value of up to 0.5%, so we can safely keep all the calculations at the SMRTA level. Details can be found in Appendix B 4. Figure 5 shows that classical FF and AI results are in good agreement, suggesting that the chosen classical FF is suited for investigating the thermal transport properties of Li 3 ClO. As recently pointed out, 43,44 higher-order scattering events can drastically reduce the value of the lattice thermal conductivity; therefore, we test the effect of the inclusion of four-phonon scattering into the computation of κ. As shown in Fig. 6, four-phonon scattering is found to have a major role in determining the value of the lattice thermal conductivity, being able to reduce κ of at least 15% at RT, the reduction being larger for larger temperatures. Details on the calculation can be found in Appendix B 5. This facts calls for a method able to include higher-order scattering: molecular dynamics (MD) simulations together with GK linear response theory (see Sec. IV D), which automatically include all the orders of interactions among phonons, are a fitting candidate for this role. 2 The choice of MD comes with an additional benefit: 2 Notice that, in BTE, the occupation numbers follow the Bose- since fast ion conduction in SSE materials such as Li 3 ClO is due to diffusing defects, lattice dynamical methods are inadequate, as they require the atoms to have fixed equilibrium positions. MD simulations are not subject to this prerequisite and, therefore, they allow us to access thermal transport in the diffusive regime where fast ion diffusion is mediated by vacancy hopping.
According to the comprehensive studies of Mouta et al 46 and Lu et al, 47 LiCl Schottky pairs-divacancies generated by the removal of neutral groups of atoms that are deposited at the surface of the material-are more likely to appear and give rise to a high Li-ion mobility than both other Schottky (Li 2 O or Li 3 ClO vacancies) and Frenkel (Li vacancies and interstitials) defects. Thus, the nonstoichiometric systems we study are of the form n = k B T / ω, according to classical equipartition. This always leads to an underestimation of thermal conductivity in MD with respect to BTE, at the same level of the theory. An account on the relationship between GK and BTE approaches to thermal transport can be found in Ref. 45. Further details can be found in Appendix B 7.
varies between 0 (perfect crystal) and 0.1. The classical FF described in Sec. IV B is used to carry out the MD simulations and sample the heat flux according to Eq. (6), which in turn is employed to compute κ as in Eq. (7). To validate this, we compare classical FF calculations to estimates of κ extracted from i) a Car-Parrinello MD simulation, on which a DFT energy flux is computed according to Ref. 17 as implemented in Ref. 48; ii) a model obtained via machine learning techniques.
The ab initio GK thermal conductivity requires a computationally demanding DFT energy flux (it requires roughly twice the computational time of the AIMD simulation the trajectory is sampled from): therefore, a single calculation at RT is carried out for the sake of checking the accuracy of the FF. We obtain κ DFT = 6.5 ± 1.0 W m −1 K −1 at T = 306 K on a 4 × 4 × 4 supercell of crystalline Li 3 ClO, in close agreement with results from classical FF (see below).
As for the machine learning model, a DeepPot-SE NN 49,50 is trained on a 3 × 3 × 3 supercell of Li 3 ClO with a LiCl pair removed. Details on the model and its validation can be found in Appendix D. NN based MD simulations are carried out for x = 0 and x = 0.1 across the whole temperature range of interest. Thermal conductivity for the NN model is computed using the methodology developed in Ref. 51. In Fig. 7 we show a comparison between the GK thermal conductivity obtained with classical-FF-and the NN-potentials for these two systems. Results from NN simulations are in close agreement with the AI ones available at room temperature and, remarkably, with those obtained from classical FFs, over a broad temperature range, thus further validating the accuracy of the latter for the purposes of the present work.
Having thoroughly verified that the classical FF closely mimics AI-quality results, we use it to perform a systematic analysis of the thermal transport properties of Li 3 ClO in a broad range of temperatures and vacancy concentrations. MD simulations are carried out on 10 × 10 × 10 supercells. The desired vacancy concentration is obtained by randomly removing the corresponding number of LiCl pairs from the supercell by means of the Atomsk code. 52 For each temperature, we employed the temperature-dependent ab initio lattice parameter computed in Sec. II A in the QHA (independently of vacancy concentration). In the equilibration phase, the canonical ensemble is sampled via the CSVR thermostat 53 for 200 ps. At this point, the thermostat is removed and a microcanonical (N V E) production run of 5 ns is carried out to collect the desired data.
The dependence of the thermal conductivity on T and x is shown in Fig. 8. The numerical data measured in the MD simulations are then fitted to a simple function, Eq. (3), as described below. The functional form of such dependence should account for the Peierls-Boltzmann asymptotic 1/T behaviour (Eucken's law), 54 holding for crystals (x = 0) at large temperatures, but should also be able to take into account the breakdown of crystalline order when vacancies are present (x = 0). The presence of vacancies, while providing access to Li-ion conduction channels, on shorter time-scales establishes an effective local disorder, that would result in a contribution to the thermal conductivity describable at different levels of accuracy 18,55-57 . The simplest significant approach is to consider the Allen-Feldman model 56 of thermal conduction in harmonic glasses, where the thermal conductivity, κ AF , takes the form Here, λ is a label for the modes, n λ (T ) is the Bose-Einstein occupation number at temperature T , and D λ is the (temperature-independent) mode diffusivity; asymptotic expansion of (1) in the high temperature limit yields i.e., the leading order is constant in temperature. Notice that MD with classical nuclei sample, by definition, classical distributions: the occupation of a mode of energy ω is thus the first order expansion of Bose-Einstein distribution, n λ (T ) ≈ k B T / ω λ , and the heat capacity per mode reduces to the Dulong-Petit result, C λ ≈ k B . All in all, we take as fitting function for κ(T, x) where C Euck (x) and C AF (x) are vacancy-dependent fitting parameters, incorporating both the Eucken and the Allen-Feldman asymptotics. Their values are reported in The fitting parameters are also reported in Table II. Error bars represent standard deviations.
Table II and shown in the inset of Fig. 8. As expected, C AF is compatible with zero in the case x = 0, while sensibly different in the other cases. It is also comforting to notice that C Euck (x) vanishes as x increases, while C AF (x) saturates. To have a clear picture of the relative importance of these contributions, we explicitly report, in Fig. 9, the decomposition of κ vs T into the Eucken and AF terms at the different vacancy concentrations inspected in MD simulations. This picture suggests that, at least in this class of SSE, ionic diffusion may contribute very little to κ, which is instead dominated by the lattice component at low temperature, and by disorder at high temperature. This was recently investigated experimentally in Ref. 58 for Ag + fast-ion conductors.
Our GK calculations show that the thermal conductivity strongly depends on the presence of defects: in the considered temperature range, even a few-% concentration of defects is able to almost halve the κ with respect to its value for the perfect crystal. This behaviour is particularly evident at low temperature, due to the suppression of the Eucken law when x increases. Nonetheless, when compared to other candidates for battery-oriented SSE, Li 3 ClO is characterised by a relatively high thermal conductivity, and meets the requirements for safe heat dissipation and management.

III. CONCLUSIONS
We conclude with a summary of our results and some perspectives on the possible use of Li 3 ClO in realistic devices. In the first part of our work we have reported stateof-the art calculations on the structural and mechanical properties of Li 3 ClO, including the zero-point lattice contribution and non-analytical corrections. According to our calculations, Li 3 ClO is mechanically stable, and characterised by a larger brittleness and anisotropy than previously thought. From the thermal transport standpoint, we extensively showed, by comparing different levels of theory, that current estimates of the thermal conductivity are fully unreliable, overestimating the value obtained via more accurate theory by several times in the temperature window at which SSE are supposed to operate. We explicitly demonstrated that, due to strong anharmonicity, not even the full solution of the Peierls-Boltzmann transport equation with three-phonon scattering is sufficient and higher order contributions must be considered, which tend to further decrease the thermal conductivity of Li 3 ClO. We performed molecular dynamics simulations and employed the GK theory of linear response and newly developed data-analysis tools to automatically incorporate all scattering orders and deal with vacancies and charge carrier diffusion as well. We observed that increasing the number of vacancies induces local disorder, which leads to a glass-like behaviour where the thermal conductivity saturates at large T rather than vanishing like 1/T according to Eucken's law. Ionic diffusion per se seems instead to affect thermal transport only marginally. Our calculations indicate that, even if the capability of Li 3 ClO to dissipate heat in a realistic device must be scaled back, Li 3 ClO should still have a good thermal conductivity in a quite large temperature range and may be safely employed as a battery-oriented SSE. We hope that our work will encourage an experimental assessment of thermal transport properties of Li 3 ClO, and will also serve as a valid methodological reference, aimed at numerical-simulation practitioners, on how to investigate the thermal transport properties of generic SSE ionic conductors through state-of-the-art theoretical and data-analysis tools.

IV. METHODS
In this work, we follow a multi-method approach, hinging on a combination of classical simulations based on force fields (FF), ab initio (AI) computations based on density-functional (perturbartion) theory, and neuralnetwork (NN) modelling of ab initio data. In all cases, the accuracy of empirical FFs and NN potentials is thoroughly benchmarked against AI data.
Electronic, structural, mechanical, and vibrational properties are studied using DFT (Sec. IV A). The thermal conductivity is evaluated both through a Boltzmann-Peierls kinetic approach (Sec. IV C) and via the GK theory of linear response (Sec. IV D) and classical MD. BTE calculations are performed both ab initio, and using classical FFs (Sec. IV B). AI data are used as a reference, while FF simulations allow us to estimate the magnitude of finite-size and anharmonic effects, whenever AI calculations would be unfeasible, as well as the role of vacancies in thermal transport in non-stoichiometric conditions.

A. Electronic and atomic structure
All the electronic-structure calculations are performed within DFT and the plane-wave pseudopotential method using Quantum ESPRESSO 59-61 and SG15 Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials. 62 Structural optimization and lattice-dynamical calculations were performed using the generalised gradient approximation (GGA) to the exchange-correlation functional in the Perdew-Burke-Ernzerhof (PBE) flavour, 63 while the electronic band structure is computed using the Heyd-Scuseria-Ernzerhof (HSE06) 64 hybrid functional. The electronic structure is obtained starting from a self-consistent calculation on a coarse kpoint grid in the Brillouin zone (BZ) leveraging a fast implementation of the Fock-exchange operator, 65 followed by a transformation to the Wannier basis, performed with Wannier90. 66 and a subsequent non-self-consistent calculation to interpolate the energy bands. All calculations are performed using a plane-wave kinetic energy cutoff of 200 Ry (such an unusually high cutoff was needed for an optimal convergence of equilibrium thermal properties, vide infra) and the Brillouin Zone (BZ) sampling is done via a Monkhorst-Pack grid 67 of 12 × 12 × 12 points displaced by half a grid step along the (1, 1, 1) direction. Vibrational properties are computed via Density Functional Perturbation Theory (DFPT). 68 dynamical matrices have been calculated on a 7 × 7 × 7 q-point grid, and Fourier-interpolated on a 192 × 192 × 192 mesh to evaluate the temperature dependence of the free energy in the quasi-harmonic approximation, 23 with the aid of the thermo pw Quantum ESPRESSO driver. 69

B. Classical molecular dynamics
Molecular Dynamics (MD) simulations have been carried out with lammps, 70,71 using different kinds of force fields. Specifically, we adopt both an empirical FF and a NN interatomic potential trained and validated on ab initio data. We anticipate that the overall very good agreement between FF and NN simulations safely justifies the extended use of the empirical FF to study ther-mal transport across a wide range of temperatures and vacancy concentrations in a computationally affordable way.
The empirical force-field, which was fitted to structural and dielectric properties from experiments and ab initio calculations, is taken from Ref. 46. We adopt the simplified version based on a rigid-ion model, as opposed to core-shell models that are necessary when addressing specifically the electronic polarisability. The rigid-ion approximation has been used by other authors to investigate charge transport. 47,72 The short-range interatomic interaction is modelled by a Buckingham potential of the form where i and j label the two atoms separated by a distance r ij , while A ij , ρ ij and C ij are fitting parameters. The long-range (Coulomb) interaction among point charges is computed via the Ewald summation technique. 73 The time-step for the dynamics is of 1 fs. A Deep Potential-Smooth Edition (DeepPot-SE) NN model 49,50,74 is employed to accurately mimic AI atomic forces in a computationally affordable way. The model is trained on a set of configurations sampled from Car-Parrinello simulations 75 enriched with additional configurations chosen by the dpgen software. 76 Further details are provided in Appendix D.

C. Boltzmann transport equation
The BTE is solved by real-space techniques with the ShengBTE code 77 and its successor, almaBTE, 78 using harmonic and higher-order interatomic force constants (IFC) from both classical FFs, so as to make a comparison with GK-MD results handy, and first-principles methods based on DFT. In the first case, IFCs are obtained in real space using phonopy. 79 Details on size effects and the choice of the parameters can be found in Sections B 1, B 2 and B 3 of Appendix B. Quadratic and cubic force constants AI IFCs are obtained from second-and third-order DFPT, using the ph.x and D3Q 42,80 components of Quantum ESPRESSO, respectively.

D. Green-Kubo linear-response theory
In the linear-response regime, the thermal conductivity, κ, is defined as the ratio between the energy flux and the negative of the temperature gradient in the absence of any convection. For solids and one-component fluids this last prescription is trivially satisfied in practical MD simulations performed in the barycentric reference frame. Within the GK theory of linear response, the thermal conductivity is given by the celebrated GK formula [13][14][15][16]19 κ where the energy flux, J E (t), reads Here Ω is the system's volume, i is the energy assigned to the i-th atom, V i and R i its velocity and position, respectively, and F ji = − ∂ i ∂Rj . For a multi-component system such as the superionic phase of Li 3 ClO, Eq. (5) cannot be applied as it is, since the prescription of vanishing convection is not automatically satisfied in a standard MD simulation: in fact, keeping the barycentre of the whole system fixed does not imply that the barycentres of each atomic species S (with S = Li, Cl, O) stay fixed. Therefore, the mass fluxes of two of the three species, say J Li (t) = MLi Ω i∈Li V i (t) and J Cl (t) = M Cl Ω i∈Cl V i (t), where M S is the atomic mass of species S, i.e. the total momenta of each atomic species, are independent and in general non vanishing. Notice that, from a more fundamental point of view, even the concept of atomic energy i appearing in Eq. (6) is intrinsically flawed, in that i) there is no a priori "correct" decomposition of the total energy of a system of interacting atoms among its constituents, and ii) the zero of the atomic energies (i.e. the energy of isolated atoms) is arbitrary. A theoretically sound solution to both these problems has been provided in Refs. 17, 19, and 81, where it is rigorously proved that these apparent inconsistencies must disappear when a measurable quantity, such as κ, is correctly calculated, leading to the formulation of so-called invariance principles of transport coefficients. 20,82 Furthermore, a multivariate technique 83,84 for the analysis of the energy flux time-series obtained from MD simulations has been developed, which allows one to compute κ for multicomponent systems, like the superionic Li 3 ClO, in an efficient and rigorous way. We redirect the reader to Refs. 83 and 84 for a thorough description of the method. Suffice it here to say that the thermal conductivity for Li 3 ClO is estimated in terms the ω → 0 limit of the power spectrum of the energy flux, after removing its coupling to the mass fluxes. This is achieved by computing the Schur complement of the matrix whose entries are the power spectra between fluxes J A , J B , i.e. S AB (ω) = ∞ −∞ e iωt J A (t) · J B (0) dt, where A, B ∈ {E, Li, Cl} refer to energy flux, Li mass flux, and Cl mass flux, respectively: Since Li 3 ClO is a three-components material, two of the three mass fluxes are independent (the third being fixed by conservation of total momentum) and as such must, in general, be decoupled from the energy flux. In the case where some mass component were non-diffusive, they would contribute nothing to the value of the decoupled spectrum at zero frequency. Nonetheless, it was shown numerically that the decoupling at finite frequency allows a better estimate of the value at ω = 0. 84 . The numerical tools needed to evaluate the power spectra of time series of J E , and J Li and J Cl have been implemented in the open-source SporTran code, 85 which we extensively adopted in the present work.

DATA AVAILABILITY
Numerical data supporting the plots and relevant results within this paper, the neural network potential and the data set used to generate it are available on the Materials Cloud Platform 86 .
Appendix A: Slack model For a crystal with n atoms per unit cell, whose average atomic mass is M , thermal conductivity in Wm −1 K −1 is estimated as where Θ D is the Debye temperature in K, γ is the Grüneisen parameter, δ is the cubic root of the average atomic volume inÅ, and M is expressed in atomic mass units. The Debye temperature is computed directly from the Debye model, by evaluating the average sound velocity obtained from the angular average of the sound velocities, which are in turn calculated solving the wave equation for each propagation direction; the other quantities are described in the main text.
with ω qs the angular frequency of vibration of the normal mode labeled by (q, s). Under the assumption that the phonon population depends on position only through the temperature, in linear response to a small temperature gradient ∇ α T along direction α = x, y, z, the occupation number (B1) acquires a contribution proportional to ∇ α T : n qsα = n qs + λ qs ∇ α n qs = n qs + λ qs ∂n qs ∂T The quantity λ qs is called phonon mean free path. From the q-gradient of the phonon frequencies, one obtains the phonon group velocities v qsα = ∂ωqs ∂qα ; these, together with the phonon energies ω qs and the out-of-equilibrium occupation numbersñ qsα , are used to define the heat current per normal mode as An expression for the thermal conductivity is retrieved after having introduced the Fourier equation of diffusive conduction, J α = − α κ αα ∇ α T , that relates the macroscopic heat flux J α = q,s j qsα /(N Ω) to the temperature gradient. Thermal conductivity reads Thus, the computation of κ requires the knowledge of the out-of-equilibrium occupation number. The statistical mechanics of an out-of-equilibrium system can be described by BTE, whose expression is v qs · ∇T ∂n qs ∂T = ∂n qs ∂t scatt. .
The right-hand side is the scattering term, whose linearized form contains all the scattering rates relating three-phonons scattering events. The inverse of the matrix of scattering rates is in turn used to compute the phonon mean free paths, which are what is needed to obtain κ. The problem of computing κ is reduced to the calculation and inversion of the scattering rates matrix. This can be demanding, and often requires additional approximations, the simplest of which being the relaxation time approximation (RTA). In RTA, only the diagonal part of the scattering rates matrix is retained: the population of each out-of-equilibrium mode is taken to be interacting with a bath of modes at thermal equilibrium. A way to go beyond this picture is to recast the inversion problem to an iterative equation to be solved self-consistently. 88

Convergence of the IFC3 with the supercell
The Interatomic Force Constant (IFC) matrices are computed via finite differences on a supercell. The second order IFC (IFC2) is calculated inexpensively on a 6 × 6 × 6 supercell. The choice of the supercell for the third order IFC (IFC3) changes substantially the computational cost of the calculations, and requires more care. The IFC3 is computed as a function of the supercell size. This relationship is studied via the so-called "triangle plot". The triangle plot for the IFC3 computed with 2, 4, 6 and 8 unit cells in each dimension is shown in Fig. 10. In the horizontal axis there is the perimeter of the triangle whose vertices are positions of the three atoms involved in each block in the IFC3 matrix, a block being the 3 × 3 × 3 tensor ∂ 3 U ∂ri∂ry∂r k . In the vertical axis there are the matrix norms of each block. A supercell of size 6 × 6 × 6 is found to be sufficient to have a IFC3 converged up to ≈ 1%.

Convergence of the IFC with the number of nearest neighbour shells
The dependence on the number of included nearest neighbour (n.n.) shells of FC3 is studied both via the triangle plot and by directly computing the thermal conductivity tensor, κ, on a coarse 12 × 12 × 12 q-points mesh. The triangle plot for the FC3 computed with 6 × 6 × 6 supercell is shown in the left panel of Fig. 11, while on the right panel there is the relative change in κ for a given number of n.n. shells, denoted by κ nns , with respect to κ 0 , i.e. the thermal conductivity tensor obtained with the highest number of n.n. shells available. Frobenius matrix norms are employed, to take into account the off-diagonal elements as well. A number of n.n. shells equal to 9 is found to be sufficient to have a thermal conductivity tensor converged well below 1%.
The convergence of κ with the number of N q of qpoints is analysed by computing the thermal conductivity tensor for different reciprocal space meshes for different temperatures. The results are then analysed in terms of the relative change in κ for a given N q , denoted by κ Nq , with respect to κ 0 , i.e. the thermal conductivity tensor obtained with the highest N q available. Convergence is deemed as satisfactory when the change is consistently lower than 1%. From Fig. 12 one can see that a mesh of 20 × 20 × 20 q-points is sufficient to provide converged results.

Relaxation time approximation vs. full solution of the BTE
Calculations for selected temperatures are carried out to address the relative importance of the full solution of the linearised BTE with respect to the RTA result. As shown in Fig. 13, it appears that RTA deviates from the full solution of at most ≈ 0.6%. Thus, it is employed for the rest of the calculations.

Four-phonons contributions
Four-phonons calculations are carried out via the Fourphonon code 89 , an extension to ShengBTE. Since four-phonon calculations are quite expensive, scattering rates and thermal conductivity for three-and fourphonon processes are performed on a coarse mesh of 10 × 10 × 10 q-points with IFCs coming from a 4 × 4 × 4 supercell, and with a fourth-order interatomic force constant matrix (IFC4) that includes contributions from the first n.n. shell only. The scale parameter for Gaussian broadening is reduced to 0.1 from the default value of 1.0.

Effect of the inclusion of NAC on thermal conductivity
The effect on κ of the inclusion/exclusion of NAC to the dynamical matrix are addressed in Fig. 14. The exclusion of NAC is found to produce a variation of up to 5% to the value of κ in the temperature range of interest.

Equipartition vs. Bose-Einstein statistics
In a BTE approach one can choose the statitstics of the phonon occupation numbers. For quantum particles, it is the Bose-Einstein distibution MD with classical particles samples the first order expansion of (B6) at high temperature i.e. the equipartition law. Since the difference between Eqs. (B6) and (B7) is always negative at finite T , the thermal conductivity computed with equipartition will be lower than the one computed with Bose-Einstein distribution at the same level of the theory. We have checked how important this difference is for Li 3 ClO, by modifying ShengBTE to allow occupations according to (B7). The results are shown in Fig. 15. As expected being Eq. (B7) a high-temperature expansion of Eq. (B6), the difference is lower at higher temperatures. This contribution alone, that is around than 5% in most of the temperature range, is not sufficient to explain the difference between the BTE-based and MD-based values of κ. Appendix C: Data analysis

Estimation of thermal conductivity from MD simulations
Thermal conductivity is estimated from the power spectrum of the energy flux S EE , after removing its coupling to the mass fluxes, S coupl. , i.e.
Details on the method are found in Refs.83 and84. Being the spectrum an even function of ω, its low frequency part needs to be quadratic. Thus, the first portion (up to ≈ 0.1 THz) of the spectrum is fitted to the quadratic An example (T = 373 K and x = 0.005) of the spectrum of the energy flux decoupled from the mass fluxes is shown in Fig. 16.

Vacancy-dependent GK thermal conductivity
GK thermal conductivity from FF based MD simulations are fitted for each value of the vacancy concentration x. The fitting function is Eq. (3). The values of the fitting parameters C Euck and C AF , obtained via standard linear regression of log(κ) versus log(1/T ), are shown in Table II. The system used for the training of the NN is a 3 × 3 × 3 supercell of Li 3 ClO with a LiCl pair removed. The Root Mean Squared Errors (RMSEs) of the trained model are 0.11 ± 0.01 meV/atom and 15 ± 3 meV/Å for the energy and the forces, respectively. The uncertainty on the RM-SEs is the standard deviation of the RMSEs of four different models that differ by the choice of the initial random seed.
The locality test 90 is performed to assess the influence of long-range interactions on the NN model, that are not explicitly considered here. For a given atom, the atoms that are within a chosen radius (in PBCs) from it are kept fixed. A random perturbation is applied to the positions of all the other atoms placed outside the sphere. The force acting on the atom at the centre of the sphere are collected for a number of different random perturbations, computed both with the NN model and ab initio. The average deviation between the NN and ab initio forces, defined as σ f = |F AI − F N N | 2 , quantifies the dependence of the atom's properties on its neighbours; hence the name of Locality Test. This procedure is carried out for different values of the cutoff and for each atomic species in the system. The results obtained with a 4 × 4 × 4 supercell are shown in Fig. 17. Force deviations are at most compatible with the errors on the training, and exhibit a mildly decreasing dependence on the cutoff of the sphere. The error is computed as the standard error on the average of 10 equivalent configurations. This is consistent with the assumption that the local environment is dominant in determining the properties of an atom, and validates the model we employ. Phonon density of states is computed and compared to the ab initio results for different supercells. The results are satisfactory and are shown in Fig. 18.