Investigating finite-size effects in molecular dynamics simulations of ion diffusion, heat transport, and thermal motion in superionic materials Open Access

Superionic (SI) materials are characterized by a matrix of atoms arranged in a (crystalline or amorphous) solid, and of one (or more) mobile species, which diffuses above some critical temperature. The interest for SI materials has largely increased during the past decades, along with the quest for good candidates in the realization of solid-state batteries, where charge carriers like lithium ions move through a solid-state electrolytic matrix.^1–4 Moreover, the superionic phases of water and ammonia⁵ have been predicted to compose a large fraction of the outer cores of ice giant planets, such as Uranus and Neptune,^6,7 and many recent theoretical^8,9 and experimental^10,11 studies have focused on transport properties of materials becoming SI at planetary conditions, to study the evolution of these celestial bodies.^12,13

The complexity and the variety of new SI materials naturally imply that, from the computational material-science standpoint, much effort is devoted to the prediction of static and dynamical properties and their microscopic description by means of atomistic simulations. In particular, molecular dynamics (MD) simulations are needed whenever dynamical properties (such as transport coefficients and correlation functions) are investigated.¹⁴ 

Due to the chemical complexity of many of these materials, ab initio MD simulations are often performed, whose computational cost currently limits the simulation box to a few hundred atoms at most. This limitation poses serious questions on the role of finite-size effects (FSE) in the characterization of the physical properties of a SI material. For instance, recent tests were run on Li₁₀GeP₂S₁₂-type SI conductors,¹⁵ where machine-learning interatomic potentials trained on ab initio calculations allowed for simulations that are inaccessible to ab initio MD.¹⁶ These calculations showed that simulation boxes containing even some hundred atoms overestimate the Li-ion diffusivity by one order of magnitude with respect to the largest size considered (1600 atoms). This may have dramatic consequences in calculations aiming to find the best SI conductors for realistic devices.^17–19 

In liquids, FSE affecting particle diffusion have been extensively investigated (see, e.g., the recent review of Ref. 20): The hydrodynamics arguments by Yeh and Hummer^21,22 suggest that, for a cubic simulation box under periodic boundary conditions (PBC) and for a given particle density, the diffusivity of the liquid can be corrected by a factor proportional to the inverse of the box side, and that the proportionality coefficient only depends on geometric factors, the temperature, and the viscosity of the liquid, which is usually largely independent of the size.^22,23 Recent works evidenced that the application of the Yeh and Hummer correction is justified also for multicomponent liquids and ionic melts.^24–26 Nevertheless, the hydrodynamics equations of an SI material are different from (and more complicated than) those of simple liquids.^27,28 For instance, transverse modes of the lattice survive, the static shear modulus is non-vanishing, and the atoms of the mobile species diffuse via hopping mechanisms that are qualitatively different from the motion of particles in a simple fluid:²⁹ the Yeh–Hummer arguments are, therefore, inappropriate to account for FSE in the diffusion of charge carriers in SI materials. Furthermore, while general trends for the FSE on heat transport in solids and liquids have been reported in the literature,^31–32 such an analysis is currently missing for equilibrium MD simulations of thermal conduction in SI materials.

This article aims at investigating of the FSE in the calculation of relevant static and dynamical properties of SI materials via equilibrium MD simulations. I focus on simple yet paradigmatic examples of type-I and type-II SI conductors, which can be effectively described in terms of largely validated empirical potentials. I restrict my analysis to systems with perfect stoichiometry. Furthermore, I consider cubic simulation boxes where the unit cell, of lattice parameter a, is equally replicated ℓ times in all the three spatial directions, to avoid additional and non-trivial effects that arise, even for simple liquids, in the case of anisotropic replications.³³ 

Section II discusses the SI materials I selected to investigate. Section III gives methodological details on the equilibrium MD simulations that I performed. Section IV provides the main results of the calculations, by analyzing the size dependence of the specific heat capacity, the mobile-ion diffusivity, the thermal conductivity, and the Debye–Waller B-factors, for each of the selected SI materials. Finally, I draw general conclusions in Sec. V.

I choose the fluorite-structure materials PbF₂, CaF₂, and UO₂, as simple, yet prototypical examples of SI materials where FSE should be particularly relevant: all the energetically equivalent regular sites of the mobile species—the anions—are occupied, and the hopping of one diffusing anion can eventually occur only with a net hopping of other anions, since anion diffusion “occurs by discrete hops between regular sites,” and anions “do not reside in a well-defined manner on the cube-centre sites” (verbatim from Ref. 29. See also Ref. 34 for an insightful analysis of anion distribution in fluorites, and Ref. 35 for a recent, comprehensive study of cooperative F dynamics in PbF₂). Such a concerted hopping mechanism may easily extend to more than one lattice constant, leading to a size dependence. Furthermore, these materials (see Fig. 1, left panels) are characterized by a continuous order–disorder transition to the SI phase with no structural change in the crystalline structure of the non-diffusive species (type II superionic materials). The diffusion mechanism depends on the specific temperature regime, and, in particular, whether the system is in the SI phase or not.³⁶ 

The finite size of the sample is known to cause a shift in the critical temperature of second-order phase transitions.³⁷ Therefore, for a given temperature and particle density, a small simulation box may be in a different thermodynamic phase with respect to a larger one, with different diffusive mechanisms and a dramatic effect on diffusion. This is likely to be the case in all those materials where the diffusion mechanisms are strongly dependent on T. As we shall see in Sec. IV, this tangling between the diffusion mechanisms (hydrodynamics) and the phase of the system (thermodynamics) is responsible for dramatic FSE on ionic transport in these materials.

I also investigate a different material, the cubic phase of silver iodide (α-AgI), as a typical example of a system where FSE on the diffusion coefficient should be less relevant (Fig. 1, right panels): in contrast with fluorite-structure materials, in α-AgI the large degeneracy of equivalent positions that one Ag ion—the diffusive species—may take inside a unit cell results in a large freedom in the choice of empty sites (empty red circles in Fig. 1) that a selected Ag ion can hop to: the temperature affects the probability that hopping occurs, but not the general mechanism of diffusion. Moreover, due to the large degeneracy of empty regular sites, there is no need for the hopping of one Ag cation to be accompanied by the hopping of other neighbor Ag cations. The α-AgI phase is superionic and is reached after a first-order phase transition (at $≈420$ K at ambient pressure³⁸) from the hexagonal, non-conducting β-phase. The sudden, discontinuous change in the ionic diffusion at the phase transition makes AgI a type-I SI conductor.

The choice of these materials is also motivated by the availability, in the scientific literature, of largely validated empirical potentials that proved to qualitatively describe the ion diffusion mechanisms as well as the static properties of these systems.³⁹ These potentials make it possible to run, at a feasible computational cost, reliably long simulations (∼ns) at different temperatures and sizes, to extract general behaviors.

For PbF₂ and CaF₂, I employ the following two-body potential, combination of Coulomb and Buckingham potentials:

with the parameters optimized in Ref. 40 (PbF₂) and in Ref. 29 (CaF₂), reported in Table I. The success of this potential in the microscopic study of ionic diffusion in these materials dates back to the 1980s. For UO₂, I employ a recently developed potential, described in Refs. 41 and 42, which combines Buckingham, Coulomb, and Morse potentials to treat two-body interactions, as well as the embedded atom method (EAM) to account for many-body interactions. I refer the interested reader to the original literature and to the Materials Cloud repository of this work for the explicit parametrization values (see “Data Availability”).

Finally, for AgI, I use the following combination of Coulomb and Morse potentials:

with the parameters of Ref. 43, reported in Table II, derived via the Chen–Möbius lattice inversion method from ab initio calculations of cohesive energies. This potential displays good agreement with experiments⁴⁴ concerning static properties of different phases, as well as Ag diffusivity in α-AgI, and it is also consistent, in a wide temperature range, with the widely used Parrinello–Vashishta–Rahman empirical potential.⁴⁵ 

All simulations are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software.⁴⁶ It has the great advantage, with respect to other MD codes, that force computation is not subject to minimum image conventions, and one can use cutoffs larger than half the simulation domain size, thanks to the inclusion of “ghost” atoms. This is particularly important for the purpose of this work, where small boxes are needed for the FSE analysis, yet the cutoff radius should be the same to avoid changes in the form of the potential. The long-range interactions are included by means of the Ewald-summation technique in MD simulations of PbF₂, CaF₂, and AgI, and with the particle-particle particle-mesh (PPPM) method for the MD simulations of UO₂.⁴⁷ The simulations of PbF₂, CaF₂, and α-AgI are run with a MD timestep of 4 fs. For UO₂, the MD time step is set to 2 fs. The trajectories (⁠$≈800$ ps), from which the mean square displacements of the atoms are computed, are sampled each 10 MD time steps. For a given material, all constant-NVT (canonical) and constant-NVE (microcanonical) simulations (N is the number of particles, V is the cell volume, T is the temperature, and E is the total energy) are run at fixed lattice constant, a, irrespective of the temperature, i.e., no thermal expansion is considered for simplicity. I employ the following lattice constants: a_PbF2 = 6.056 Å (value at T = 792 K in Ref. 40), a_CaF2 = 5.712 Å,²⁹ a_UO2 = 5.65 Å (online material of Ref. 42), and a_αAgI = 5.37Å.⁴³ The NVT simulations are run with a Bussi–Donadio–Parrinello stochastic-velocity-rescaling (SVR) thermostat,⁴⁸ as implemented in LAMMPS. The temperature damping parameter of the SVR thermostat is set to 100 MD timesteps. Further details on the equilibration procedures and on the cutoffs employed are reported, for the sake of reproducibility, in the input scripts of the simulations, which are all available in the Materials Cloud repository of this work (see “Data Availability”).

I proceed by investigating FSE for physical quantities important for superionics, namely, the specific isochoric heat capacity (Sec. IV A), the diffusivity of the mobile species (Sec. IV B), the thermal conductivity (Sec. IV C), and the Debye–Waller B-factor of the non-diffusive ions of the solid matrix (Sec. IV D).

The isochoric molar specific heat capacity, c_V, is obtained from the finite-difference derivative of the average energy with respect to the temperature and displayed vs T in Fig. 2 for the four materials considered. Let us first focus on the fluorite-structure materials (first three panels of Fig. 2). For sufficiently large cells (ℓ ≥ 2, i.e., N ≥ 96), the heat capacity clearly displays a peak. Experimentally, a peak in the heat capacity—at some high critical temperature, yet below the melting point—has been observed via heat-content measurements of fluorite-structure materials, and associated with a transition that is not of the first order; such anomaly accompanies a sensible onset of electrical conduction, indicating a transition to the superionic phase.^49,50 At a large enough size, the calculated c_V is in fairly good agreement with experiments. As mentioned in Sec. II, it is known since the late 1960s that the effect of a finite size is to broaden a second-order transition and to shift the (pseudo)critical temperature T_c(ℓ) with respect to its thermodynamic-limit value T_c(∞).⁵¹ Whether the shift is positive or negative depends, among other factors, on the boundary conditions: usually, in PBC, T_c(ℓ) > T_c(∞), as a result of extra “communication” via paths that encircle the torus (verbatim from Ref. 51). This is, in fact, the behavior observed in Fig. 2, where the peak shifts toward lower temperatures as the size of the simulation is increased, in agreement with the existing literature.⁵² In the smallest cell, ℓ = 1 and N = 12, no peak is observed and, in the temperature range that I consider, the heat capacity is always sensibly lower than the one obtained with larger simulation boxes. I remark that, in agreement with Ref. 42, the melting of UO₂, predicted for this potential at $≈3100$ K⁴¹ via a moving interface method, is not observed in these simulations, where the lattice parameter is kept fixed to its value at $≈2600$ K⁵³ (see Sec. III B).

In striking contrast with fluorite structure materials, for α-AgI, no significant size effect is observed, in line with its intrinsic superionic structure. The link between FSE and the onset of a superionic phase transition is clearly highlighted from the analysis of the diffusivity of the mobile species, as discussed below.

The diffusion coefficient of the mobile species is the most characterizing quantity of superionic conductors: its behavior at the superionic transition dictates the classification of SI materials (see Fig. 1). Ion diffusion is usually described as an Arrhenius-like process,

where k is Boltzmann’s constant, A is a prefactor with the dimensions of a diffusivity, and E_a is the activation energy of the hopping mechanism leading to particle diffusion. In SI materials, in general, for a specific regime of diffusion, A and E_a weakly depend on the temperature: it makes sense, therefore, to plot the logarithm of the diffusion coefficient against the inverse temperature (in the so-called Arrhenius plot) to highlight significant changes in the diffusion mechanism. E_a represents the slope of the Arrhenius plot, and A is its intercept.

Figure 3 displays the Arrhenius plots of the diffusivity of the mobile species for the materials considered in this work, at different system sizes. The reported values for D (markers) are obtained from the slope of the mean square displacement (MSD) of the mobile species at large enough time,

where N_d is the number of atoms of the diffusive species and r_i(t) is the position of the ith atom of the mobile species at time t. A linear fit on the MSD at large t was used to obtain D. A block analysis, with four blocks of $≈170$ ps each, is performed to extract the uncertainty on the MSD. This uncertainty propagates to the asymptotic-time slope of the MSD and eventually to the diffusion coefficient. The uncertainty on D is reported as the shaded area in the Arrhenius plots. The reported D is computed in the “laboratory” reference frame, where the center of mass of the non-mobile species is fixed, although the MD simulations are performed in the barycentric reference frame, where the total momentum vanishes. The relations between the different, reference-frame dependent definitions of the diffusion coefficients are described in the supplementary material, Sec. S2 A. To check whether the SVR thermostat affects the estimate of D, I also run NVE simulations, previously equilibrated at the target temperature. Diffusion coefficients extracted from NVE and NVT simulations display fully compatible Arrhenius plots for each simulation cell size, as reported in the supplementary material, Fig. S9.

In PbF₂, CaF₂, and UO₂, I find that D(T) strongly depends on the system size. For instance, the diffusivity for N = 96 is greater (smaller) than the large-N values at low (high) T, in agreement with previous observations.²⁹ Furthermore, the kink in the Arrhenius plot at T_c, characterizing type II materials, is observable only for N ≥ 324 (i.e., ℓ ≥ 3). This is quantified by the values of E_a(ℓ) (the slope of the Arrhenius plot) extracted from the fit of D to Eq. (3) below and above the size-dependent critical temperature, T_c(ℓ), and reported in Table III: only for N ≥ 324, a sensible difference between E_a(T < T_c) and E_a(T > T_c) is observed. Notice, as well, that the kink in the Arrhenius plot, associated with the SI transition, moves toward lower temperatures as the size is increased, in accordance to the existing literature,⁵⁴ and with the shift in T_c extracted from the maximum of the heat capacity (Fig. 2).

As shown in the supplementary material, Figs. S4 and S5, all these results are in qualitative agreement also with frozen-matrix simulations,⁵⁵ where the non-diffusive ions of the solid matrix are kept fixed to their equilibrium position, and with simulations employing a short-range version of the Coulomb interaction.⁵⁶ This agreement confirms the proposed picture whereby the FSE on the diffusivity of mobile ions are mainly imputable to geometric factors (such as the degeneracy of empty sites and the extra interatomic communication paths occurring in PBC), rather than to the details of the potential or the vibrations of the solid matrix. This is also justified by a set of simulations run on non-stoichiometric lead fluoride, where a size-independent concentration of empty sites is generated by removing a set of randomly selected F⁻ ions accordingly: the presence of available empty sites favors hopping on a more local scale than in systems with perfect stoichiometry, where, below the SI transition, all the regular sites are occupied and the hopping of one F⁻ anion can only occur with a net hopping of other F⁻ anions, and results in a drastic reduction of FSE on fluorine-ion diffusion, as shown quantitatively in  Appendix B.

In striking contrast with fluorite-structure materials, in α-AgI, the values of the diffusivity of Ag ions at different cell sizes (even for very small cells) are all consistent with each other. The activation energies range from 93 ± 3 to 100 ± 7 meV (without any particular trend connected to the box size) and agree with the experimental value E_a = 94.97 meV of Ref. 44. As suggested above, I ascribe the size independence of the diffusivity to the large degeneracy of equivalent Ag ion sites within the conventional unit cell of α-AgI: in this superionic phase, Ag-ion hopping occurs on a smaller length scale than in fluorite-structure materials, where, instead, the hopping of one ion can occur only if accompanied by the concerted hopping of other ions, since—apart from the short transient of the jump—all the regular sites are occupied, as confirmed also in the literature (see, e.g., Refs. 40 and 57).

Figure 4 shows the temperature and size dependence of the thermal conductivity, κ, for the materials considered, extracted from NVT MD simulations according to the Green–Kubo (GK) theory of linear response for multicomponent, diffusive systems, as described in  Appendix A. Even if a mode-based analysis^58,59 of heat conduction would be required for a quantitative assessment of the role of disorder on phonon propagation, here, I employ the GK formula and MD simulations to include the role of diffusing ions, which, by construction, is not considered in mode-based calculations.

Let us first focus, again, on fluorite-structure materials. At large T, κ is almost independent of T due to the increasing disorder that accompanies the diffusion of mobile ions, suppressing phonon propagation. This is a manifestation of the Allen-Feldman regime,⁶⁰ and a typical feature also observed in other SI materials, such as solid-state electrolytes.⁶¹ The FSE of κ are system dependent, the values for PbF₂ being almost converged for N = 96, in contrast with CaF₂ and UO₂. Nonetheless, a general trend can be observed. The N = 12 (ℓ = 1) cell dramatically fails in describing heat transport and displays a non-physical, slight increase of κ with temperature. For ℓ > 1, the Allen-Feldman limit is achieved at lower temperatures for smaller simulation boxes, in agreement with the larger diffusivity at low T observed for small size, and the higher degree of disorder in smaller cells.⁶² Furthermore, the presence of a single defect is not supposed to strongly perturb a large system, while it would dramatically affect a few-atom cell, suppressing phonon propagation in favor of the Allen-Feldman limit. The activation of diffusion with the onset of significant disorder is much more relevant for heat transport than the actual SI phase transition: at T = T_c(ℓ) (vertical dotted lines), no particular feature of κ(T) is, in fact, observed.

Overall, the large-cell values are in good agreement with existing literature: κ_PbF2 = 1.4 Wm⁻¹ K⁻¹ at 300 K, Ref. 63 (experimental); κ_CaF2 = 1.46 ± 0.29 Wm⁻¹ K⁻¹ at T = 1694 K, Ref. 65 (numerical simulation, same density, and potential); κ_UO2 = 1.5 Wm⁻¹ K⁻¹ at T = 3000 K, Ref. 64 (numerical simulation, close density but different potential). I remark that, in UO₂, for T ≳ 2000 K, the growing contribution of electrons to heat conduction must be added for a comparison with experiments.

The thermal conductivity of α-AgI is almost constant in the considered temperature range. This indicates that phonon propagation is always suppressed by disorder in favor of the AF regime in this intrinsically superionic phase. All simulations with N ≥ 108 are fully compatible. For the small cell, N = 32, κ is lower, though only by 10%–15%, than the converged value. The results are in good agreement with the experimental value of 0.17 Wm⁻¹ K⁻¹ at T ≈ 500 K.⁶⁵ 

FSE also affect the dynamics and thermal vibrations of the solid matrix, which I investigate in this section in terms of the mean-square-displacement of the non-diffusive (n.d.) ions with respect to their equilibrium position $⟨u2⟩n.d.≡⟨[r(t)−req]2⟩n.d.=12limt→∞⟨r(t)−r(0)2⟩n.d.$⁠, which is a time-independent quantity. Furthermore, following a standard convention, I shall employ the so-called B-factor, $B≡8π23u2n.d.$⁠, entering the Debye–Waller factor that dictates the attenuation of x ray or neutron scattering in experiments.

For SI phases, MD simulations are needed to compute the B-factor, since normal-mode-based approaches⁶⁶ cannot be applied. Nonetheless, the values obtained from MD simulations are known to slowly converge with size. Early calculations for cubic hard-sphere solids show FSE corrections that follow a 1/ℓ law.⁶⁷ This law stems from the minimum frequency that can be sampled in a finite box size. An extensive derivation of this FSE is provided in the supplementary material, with an example on FCC solid argon. Interestingly, the same trend is observed in my simulations on SI materials. Figure 5 shows the B-factor of PbF₂, as a function of 1/ℓ, at different temperatures. The linear behavior is evident, and a linear fit of the data can be used to extrapolate the B-factor for infinite size, B(∞) (blue crosses). Note that the B-factor is here re-scaled by T and other geometric factors so that the slope represents the inverse of an effective elastic modulus of the material (see the supplementary material, Sec. V A). Only the smallest box considered (i.e., the conventional unit cell, ℓ = 1, with N = 12 atoms) deviates from the 1/ℓ behavior, not surprisingly, as it fails in describing most of the quantities analyzed so far, and where higher order FSE may enter.

The 1/ℓ behavior is observed in all the systems considered. An exception is the temperature range between 2600 and 2900 K for UO₂: in fact, the drastic change in the elastic properties of UO₂, and, therefore, of B, at the SI phase transition, occurs at different critical temperatures for different sizes (see Sec. IV A). Figure 6 displays the B-factor as a function of temperature for different systems and sizes. The blue, shaded area represents B(∞) with its uncertainty. Note that significant FSE are here observed not only in fluorite-structure materials but also in α-AgI. In fact, the very reason why the B-factor exhibits FSE—the minimum mode frequency that can be probed for a given cell—is very general and not system dependent. Once again, these calculations show that wrong results may be obtained even from an accurate description of interatomic interactions in a MD simulation (e.g., by computing forces ab initio), whenever FSE are not correctly accounted for.

Finite size effects (FSE) in superionic materials are, in general, not negligible and strongly depend on the specific system and property examined. In materials that possess intrinsically a high availability of degenerate hopping sites, such as α-AgI, the FSE on the diffusivity—the key quantity of superionic conductors—are weak. Things change, instead, whenever the mechanism of diffusion strongly depends on the temperature range, like in fluorites. These systems display an order–disorder transition to the superionic regime without a net change neither in the lattice structure nor in the available hopping sites, which coincide with those occupied by the mobile ions in the non-superionic phase. This work shows that, in these materials, the diffusivity is strongly affected by the system size, even qualitatively: the change of activation energy of the diffusion process at the SI critical temperature, resulting in the typical kink in the Arrhenius plot of the diffusivity, can be observed only above a certain threshold size. The order–disorder critical temperature is also largely affected by FSE, being larger for smaller sizes, as indicated by the shift of the maximum of the specific heat capacity. Therefore, in these materials, tangling between hydrodynamics (diffusion processes) and thermodynamics (the specific phase of a material) is responsible for changes in the diffusivity of even some orders of magnitude, at a given temperature, depending on the size. The thermal conductivity is also affected by FSE, mainly due to the role of disorder in hindering phonon propagation: this effect is larger at smaller simulation boxes, where the Allen-Feldman regime is reached at lower temperatures. The transition to the SI phase does not seem to strongly affect κ, instead. In general, an a priori determination on the minimum size that is sufficient to obtain satisfactory results is hardly feasible, due to strong dependence of FSE on the specific superionic material. My analysis suggests, a posteriori, that, for fluorite-structure materials, a minimum of N = 324 atoms is needed to correctly capture at least the main qualitative features of particle diffusion and heat transport, while for the α-AgI phase, even a relatively small cell of N = 108 atoms is sufficient to obtain a quantitative convergence of all the analyzed properties.

FSE also affect the Debye–Waller B-factor incorporating the thermal motion of the non-diffusive ions of the solid matrix. The 1/ℓ (or 1/N^1/3) behavior, predicted for simple solids, is observed also for the SI materials considered. This is motivated by the minimum frequency of oscillation, which can be captured by a simulation of a given size, irrespective of the specific system or phase considered.

Increasingly reliable interatomic potentials have provided, during the past few years, a systematically more accurate description of SI materials. However, this work clearly shows that accurate potentials alone are not sufficient to converge the calculation of many key properties, if not complemented with a full treatment of FSE. From the application point of view, particular attention must be paid in calculations aiming to compare different candidates for realistic devices, such as solid-state batteries. A given simulation size may be sufficient for some superionic materials but not for others.

See the supplementary material for more details about (i) the calculation of the heat capacity, (ii) the calculation of the diffusivity (role of the reference frame, of modified Coulomb interactions and of vibrations of the solid matrix), (iii) a comparison of the results obtained with NVT and NVE simulations, (iv) the calculation of the thermal conductivity, (v) the derivation of the 1/ℓ law for FSE of the B-factor and a simple application to solid argon, and (vi) the heat capacity of the defected structure of  Appendix B.

I thank Michele Ceriotti, Loris Ercole, Alfredo Fiorentino, Lorenzo Gigli, and Paolo Pegolo for insightful discussions and fruitful comments on the manuscript. I acknowledge funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Action IF-EF-ST (Grant Agreement No. 101018557; TRANQUIL).

The author has no conflicts to disclose.

The data that support the plots and relevant results within this paper are available on the Materials Cloud platform.⁶⁸ See DOI: https://doi.org/10.24435/materialscloud:jy-tw.

The thermal conductivity, κ, is the proportionality coefficient between the energy flux and the (negative of the) temperature gradient in the absence of any convection. For a two-component system, like the SI materials studied in this work, characterized by the presence of one diffusive species,^61,69 the Green–Kubo theory of linear response allows to extract κ from MD simulations as the zero frequency component of

where

is the power spectrum of the fluxes J_A(t) and J_B(t). The flux J_e(t) is the total energy flux, here computed via the compute heat/flux command of LAMMPS, while $Jd(t)=1V∑i=1Ndṙi$ is the convective flux of the diffusive species. Both J_e(t) and J_d(t) have been computed in the barycentric reference frame. A multivariate technique^72,73 for the cepstral analysis of time-series of the energy flux obtained from MD simulations has been employed in this work,^70,71 allowing one to compute κ, efficiently and rigorously, for multicomponent superionic materials. As for the diffusivity, I checked any possible influence of the thermostat on the value of κ, by repeating several simulations in the NVE ensemble, previously equilibrated at the desired temperature. The results from NVT and NVE are fully compatible (see the supplementary material, Fig. S10).

The second term between square brackets in Eq. (A1) represents the contribution to heat flux due to convection. It must be removed from the first term, S_ee, to correctly calculate κ, i.e., the coefficient of thermal conduction.^64,69,74,75. Its effects on the thermal conductivity of CaF₂ at different temperatures and system sizes are shown in Fig. 7: the full calculation, employing Eq. (A1) (solid lines and error bars), coincides with the single component calculation (shaded areas), that is, $κ=V6kT2See$⁠, only at low temperatures, where diffusion is negligible. Note that the departure from the multicomponent value occurs at lower temperature for small sizes, in agreement with overestimation of D in small boxes at low T.

I investigated the role of a higher number of accessible empty sites, to which a F⁻ ion can jump, by considering a strongly defected lead fluoride system where, for each simulation, ℓ³ fluorine ions are randomly selected and removed from the ℓ × ℓ × ℓ supercell with perfect stoichiometry, ℓ being, as usual, the number of replicas of the conventional cell. In this way, the concentration of F⁻ vacancies is set to 1/8, and it is independent of ℓ. Therefore, the simulation cells possess 4ℓ³ lead ions and 7ℓ³ fluorine ions. The selection of the F⁻ ions to remove, and the generation of resulting defected cell, is performed with the code ATOMSK.⁷⁶ A different seed for the random selection is used for each temperature and size. After a long initial equilibration (400 ps), a production run in the NVT ensemble is performed with LAMMPS. A uniform, neutralizing, background charge distribution is implicitly applied in the simulations employing Ewald’s method.

Although such a high concentration of F⁻ vacancies is hardly attainable in actual experimental samples, this model serves as a test bench to validate the general picture whereby the presence of empty sites tends to reduce FSE, as extensively discussed in the main manuscript. In fact, as shown in Fig. 8, in these simulations, convergence in the Arrhenius plots is reached already for the ℓ = 2 cell (N = 88 atoms), in striking contrast to cells with perfect stoichiometry. Another interesting difference with the perfect stoichiometry case is that the activation energy below the transition to the superionic phase is significantly lower than that above the SI transition (whose value is itself quite close to the one for perfect stoichiometry): the presence of empty sites facilitates F⁻ ion diffusion below the transition to the SI phase, while it becomes much less relevant above it, where the system is globally disordered.