Mechanism of Charge Transport in Lithium Thiophosphate

Lithium ortho-thiophosphate (Li3PS4) has emerged as a promising candidate for solid-state electrolyte batteries, thanks to its highly conductive phases, cheap components, and large electrochemical stability range. Nonetheless, the microscopic mechanisms of Li-ion transport in Li3PS4 are far from being fully understood, the role of PS4 dynamics in charge transport still being controversial. In this work, we build machine learning potentials targeting state-of-the-art DFT references (PBEsol, r2SCAN, and PBE0) to tackle this problem in all known phases of Li3PS4 (α, β, and γ), for large system sizes and time scales. We discuss the physical origin of the observed superionic behavior of Li3PS4: the activation of PS4 flipping drives a structural transition to a highly conductive phase, characterized by an increase in Li-site availability and by a drastic reduction in the activation energy of Li-ion diffusion. We also rule out any paddle-wheel effects of PS4 tetrahedra in the superionic phases—previously claimed to enhance Li-ion diffusion—due to the orders-of-magnitude difference between the rate of PS4 flips and Li-ion hops at all temperatures below melting. We finally elucidate the role of interionic dynamical correlations in charge transport, by highlighting the failure of the Nernst–Einstein approximation to estimate the electrical conductivity. Our results show a strong dependence on the target DFT reference, with PBE0 yielding the best quantitative agreement with experimental measurements not only for the electronic band gap but also for the electrical conductivity of β- and α-Li3PS4.


I. INTRODUCTION
The growing demand for portable electronic products and electric vehicles has stimulated the creation of energy storage systems that offers better safety and higher energy density than current Li-ion battery systems [1].While commercial Li-ion batteries use organic liquid electrolytes and additives to achieve a high working voltage [2,3], these materials pose safety concerns due to their flammability and susceptibility to thermal runaway [4,5].To address these issues, researchers are developing all-solid-state batteries (ASSBs) with inorganic solid electrolytes (SEs) to provide a sustainable solution for energy storage, exploiting their expected longer lifespan and improved energy efficiency [6,7].Many families of SEs have been considered and studied during these years [8][9][10].Sulfides are recognised as uniquely promising materials due to their remarkable mechanical stability and room-temperature ionic conductivity [11][12][13][14][15][16].In particular, the family of lithium thiophosphates (LPS), with its archetypal Li 3 PS 4 compound, is widely recognised as one of the most promising family of sulfide electrolytes and it has been the subject of many experimental and computational studies [1,[17][18][19][20][21][22][23][24][25][26][27].
Li 3 PS 4 has three main polymorphs: α-Li 3 PS 4 (with space group Cmcm [28]), β-Li 3 PS 4 (P nma) and γ-Li 3 PS 4 (P mn2 1 ).Whereas the γ polymorph is the most stable at room temperature, it also exhibits low roomtemperature ionic conductivity (≈ 3 × 10 −7 S cm −1 [20]).The system transforms into the metastable β- † These authors equally contributed to this work * michele.ceriotti@epfl.chpolymorph at 573 K and then into the α-polymorph at 746 K [20].Despite their great relevance, in the past years computational studies have been limited by the use of empirical potentials [23,24] and ab-initio molecular dynamics (AIMD) [22], based on density functional theory (DFT) with generalised gradient approximation (GGA) [29,30].The former can provide useful mechanistic insights, but fail to correctly predict the activation energies of the conductive phases [23] and are inherently limited in their accuracy and transferability.Quantum mechanical approaches, on the other hand, are more accurate but they are burdened by a higher computational cost, which hinders their applicability to realistic systems.For example, recent studies based on AIMD-PBE simulations attributed the superionic conductivity of glassy 75 Li 2 S-25 P 2 S 5 and that of bulk β-Li 3 PS 4 to the presence of fast cation-polyanion correlations -the so-called paddlewheel effect [31,32].While providing evidence of this effect, the simulations carried out in these works are clearly limited in the simulation times and system sizes they can achieve, potentially leading to unphysical outcomes.
In the last decade, the advent of machine learning has allowed the construction of interatomic potentials possessing quantum mechanical accuracy at a cost that is only marginally higher than that of classical force fields [33][34][35][36][37][38][39][40][41][42][43][44][45].Machine learning potentials (MLP) rely on the construction of physically-motivated representations to predict a given target property.In particular, representations of atomic configurations should preserve key physical symmetries: global translational and rotational invariance, as well as invariance with respect to the permutation of atoms of the same chemical species [46].Among the numerous potential representations, the The γ structure has all tetrahedra aligned along the [100] direction for all this set of planes, while the β structure has the tetrahedra aligned along both the [100] direction and the [ 100] across each plane.Finally, the α phase has a staggered ordering: the tetrahedra are aligned along [100]/ [ 100] for the planes that are numbered with even/odd numbers.
Smooth Overlap of Atomic Positions (SOAP) [47] used in combination with appropriate regression schemes has facilitated the development of ML potentials for simulating a variety of materials properties via extensive finitetemperature thermodynamic sampling [25,26,[48][49][50][51][52][53][54][55][56][57][58] Notable examples of the use of MLPs to study the ionic conductivity in solid-state electrolytes are the Gaussian Approximation Potential (GAP) for lithium thiophosphate developed by Staacke et al. [26] and the Deep Neural Network (DNN) for Li 10 GeP 2 S 12 -type compounds developed by Huang and coworkers [59].Both these studies were able to characterize the diffusion properties of their respective target compounds, and overcome some known limitations of AIMD, namely the small size and the short simulation times that are accessible by this type of modelling.Despite these important breakthroughs, two main aspects are still missing to provide a comprehensive study of transport properties in this class of materials.First, the accuracy of the aforementioned potentials is limited to the GGA level of theory, due to the choice of the reference DFT functional (PBE and PBEsol) for the calculation of the training set structures.While this is a standard choice for performing first-principles calculations in solids, the relatively small number of reference single-point calculations (usually a few thousands) that are needed to reach the desired target ML accuracy enables the use of more accurate references, like meta-GGA and hybrid functionals [60,61].To our knowledge, no systematic study comparing different DFT references exists up to date for this class of materials.Secondly, these studies neglect the contribution of inter-ionic correlations to the electrical conductivity, and its relation with polyanion rotations.
In this work, we train three MLPs to investigate the physical mechanisms of charge transport in Li 3 PS 4 and their effect on the electrical conductivity in its stable polymorphs.Each potential is trained over datasets computed at a different level of theory: GGA, metaGGA and hybrid functionals.In particular, we use the Perdew-Burke-Ernzerhof functional revisited for solids (PBEsol) [29,62], the regularized version of the strongly constrained and appropriately normed (r 2 SCAN) functional [63] and the PBE0 functional [64].We explore the temperature dependence of the ionic conductivity of Li 3 PS 4 showing that different functionals predict different critical temperatures for the onset of the conductive regime, which is roughly associated with the onset of a structural phase transition.We also elucidate the importance of including the effects of the interionic correlation in the conductivity, by computing it with the full Green-Kubo (GK) theory of linear response [65,66], instead of the Nernst-Einstein approximation commonly employed in the literature.Overall, we find that the PBE0 functional gives the best quantitative agreement with existing experimental measurements of the ionic conductivity of β-Li 3 PS 4 .Furthermore, we relate the onset of the superionic phase of the Li 3 PS 4 compound with the PS 4 flipping dynamics and find that discrete P-S flips induce a structural phase transition from the non-conductive γ to a mixture of the β and α structures, that cannot be fully resolved at the size and time scale of these simulations.This structural change determines a drastic decrease of the slope of the Arrhenius curve and thus a significant reduction of the activation energy of Li-ion diffusion (by a factor of six compared to the γ phase).Finally, we detect a second transition to a disordered phase with freely rotating polyanions at even higher temperatures that we attribute to the melting of the PS 3− 4 sublattice.Both the transition to the conductive phase and the Li 3 PS 4 melting appear as peaks of the heat capacity and are thus associated with separate first-order phase transitions of this material.

A. Training set construction and validation of the ML models
We construct the training set for the ML models in an iterative fashion.A starting dataset of structures is generated by running NVT Car-Parrinello Molecular Dynamics [67] with the PBEsol functional [62] as provided by the Quantum ESPRESSO package [68][69][70][71], for a set of selected temperatures (250 K, 500 K and 1000 K) and volumes.This initial dataset is then computed via more converged single-point DFT calculations with the PBEsol functional.Further details can be found in the Supporting Information.
As a next step, we fit a preliminary MLP on this dataset and run finite-temperature Molecular Dynamics (MD) with i-PI [72] over the entire temperature range of interest (between 200 and 1000 K).Among the uncorrelated structures generated in the resulting trajectories, only a subset consisting of the most diverse ones according to the Farthest-Point Sampling (FPS) method [73] is selected and recomputed using DFT.This active-learning loop, consisting of the regression of an MLP, MD simu-lations and re-calculation of a set of structures by DFT, allows us to extend the dataset until the model is deemed sufficiently accurate and robust.In order to generate datasets for the ML-SCAN and ML-PBE0 models, we select via FPS a subset of snapshots out of the whole PBEsol dataset and we use a two-level machine learning (2LML) scheme [74] to train accurate potentials from a minimal number of expensive r 2 SCAN or PBE0 calculations.The 2LML is a specific case of the general multilevel machine learning scheme, and consists in: training a ML model on the large PBEsol dataset, then computing energy and forces at the r 2 SCAN (PBE0) level and finally training a new ML potential on the difference between the ML-PBEsol predicted energies and forces and the true r 2 SCAN (PBE0) references [74,75].The final datasets consists of 2400 structures for the ML-PBEsol model, 740 for the ML-r 2 SCAN model and 790 structures ML-PBE0 model.Within these datasets, a subset of 100 randomly selected structures is used as a test set for the ML-PBEsol model, while a subset of 40 structures is used for both the ML-r 2 SCAN and the ML-PBEsol model.PBEsol calculations are performed with Quantum ESPRESSO, while r 2 SCAN and PBE0 calculations are performed with VASP [76][77][78].The training of all the models is performed targeting the cohesive energies to avoid offset issues induced by different pseudopotentials.Figure 2 shows the parity plot for the forces of the three models over their respective test sets.Table I contains the root-mean-square-errors (RMSE) for all models, showing that our model can achieve errors similar (or better) than those obtained in other similar works [25,26,59].The learning curves for each of these models are reported in the Supporting Information (Sec.I).The Supporting Information also reports results from kernel principal component analysis [79] (Sec.II) and the newly introduced local prediction rigidity [80] (Sec.III) to check the distribution of the environments in our dataset and along the MD trajectories, and to verify that our training set can reliably represent the complex local environments that occur during PS 4 flips.Dynamical properties, like the mean square displacement and atomic diffusivity of Li ions, also appear to be properly reproduced by our ML potentials, as we have directly tested via MD simulations of the α phase at high temperature, obtained with PBEsol ab initio potential and with its corresponding ML model, as reported in Sec.IV of the Supporting Information.

B. Collective variables for Li3PS4
The three main polymorphs of Li 3 PS 4 are differentiated by the relative orientation of the PS 4 tetrahedra.In order to distinguish them and identify phase transitions in MD simulations, we construct two collective variables (CVs), based on the alignment, along the [100] direction, of the tetrahedra of the (010) planes (see Fig. 1).To this aim, we first compute the polar angle θ SP , spanned by the vector r SP ≡ r S − r P that connects, for any PS 4 group, a given S atom with the central P atom, with respect to the x-axis shown in Fig. 5: One can then define an order parameter that measures the average alignment of PS 4 tetrahedra within each (010) plane, as follows: where C = 1, . . ., 6 labels the (010) planes as shown in Fig. 1, N C P = 16 is the number of P atoms in each plane, and ⟨s, p⟩ represents the S atoms that are first nearest neighbours of P atom p.In other words, the outer sum runs over P atoms that belong to a given crystallographic plane, while the inner sum runs on the four S atoms that belong to the tetrahedron centered around atom p. Whenever one tetrahedron is perfectly aligned along ±x, the cosine of one P-S angle is close to ±1, while the remaining three have a value of cos(109.5 • ) ≈ −0.3338 or cos(70.5 • ) ≈ +0.3338 for +x and −x respectively, due to tetrahedral symmetry.We thus raise the cosine to the fifth power in Eq. ( 2), so that the result of the inner sum will be approximately equal to ±1.While any odd power greater that 1 would serve this scope, 1 the fifth power gives the best results on preliminary tests.
Since the final aim is to define a global measure of the relative alignment across planes, so as to capture the staggered ordering of the α structure, we construct two intermediate order parameters s even and s odd , by averaging s C for even and odd values of C. We finally construct the following pair of collective variables: 1 This observation stems from the symmetry of the tetrahedra, ⟨s,p⟩ cos (θsp) = 0 In Sec.III B we demonstrate that s 1 completely distinguishes the three polymorphs, as it takes on values close to +1 when the structure is similar to the perfect γ, −1 when it is similar to α and 0 when the structure is close to β or fully disordered.s 2 instead measures a global alignment of the entire structure and does not resolve the β and α structure, as they both contain mixed PS 4 orientations in equal proportions.Still, this CV is meaningful when combined with s 1 as it carries information about the relative number of tetrahedra that are aligned along the positive and negative direction of the x-axis.For instance, it can distinguish between two perfect γ structures that are mirror-symmetric with respect to a (100) plane.

C. Green-Kubo theory
The Green-Kubo (GK) theory of linear response [65,66] provides a rigorous and elegant framework to compute transport coefficients of extended systems in terms of the stationary time series of suitable fluxes evaluated at thermal equilibrium with MD.For an isotropic system of N interacting particles, the GK expression for the electrical conductivity reads: where k B is the Boltzmann constant, T the temperature and Γ t indicates the time evolution of a point in phase space from the initial condition Γ 0 , over which the average ⟨•⟩ is performed.J q is the charge flux, that can be easily computed from MD, knowing the velocities of the atoms, v i , and their charges, q i : Here, the sum runs over all the atoms, e is the electron charge, and the q i are equal to the nominal oxidation number of the atoms [81]: in the absence of electronic conductivity due to conduction electrons or polaronic states, the overall electrical conductivity coincides with that obtained from Eq. ( 6) using integer, time-independent ionic charges [82].
A commonly used approximation of Eq. ( 5) is the Nernst-Einstein (NE) equation: where D Li and N Li represent the diffusion coefficient and the total number of the lithium atoms respectively.Equation ( 7) is widely used in practice to estimate the ionic conductivity, due to the high statistical accuracy with which atomic diffusion coefficients can be computed from numerical simulations [83].Nevertheless, its application to solid-state-electrolytes (SSEs) is burdened by systematic errors [84]: in fact, the large interatomic dynamical correlations, both between carriers (Li + ) and the solid matrix (PS 4 3− ) and among the carriers themselves, which is typical in systems with a high carrier concentration like SSEs, is completely neglected by Eq. (7).The discrepancy between σ and σ NE can be quantified by the Haven ratio [84][85][86]: We redirect the reader to Sec.III D for a thorough comparison between σ and σ NE in the different polymorphs of Li 3 PS 4 .
From a methodological standpoint, Eq. ( 5) can be expressed in an equivalent formulation, called the Helfand-Einstein (HE) formula, which exhibits better statistical behaviour [87]: The Li diffusivity appearing in Eq. ( 7) is obtained from the asymptotic slope of the mean square displacement of the Li ions: In this case, care has to be taken to compute D Li in the reference frame where the solid matrix is fixed to avoid non-physical contribution to the calculations of the electrical conductivity.These spurious effects arise for simulations run in the barycenter reference frame, where the position of the center of mass of the entire system is fixed.In practice, the difference between D Li computed in these two reference frames vanishes when the box size is increased [88].Due to its very general formulation, the GK expression of Eq. ( 5) can be used to investigate, with minimal variations, other characteristic properties of Li 3 PS 4 .In Sec.III C we will characterize the rotational properties of the PS 4 polyanions at high temperature, by computing a rotational diffusion coefficient as follows: In this equation j runs over every P-S bond of each PS 4 tetrahedron and ω j (t) represents the time series of its angular velocity: where r j P,S (t) and v j P,S (t) are the positions and the velocities of the P and the S atom belonging the j-th bond.

D. ML-MD computational details
We use the MLPs constructed in Sec.II A to investigate the physics of Li 3 PS 4 via constant-temperature MD simulations using a combination of i-PI [72,89], lammps [90,91] and librascal [92].In order to simulate the phase transitions and the charge transport in Li 3 PS 4 , we perform MD simulations in the NpT ensemble of a quasicubic 768-atom cell in all stable α, β and γ phases with a constant isotropic pressure of p = 0 atm for a set of temperatures between 200 K and 1000 K.The system's center of mass is kept fixed during the simulations.A generalized Langevin equation (GLE) thermostat [93,94] is used to equilibrate the cell volume, while a stochastic velocity rescaling (SVR) thermostat [95] is used to thermalize the velocity distribution of the atoms without affecting significantly the dynamical properties.The characteristic times of the barostat, the SVR thermostat and the MD timestep are set to 1 ps, 10 fs and 2 fs, respectively.We run these simulations long enough to ensure statistical convergence of the ionic conductivity (see Sec. II C).Specifically, we run the weakly conductive simulations of the γ phase for ∼ 6 ns, the β phase for ∼ 4 ns the α phase for ∼ 2 ns.As discussed in Sec.III B, this setup ensures that the system can sample configurations within the range of stability of each phase, without explicitly requiring a quantitative prediction of the temperaturedependent phase diagram of Li 3 PS 4 .A validation of the setup via a heat-quench simulation in the NST ensemble is found in the Supporting Information (Sec.VII).

III. RESULTS
In this section, we compare the results obtained with the PBEsol, PBE0 and r 2 SCAN functional and the corresponding ML models, as described in Sec.II.In subsection III A, we compute the electronic band structure of the β polymorph using DFT and show that the band gap predicted by the PBE0 functional is in good agreement with experiments, while PBEsol and r 2 SCAN considerably underestimate it.In the following subsections, we investigate the finite-temperature predictions of the ML models.First, in subsection III B, we analyse the MD simulations of the Li 3 PS 4 polymorphs using the collective variable introduced in Sec.II B and discuss the onset of a phase transition from the γ structure to a structure with a mixed α and β arrangement by increasing the temperature.In subsection III C, we investigate the rotational dynamics of the PS 4 tetrahedra and relate the occurrence of phase transitions with the thermal activation of correlated polyanion flips.Furthermore, we detect the presence of a high-temperature liquid phase characterized by freely rotating polyanions.In subsection III D, we compute the ionic conductivity from MD NpT simulations using both the NE and the HE expressions introduced in Sec.II C and the Haven ratio of Li 3 PS 4 .Finally, in subsection III E, we discuss the role of spatial correlations and their effect on the calculated ionic conductivity.
A. Electronic band structure The Generalized Gradient Approximation (GGA) functionals offer a good compromise between accuracy and computational efficiency, making them a practical choice for a broad range of materials and systems.It is known, however, that they often fall short when it comes to accurately characterizing critical electronic properties, such as the electronic band gap, which is frequently underestimated in GGA, and the density of states [98][99][100].In order to solve this problem, different functionals have been developed, such as meta-GGA [60] and hybrid functionals [61], that offer more accurate approximations of the exchange-correlation functional and are able to better describe long-range electron-electron interactions.These new functionals have enabled in more accurate predictions of electronic properties in a variety of different materials [101][102][103][104][105][106] and notably solid-state electrolytes [107][108][109].In Fig. 3 we compare the band structure and the density of states (DOS) for the β phase computed with the PBEsol (GGA), r2 SCAN (meta-GGA) and PBE0 (hybrid) functionals 2 .Since in the β phase the Wyckoff sites of the Li atoms have partial occupations, we perform the calculation using one (B3C1 [108]) of the known configurations with minimum energy, since the electronic bands and the gap are only weakly dependent on this choice [108].Furthermore, in Table II we compare the band gaps predicted by the different functionals in the γ and β phases and recent experimental measurements.We note that both PBEsol and r 2 SCAN considerably underestimate the electronic band gap, while PBE0 shows a remarkably good agreement3 , thus further motivating the use of this functional as a reference for the training of a dedicated ML model.

B. Structural phase transitions
We use the pair of CVs s 1 and s 2 , introduced in Sec.II B, to investigate the presence of structural phase transitions appearing in the MD trajectories.As anticipated, s 1 characterizes the mutual orientation of adjacent (010) planes (i.e., even and odd numbers in Fig. 1) along the [100] direction.Consequently, s 1 = −1 for the α phase, where the PS 4 orientation of adjacent planes is antiparallel, s 1 = 0 for the β phase (each plane has no net orientation of PS 4 units), and s 1 = +1 for the γ phase, where adjacent planes share the same orientation of PS 4 tetrahedra along the positive x axis.Conversely, s 2 measures the global orientation of PS 4 tetrahedra and vanishes for both the α and β phases, while is equal to 1 for the γ phase.Figure 4 displays, with red dots, the evolution of the CVs across a set of MD simulations run with the ML-PBE0 model at different temperatures T , and initialized in the γ phase.The green markers of three different shades represent reference points sampled from MD simulations in the α, β and γ phases below T c = 750 K, and are used as a guide to interpret the T -dependent results.For T < T c , the red dots are all concentrated in a region around (s 1 , s 2 ) = (1, 1), which is typical of the pure γ phase, the small deviations being due to the thermal motion of the atoms only.As T is raised above ≈ T c , a structural transition occurs, and the CVs approach a region in between the α and β phases.Although the total timescale and size of the simulations are not sufficient to allow for a complete transition of the γ phase to a specific polymorph, but only to intermediate configurations, our MD simulations capture the microscopic driving mechanism, i.e., the onset of concerted reorientations of PS 4 tetrahedra.Figure 5 shows a typical example of this phenomenon, occurring during a ≈ 15-ps segment of a MD trajectory: the starting (ending) configuration is depicted in red (blue).The color fades from red to blue (with an RWB scheme) in a continuous manner as the transition occurs.The trajectories of the S atoms lying at the vertices of the PS 4 tetrahedra clearly indicate that the reorientation of the entire row occurs coherently, and not as a collection of individual, decorrelated flips.Figure 5 also shows that the transition is purely orientational, and does not occur through the hopping of S anions between adjacent PS 4 groups.Additional information and a comparison of this mechanism with an heat-quench simulation showing a similar behavior can be found in the Supporting Information (Secs.VII and VIII).

C. PS4 rotational dynamics and heat capacity
Further insights into the rotational reorientation of PS 4 planes and their relation with structural transitions can be obtained by a direct investigation of the rotational dynamics of PS 4 tetrahedra, and specifically those that form [100] rows.Figure 6 shows the dynamics of the polar angles θ SP , as defined in Eq. ( 1), for a set of four tetrahedra (row panels) that belong to the same [100] row in NpT simulations run with the ML-PBE0 model4 .We also compare three trajectories initialized in the γ phase and equilibrated at T = 725 K, i.e. just below T c = 750 K (left column); at T = T c (central column); and above melting, at T = 900 K (right column).The four lines in the plots correspond to the dynamics of the four bonds that constitute each PS 4 tetrahedron.
At 725 K, only small angular fluctuations occur, with no reorientation of the tetrahedra.The average angles of the P-S bonds define the mean orientation of each tetrahedron: as one of the P-S bonds always has an average value that is close to 0, the orientation is along the positive x-axis (+x) during the entire simulation time, which is typical of the γ phase.Instead, the angles of the other three bonds oscillate around θ 0 = 109.5• (black dashed FIG. 6. Sketch of the rotational dynamics of the PS4 groups: low temperature (panel a), where only small librations with respect to the initial configuration occur, at intermediate temperature (panel b) where PS4 flips determine the structural transition observed in Fig. 4 and at high temperature (panel c) where the system is melted.The lower plots show the time evolution of the polar angle θSP as defined in Eq. ( 1) (angle with respect to the x-axis) for a set of three distinct tetrahedra forming a [100] row.Each panel represents the dynamics of the four PS bonds forming each tetrahedron.Rows correspond to different tetrahedra, while columns correspond to different NpT trajectories at T = 725 K, 750 K and 900 K. Horizontal dashed black lines indicate the position of the ideal tetrahedral angles at 70.5 • and 109.5 • , while grey dotted lines mark the extremes of the domain of θSP (i.e.0 • and 180 • ).line), corresponding to the P-S bond angles in the perfect tetrahedral geometry.The jump observed in the central panel of the left column corresponds to a rotation of the tetrahedron, that does not involve a re-orientation towards the negative x-axis.Instead, at the transition temperature (T c = 750 K), simultaneous flips of one P-S bond from 0 to 70.5 • and another P-S bond from 109 • to 180 • correspond to the reorientations of the tetrahedra from +x to −x, consis-tently with the mechanism shown Fig. 5. Notably, these flips occur at the same instants of time for every tetrahedron in a row (see, e.g., central column of Fig. 6 between 1 and 2.5 ns), confirming that they are highly correlated across [100] rows.This effect is at the basis of the phase transition observed in Fig. 4 as it modifies the relative orientation of tetrahedra across (010) crystallographic planes.It is crucial to note, in this respect, that the spatial correlation of these flips at the transition ex- tends up to the edge of the simulation box, potentially leading to finite-size effects.Specifically, we expect this transition to manifest in larger boxes by nucleation of ordered clusters with opposite orientation, as a result of the formation of defects with sudden changes of the PS 4 orientation.We also note from panels b) of Fig. 6 that the time lag between subsequent flips is of the order of 1 ns.As we will see in Sec.III D, the presence of this long time scale will be important to elucidate the mechanism of charge transport in this material.At 900 K, the dynamics of the tetrahedra changes dramatically and a second phase transition occurs.In particular, all PS bonds across every tetrahedron span the entire range of angles between 0 • and 180 • .This suggests that, unlike the phases observed at lower temperatures, the tetrahedra are freely rotating in the simulation box.A direct inspection of the simulations indicates that this behaviour is accompanied by a melting of the system into a mixture of Li + cations and PS 4 3− anions (see also the P-P and P-S radial distribution functions shown in the Supporting Information, Fig. S9).
The onset of these two phase transitions-from the nonconductive γ phase to the superionic β/α hybrid phase and the melting of Li 3 PS 4 -can be quantitatively investigated for all ML models by computing the temperature dependence of a set of relevant quantities.Panel a) of Fig. 7 shows, for each of the ML models, the isobaric heat capacity C p (T ) computed from the finite-difference derivative, with respect to T , of the mean enthalpy collected in the NpT simulations.Panel b) shows the temperature dependence of the mean squared fluctuation of the polar angle θ SP as defined in Eq. 1 and further averaged over every P-S bond.Panel c) shows the PS 4 rotational diffusion coefficient D ω , defined in Eq. 11, and panel d) the linear diffusion coefficient of the P atoms, D P .
C p (T ) displays two distinct peaks, characteristic of the phase transitions observed in Fig. 6, while the associated critical temperatures depend on the specific functional.We can characterize more clearly the position of these peaks by analysing their relation with the microscopic quantities shown in panels b), c) and d).Since the first transition is associated with discrete PS 4 flips, the increase of C p is accompanied by a sudden change of slope of the angular fluctuations.Conversely, the transition to the molten phase occurs with a dramatic increase of both D ω and D P by one and two orders of magnitude, respectively.In other words, the action of thermal fluctuations at this high temperature destroys the periodic arrangement of P atoms, while the tetrahedra are still intact and can freely rotate by a rate given by D ω .Notably, the transition temperature to the molten salt as predicted by the ML-PBE0 model is in agreement with a previous experimental measurement of the binary phase diagram of β-Li 3 PS 4 -Li 4 GeS 4 solid solutions obtained through differential thermal analysis [110].More specifically, the transition point upon heating is reported to be 600 • C when the concentration of β-Li 3 PS 4 is equal to 98% (P-rich regime), which is compatible with our PBE0 estimate.
We stress that the angular deviations of panel b) can be defined only with respect to a local equilibrium for each P-S bond and are thus meaningful only in the low-T phase.Conversely, the diffusion coefficients of panels c) and d) are physically meaningful when the simulations sample sufficiently many configurations with displaced P atoms and rotated PS 4 anions.They are thus not well defined if the MD simulations are not fully ergodic [111].In Fig. 7, we display each of the quantities with solid lines in the regions where they are well defined, otherwise they are shown with dotted lines.The phase transitions investigated so far have strong implications on electrical conduction, as we describe in the following subsection.

D. Ionic conductivity and Haven ratio
The ionic conductivity, σ, is the crucial property to identify promising solid-state electrolytes.As discussed in Sec.II C, the GK theory of linear response, in its HE formulation, gives us an efficient and statistically robust method to obtain an estimate of σ from equilibrium MD simulations at any target temperature.
The upper panels of Fig. 8 shows the temperature dependence of the ionic conductivity at zero pressure for a set of NpT simulations that start from the ideal α, β and γ polymorph.Imperfect ergodicity, and the constraints on cell shape, make simulations dependent on the initial conditions.Even though simulations can only be considered converged within the stability range of each phase and target functional, results outside this range still report useful information about their behavior when metastable.Thus, results are shown for every target functional and for all temperatures where the ionic conductivity is non-zero within the errorbars.
All the ML models predict the γ/α phase to be the least/most conductive, while the β phase has intermediate values of σ over a wide temperature window (up to 600-800 K depending on the model).This result is in agreement with previous computational studies [26,112] and experimental measurements of the ionic conductivity [20,28] on the known crystalline phases of Li 3 PS 4 .The negative slope of the profiles of σ with respect to the inverse temperature is typical of the Arrhenius plots.In fact, since lithium-ion diffusion occurs via thermal activation, one can mathematically relates the Nernst-Einstein conductivity σ NE with the activation energy of the Li-hopping process.In particular, higher negative slopes correspond to higher activation energies.The γ (α) phase is thus not only the least (most) conductive phase but also the one that has the largest (smallest) activation barrier for lithium diffusion.
The high-temperature ends of the conductivity profiles of Fig. 8 (for T < 670 K for ML-PBEsol, T < 880 K for ML-r 2 SCAN and T < 770 K for ML-PBE0) give us additional insights.While the curves related to the three phases clearly span different conductivity ranges and show different slopes at low temperatures, this difference is not noticeable any more at high temperatures and the conductivity and the activation energies all approach the values of the α phase, with a characteristic kink that is most visible for simulations initialized in the γ phase.The critical temperature at which this kink occurs depends on the reference functional.This is a clear effect of the structural transition studied in Sec.III B and further investigated in Sec.III C by an analysis of the PS 4 rotational dynamics.In other words, the PS 4 flips that induce the transition from the γ polymorph to the partly-ordered β-α structure are responsible for the changes of the ionic conductivity due to the larger availability of hopping sites for Li-ions [20], as well as a reduction of the Li-hopping activation barrier.This conclusion is also highlighted by simulations of the β and α phases at low temperatures, that show large σ although PS 4 flips are suppressed and librations are weak.To quantify the reduction of the activation energy due to the structural transition, in Table III, we report the activation energies that are fit to the Nernst-Einstein conductivities below and above the transition temperatures of each model (see Sec. X of the Supporting Information for additional details).The effect of the transition is remarkable, as we observe a reduction of the activation energies of up to a factor of 6 depending on the reference DFT functional.Furthermore, the activation energies are very small above the transition, with values ranging between 0.25 and 0.32 eV.These values are very close to the values observed for the α phase and smaller than the β phase (see Table IV), indicating a superionic behavior.In contrast, we note that the transition to the molten salt, that we observed for T > 800 K in Sec.III C has practically no effect on the conductivity profiles.
This analysis, combined with the results of Secs.III B and III C, also allows us to rule out any paddle-wheel effect, whereby PS 4 motion is time-correlated with Lihopping and increases Li-ion diffusion.In fact, due to the different rates of PS 4 flipping (≈ one every ns) and Li-ion hopping (≈ one every ps) even at large temperature right below melting, the two mechanisms cannot be related.This point is strengthened in Fig. S12 of the Supporting Information, showing the fast and linear increase of the mean square displacement of lithium ions in the simulation at 750 K of Fig. 6 over a 30-ps time window.Instead, the interaction between Li-ion diffusion and PS 4 libration, which share the same time scale, may be important, as we have directly inspected in Sec.III E by analyizing the local contributions to σ stemming from the interaction between Li and S ions.
The lower panels of Fig. 8 show the temperature dependence of the Haven ratio, H R , computed using Eq. ( 8).As anticipated in Sec.II C, this coefficient quantifies the discrepancy between σ and σ NE .H R < 1 for almost every temperature and only for the γ phase it approaches one at low temperatures, where the system is weakly conductive.This indicates the presence of inter-atomic correlations in the ionic conductivity, both between carriers (Li + ) and the solid matrix (PS 4 ), that cannot be captured by the NE approach (see Eq. ( 7)), as the latter only estimates the conductivity based on the self-diffusion of the lithium ions.This effect is most pronounced in the α phase, where H R ≈ 0.4 below melting, meaning that the NE estimate underestimates the GK conductivity by more than a factor of two.At high temperatures, the Haven ratio slightly increases, indicating that the material becomes disordered.This might be a result of the phase transformations of the solid matrix, which we expect to weaken the interionic correlations.Still, the Haven ratio never exceeds 0.8, even at 1000 K for any of the ML models studied.
In Table IV we quantitatively compare the activation energies and the ionic conductivities predicted by our ML models with recent experimental measurements [28,96] III.Activation energies for Li-ion diffusion for the simulations initialized in the γ phase for temperatures below the phase transition temperature (T > Tc) observed in Sec.III C and above Tc, see the SI for details on the computation of the activation energies.The last column is the ratio between the EA(T < Tc) and EA(T > Tc) and quantifies the reduction of the Li-diffusion activation energy due to the phase transition.
and with computational studies based on AIMD at the PBE level [22,114,115].The first two columns of Table IV compare the activation energies in the β and α phase, while the remaining three report the estimates of the ionic conductivity of the β phase at 298 K 5 and 500 K and the ionic conductivity of the α phase at 298 5 This value is extrapolated from a fit of the temperature profiles of the β structure (see Fig. 8) for each ML model, see SI for details on the computation of the activation energies.
K. The ML-PBEsol model predicts activation energies in agreement with previous AIMD for both the β and the α structure.In the case of the β phase, the MLr 2 SCAN and ML-PBE0 predict an activation energy that is greater than the one computed with the ML-PBEsol model, but overall the values are close to the experimental results.For the α phase the activation energies are particularly close to the experiment.These last results are remarkably good in particular when comparing them with the prediction of the classical empirical potential recently introduced by Forrester et al [23].For this model, the activation energy for stoichiometric α-Li 3 PS 4 is much larger (0.40 eV) than the experimental result and similar to the value of the β phase.The empirical potential also predicts the α phase to be only slightly more conductive than the β phase for all temperatures.
In conclusion, our analysis shows that the ML potentials are the only possible solution to accurately predict the properties of Li 3 PS 4 , given the unreliability of empirical potentials and the prohibitive cost of ab initio simulations, in particular at PBE0 level.TABLE IV.Comparison between the predicted activation energies and conductivities (computed at 500 K and extrapolated at room temperature; see SI for details) of the β and α phases with both experimental references [28,96] and previous ab initio MD studies [22,114,115].

E. Spatial correlations
In this last section, we investigate the role of local correlations in determining the full ionic conductivity by computing the spatial dependence of the integral of the partial cross-correlation functions, I LiA , between lithium atoms and other atomic species A = Li, P, S.However, before we start this analysis, it is necessary to make a few technical considerations.While the total conductivity does not depend on the frame of reference, due to the charge neutrality of the simulation cell, the value of any partial correlation does depend on it, as we show in the SI.For instance, in the reference frame of the matrix, all the Green-Kubo integrals of the partial correlation between Li and the other species (P and S) vanish.In contrast, in the reference frame of the center of mass of the entire system the solid matrix recoils due to the diffusion of the center of mass of Li atoms (see SI).To carry out our analysis, we choose the reference frame where the center of mass of the PS solid matrix remains stationary, e.g.
i∈P,S v i = 0 at every instant.This choice is motivated by the fact that this reference frame is in principle the same in which the lithium diffusivity, entering the NE relation, Eq. ( 7), should be computed.Moreover, this is the most natural choice for a battery, since in this reference frame the solid matrix of the battery is not moving.
To study the spatial dependence of the integral of the correlation functions we perform a Gaussian kernel density estimation (KDE) [116] of the correlation functions: where i runs over all Li atoms and j runs over all the atoms of type A; w ij (t, r) is a Gaussian weight with width ς, which we set to 0.33 Å; wij (r) is the time average of w ij (r) over the whole trajectory.The connection between Eq. ( 13) and the ionic conductivity of Eq. ( 5) is readily established.In fact, within the KDE, where g LiA (r) is the radial distribution function between the Li and the species A. Therefore, the integral over r of I LiLi (r) plays the role of the correction to the NE relation that is needed to recover the full GK conductivity: Fig. 9 shows I LiA (r) for different values of r.The most relevant correlations come from the Li-Li and Li-S interactions, while the correlations with the P are compatible with zero for any r.As expected, for large r all the correlations become compatible with zero.The oscillations in the correlations functions follow, as expected, the peaks in the radial distribution function g LiA (r) (lower panel), and show that only the first shell of Li atoms and the first two of S atoms are important for the conductivity.The correlation of the first shell of Li is strongly positive in agreement with Ref. [117], where it was shown that the hopping of Li atom in one direction facilitates the movement of the other Li in the same direction.

IV. CONCLUSIONS
In this work, we have presented a computational study of lithium ortho-thiophosphate via multiple machine learning models, targeting DFT references of increasing accuracy, and elucidating the critical role of PS 4 flips and phase transitions in determining the observed superionic behavior of Li 3 PS 4 .We find that all the ML models predict two distinct phase transitions.First, we observe a transition from the γ to a partly ordered phase with both β and α alignments, that cannot be fully resolved at the time and length scale of the MD simulations.While this limitation does not allow to provide a prediction of the full phase diagram of Li 3 PS 4 , we identify the presence of collective PS 4 flips as the driving mechanism.Secondly, we observe the melting of the system into its constituent Li + cations and PS 3− 4 polyanions at elevated temperatures (>800 K).Both these transitions are associated with drastic changes in the rotational dynamics of the PS 4 groups as a function of the temperature and appear as distinct peaks of the heat capacity.
We also compute the ionic conductivity of Li 3 PS 4 in all its stable polymorphs and elucidate the importance of including the effects of interatomic correlations, by computing it with the full Green-Kubo theory of linear response and with the Nernst-Einstein approximation.We find that the interionic correlations account for considerable deviations between the NE and the GK estimates, as quantified by a Haven ratio that is smaller than one in every polymorph at all temperatures, except for the weakly conductive γ structure at low temperatures.Notably, the Haven ratio reaches values of 0.4 in the highly conductive α-phase, suggesting that a pure NE approach can result in the underestimation of the conductivity by more than a factor of two.From a spatially-resolved analysis performed in the reference frame of the solid matrix, we find that most of these interionic correlations come from the first shell of Li-Li neighbors, thus indicating that a concerted Li-ion hopping is a key aspect of charge transport in this material, in agreement with Ref. [117].
Finally, we investigate how the observed phase transitions of Li 3 PS 4 affect the ionic conductivity.We find that the occurrence of correlated PS 4 flips results in a dramatic decrease of the activation energy (up to a factor of six) when the system transitions from the γ to a mixed β-α phase.Furthermore, we show that subsequent PS 4 flips occur at the time scale of nanoseconds, that is much larger than the typical time laps between two subsequent lithium ion hoppings, in agreement with Refs.[117].We thus conclude that the sudden change in the PES of the lithium ions that is due to the rearrangement of the PS 4 tetrahedra is the physical mechanism for the observed superionic behavior of Li 3 PS 4 .Crucially, this mechanism is fundamentally different from the one proposed in Ref. [31].There, a characteristic paddlewheel effect was observed in AIMD-PBE simulations at elevated temperatures and invoked to explain the neutron diffraction measurements showing a polyanion reorientational disorder.The AIMD simulations were however limited in size and thus showed much larger finite-size effects than the ones we observe in this work, including an artificial stabilization of the solid phase up to temperatures (1200 K) that are much larger than the nominal melting point of Li 3 PS 4 .We also stress that the paddlewheel effect itself was observed in AIMD simulations at these very high temperatures, thus likely making it an artifact of the small simulation box that was used there.The ML models that we present in this work overcome these limitations and offer a more natural interpretation of the experimental results.This finding also suggests additional directions of research in the quest for a promising solid electrolyte and potentially a way to design new target compounds.In particular, we expect that tentative modifications of Li 3 PS 4 to stabilize its superionic phases at room temperature, by e.g.atomic substitution and amorphization, should be accompanied by a reduction of the polyanion rotational free energy barrier, that limits the spatial extension of the PS 4 fluctuations.Further developments of this work thus imply a detailed thermodynamic study of the phase transitions observed here and a comparison of multiple different SE compounds, with the aim of suggesting a target compound for experimental synthesis.
While providing these useful mechanistic insights, our ML models show a remarkable agreement with experiment in the prediction of a number of independent quantities.Specifically, the PBE0 functional provides the best agreement on the prediction of the electronic band gap, while the associated ML-PBE0 model reaches overall the best accuracy on the prediction of lithium activation energies in the β and α phase, the ionic conductivities at 298 and 500 K and the melting temperature.In particular, our results proved to be much more accurate than empirical potentials [23], that are often used to overcome the high computational cost of AIMD.These results and the observed dependence of the finite-temperature predictions of the ML models on the DFT reference indicate the necessity of using more accurate functionals for the description of transport properties in solid electrolytes, than state-of-the-art GGA functionals.Machine learning becomes then a necessary step in modelling this class of materials, as ab-initio studies with the PBE0 functional are far beyond reach because of their very high computational cost.
In conclusion, we have shown how the use of machine learning potentials for a prototypical solid electrolyte captures the mechanisms of the transition to its supe-rionic phase and the quantitative values of the ionic conductivity, while also allowing us to investigate the role of inter-ionic correlations.This work thus opens up a new frontier in the exploration of superionic materials as it allows their large-scale simulations at hybrid-DFT accuracy for hundreds of nanoseconds.This will offer crucial insights into the fundamental properties of solid electrolytes, as well as guidance for the experimental realization of new candidate compounds.

FIG. 1 .
FIG.1.Sketch of the α, β and γ phases of Li3PS4, showing the difference of the relative alignment of the PS4 tetrahedra along reference (010) crystallographic planes (here numbered between 1 and 6 for clarity).The γ structure has all tetrahedra aligned along the[100] direction for all this set of planes, while the β structure has the tetrahedra aligned along both the[100] direction and the[ 100] across each plane.Finally, the α phase has a staggered ordering: the tetrahedra are aligned along[100]/[ 100] for the planes that are numbered with even/odd numbers.

FIG. 4 .
FIG. 4. Evolution of the collective variables of γ−Li3PS4(red points) sampled over a set of MD trajectories generated with the ML-PBE0 model.Green markers of three different shades represent a sample of reference points extracted from all MD simulations in the α, β and γ phase below Tc = 750 K where no phase transitions are observed.The purple markers in the topmost panel indicate the CVs for the ideal crystalline structures.

FIG. 5 .
FIG.5.Transition corresponding to the re-orientation of one line of PS4 tetrahedra.Tetrahedra colored in red and blue correspond to snapshots of the solid matrix taken over the transition for one trajectory at 750 K starting from the γ structure and run with the ML-PBE0 model.The trajectory of two vertices of each tetrahedron is displayed with lines, that are colored with an RWB scheme and with a smoothening window of 2.5 ps.The red end of these lines corresponds to t0, while the blue end to t1 = t0 + 15 ps.

FIG. 7 .
FIG. 7. a) the heat capacity (Cp), b) the fluctuations of the P-S polar angle averaged over all bonds of every PS4 tetrahedron in the simulation box ( ⟨∆θ 2 ⟩), c) the PS4 rotational diffusion coefficient (Dω) and d) the linear diffusion coefficient of the P atoms (DP) as a function of the simulated temperature.Solid/dashed lines represent each quantity in the temperature regime where it is well/ill-defined.The dashed black line of panel a) represents the Dulong-Petit limit, 3kB.

FIG. 8 .
FIG. 8. Temperature dependence of σ and Haven Ratio.Panels a), b) and c) show a comparison between the ionic conductivities predicted by the ML models as a function of the inverse temperature, computed via the Green-Kubo relation.Panels d), e) and f) show the behaviour of the Haven ratio, HR = σNE/σGK, as a function of the inverse temperature.Error bars obtained from standard block analysis over eight blocks are displayed.Vertical dashed lines indicate the experimental stability boundaries for the three phases.
FIG. 9. (Upper panel) Integrals of the cross correlation functions as defined in Eq. (13), for different pairs of atomic species and as a function of R. All the velocities are computed in the system of reference of the center of mass of the solid matrix of the battery.The shaded area indicate the uncertainty on the mean value, obtained from block analysis on 10 blocks.(Lower panel) Radial distribution function between the Li and all the other atomic species.The data reported in the figure are obtained from 190 ps of a simulation of the α-phase at 650 K.

TABLE II
[97]ergy Band Gap with the different models for the γ and β phase.The PBE0 values are in best agreement with the reported experimental value from Ref.[97].