Computational Characterization of β-Li3PS4 Solid Electrolyte: From Bulk and Surfaces to Nanocrystals

The all-solid-state lithium-ion battery is a new class of batteries being developed following today’s demand for renewable energy storage, especially for electric cars. The key component of such batteries is the solid-state electrolyte, a technology that promises increased safety and energy density with respect to the traditional liquid electrolytes. In this view, β-Li3PS4 is emerging as a good solid-state electrolyte candidate due to its stability and ionic conductivity. Despite the number of recent studies on this material, there is still much to understand about its atomic structure, and in particular its surface, a topic that becomes of key relevance for ionic diffusion and chemical stability in grain borders and contact with the other device components. In this study, we performed a density functional study of the structural and electronic properties of β-Li3PS4 surfaces. Starting from the bulk, we first verified that the thermodynamically stable structure featured slight distortion to the structure. Then, the surfaces were cut along different crystallographic planes and compared with each other. The (100) surface is confirmed as the most stable at T = 298 K, closely followed by (011), (010), and (210). Finally, from the computed surface energies, the Wulff nanocrystals were obtained and it was verified that the growth along the (100) and (011) directions reasonably reproduces the shape of the experimentally observed nanocrystal. With this study, we demonstrate that there are other surfaces besides (100) that are stable and can form interfaces with other components of the battery as well as facilitate the Li-migration according to their porous structures.


Introduction
In the last years, there has been an increase in the need for more efficient batteries, mainly driven by the transition towards electric car mobility and renewable energy sources. Research leads thus towards more efficient, portable batteries with a large energy storage capacity. Safety and reliability have become, more than ever, key aspects.
In this, all-solid-state-lithium batteries (ASSLBs) have shown good results because, in addition to improved safety and lifetime, they promise a good energy density, packaging, and operating temperature range. The main change concerning traditional batteries is the introduction of one (or eventually more) solid-state electrolyte (SSE), which main classes of SSE were recently reviewed in [1]. In general, ideal SSEs feature good ionic conductivity, high chemical stability, large electrochemical window, and good mechanical strength. Moreover, they should be environmentally benign and have low production costs. However, according to Homma et al. [2], a good SSE system should avoid nonmetallic elements such as germanium and/or silicon. As a matter of fact, these elements might easily be reduced and oxidized during the electrochemical reaction, yielding higher electrochemical instability. This is an important condition for the solid electrolyte in solidstate batteries.

Materials and Methods
All the calculations were performed with the CRYSTAL17 package [7], within the Density Functional Theory (DFT) method as developed for periodic systems in a basis set of localized Gaussian-type atomic functions, centered on each atom. Computational parameters were chosen and optimized by performing a preliminary analysis on the β-Li 3 PS 4 bulk. In particular, the hybrid PBE0 functional [8,9] was adopted and the 6-11G [10], 86-311G* [11], and 85-211dG [12] basis sets were used to describe the Li, S, and P atoms, respectively.
The DFT integration was performed within a grid containing 99 radial and 1454 angular points, as specified by the XXLGRID keyword [7,13]. The accuracy of the truncation criteria for the bi-electronic integrals, Coulomb and HF exchange series, is controlled by a set of five thresholds for which the strict values of [8,8,8,8,16] were adopted. In the self-consistent field (SCF) procedure, the shrinking factor for both the diagonalization of the Fock matrix and the calculation of the energy was set to 4, corresponding to 4 independent k-points in the irreducible part of the Brillouin zone. The total and projected density of states (DOS) and the band structure were plotted using the same k-point sampling as in the SCF. The vibrational frequencies at the Γ point were computed within the harmonic approximation by diagonalizing the mass-weighted Hessian matrix [14,15].
Topological analysis of the electron density, ρ(r), was conducted according to the quantum theory of atoms in molecules and crystals (QTAIMC) as developed by Bader [16] and implemented in the CRYSTAL code [13,17], adopting the default values for all the computational parameters [18]. In the context of our investigation, it provided information on the bonding framework and helped highlight significant differences in the electronic structure of similar crystals. It consists of a few steps. First, the points where the gradient of the density vanishes, ∇ρ(r) = 0, i.e., the ρ(r) critical points, are found. Then, they are classified at local minima, maxima, and saddle points according to the sign of the second derivatives of the density [19]. Among these, the bond critical points (BCP) represent the saddle in the charge density, being a minimum along the atom-atom direction and a maximum in the two perpendicular directions. Interesting enough, chemical bonds can be classified according to the value of the Laplacian, ∇ 2 ρ(r), the potential energy density, V(r), the positive definite kinetic energy density, G(r), and the bond degree, H(r)/ρ(r), with H(r) = V(r) + G(r), as calculated in r = r BCP [20]. Put simply, interactions can be subdivided in (i) covalent, if the Laplacian is negative and H(r) and V(r)/G(r) > 2, i.e., if there is an excess of potential energy at the BCP; (ii) transit, if the Laplacian is positive, the bond degree is close to zero and 1 < V(r)/G(r) < 2; (iii) ionic and multipolar (i.e., hydrogen bonds and van der Waals interactions) characterized by positive Laplacian and H(r) and V(r)/G(r) < 1 due to the dominance of kinetic energy at the BCP [21,22].
As the β-Li 3 PS 4 bulk contains 4 units of Li 3 PS 4 in the reference cell (32 atoms in total), the number of layers of each surface was chosen to be multiple of 4. Therefore, starting from the minimal 32-atom unit, slabs with an increasing number of layers were considered. For each surface, different terminations along the ±z-direction were analyzed in order to consider different types of low coordination atoms and surface morphologies (see Table S4). The results for each surface, considering the different thicknesses and possible terminations, are discussed in Topic S2 of the Supplementary Information (SI) and shown in Figures  are the energies of the non-optimized surfaces with and without ghost functions. In E form , the basis set superposition error (BSSE) is incorporated, whose value is related to the atoms that belong to different structures approaching each other, interacting, and the overlap among their basis functions produces a not-genuine stabilization that is inversely proportional to the quality of the adopted basis set. In this work, a posteriori counterpoise (CP) method [23] was applied.
The thermodynamic stability was evaluated by computing the Gibbs free energy at where E elec is the electronic energy, E 0 the zero-point vibrational energy, E T vib is the thermal contribution to the vibrational energy at T = 298 K, PV is the pressure versus volume, T = 298 K, and S vib is the vibrational entropy contribution.

Bulk Analysis
As already mentioned, the performance of the β-Li 3 PS 4 phase as an electrolyte is probably linked to a residual configurational disorder on the atomic scale. In the Pnma structure, the lithium atoms can occupy three distinct positions, corresponding to Wyckoff 8d, 4b, and 4c symmetry, for a total of 16 sites in which to place the 12 Li of the reference cell. At the macroscopic level, this translates into fractional occupation numbers, which experimentally, at finite temperature, proved to be 1.0 for 8d, 0.7 for 4b, and 0.3 for 4c. Since from the computational point of view it is necessary to maintain the symmetry of the space group, we have adopted the most probable configuration, occupying positions d and b, in line with most of the other theoretical studies [24,25].
First structure optimization of the Pnma crystal led to lattice parameters a = 12.96 Å, b = 8.08 Å, and c = 6.25 Å, in good agreement with experimental data [2], with a deviation of 1.09%, −1.07%, and 2.12%, respectively. However, in order to get rid of two imaginary frequencies, both related to the lithium atom in 4b, we released the symmetry constraints and performed a scan of geometries along the normal modes corresponding to the imaginary values [13]. The structure underwent structural deformation and turned into a Pn2 1 a crystal, a subgroup of Pnma, with lattice parameters a = 12.91 Å, b = 8.14 Å, and c = 6.23 Å, corresponding to a deviation of 0.70%, −0.97 %, and 1.80% from the experimental values of Homma et al. [2]. Its electron energy per cell was 0.09 eV lower than that of the Pnma system and the computed frequencies were then all positive. According to our results, in the Pnma, Li atoms have coordinates similar to those measured by Homma et al. [2] and computed by Yang and coworkers [24], while in the Pn2 1 a model they occupy 4a Wickoff positions, similar to those in the distorted B3C1 structure recently calculated by Lim et al. [26], see Table 1. It is worth mentioning that the positions of P and S atoms are practically the same in the two models, and there are also other more than reasonable similarities given the substantial equivalence between the two structures.
From an analysis of the chemical bonding profile along the Li-S and P-S directions, using the QTAIMC indicators, we see that the bonding framework is practically the same in the Pnma and Pn2 1 a models, see Table 2. The four almost equivalent P-S bonds can be classified as covalent interactions based on the topological values in r = r BCP : i.e., a high charge density, the negative values of the Laplacian, and a V/G ratio greater than 3. The four Li-S bonds have a so-called transitory nature, sharing features with both the ionic (positive values of the Laplacian in r BCP ) and covalent interactions (cylindrical shape of the charge density around the bond axes and directionality). The small differences found in the topological indicators are due to slight discrepancies in the first coordination spheres of Li atoms at different Wyckoff positions. Based on this analysis, it can be concluded that the Li-S interactions are weak and can be easily broken, leaving the lithium ion free to migrate within the porous material. This obviously plays an important part in Li diffusion. Table 2. Value of several local quantities at the bond critical points in Pn2 1 a β-Li 3 PS 4 structure: the bond distance (Å), the electron charge density (ρ(r)) (Å −3 ), the Laplacian of the density (∇ 2 ρ) (Å −5 ), the ratio between the potential energy density and the kinetic density (|V|/G), the bond degree (H/ρ) (a.u.) (i.e., the ratio between the total energy density and electron density), and the ellipticity ε. The indexes are related to the atom label. The computed mechanical properties for both the Pnma and Pn2 1 a are close to literature data [24], see Table S2. In particular, the β-structure is considered a soft material and this is a property of great interest in the case of SSE devices as the volume of the various components varies during battery operation, and any rigidity would make the structure unstable. A close comparison shows that the Pn2 1 a presents Young and Bulk moduli greater than those of the Pnma, highlighting a slightly higher stiffness, but it remains an extremely deformable system.
Since the differences between structures are minimal and boil down to thermodynamic stability, in the following we will refer only to the Pn2 1 a, being a global minimum in the potential energy surface.
The Li 3 PS 4 Pn2 1 a model is an insulator with an indirect band gap of 4.75 eV between the Γ and Z points. This value, computed at the PBE0 level, is in fairly good agreement with the one measured by Rangasamy et al. [27] of 5 eV and, as expected, is slightly overestimated when compared with those calculated in a previous theoretical study at the GGA-PBE and HSE06 [24] level. The band structure, reported in Figure 1, is very similar to those proposed by Lim and co-workers for the B3C1 and B4C0 systems [26]. From a qualitative point of view, the presence of wide, well-dispersed bands indicate small effective masses, band velocities, and large electron mobility. The projected density of states (DOS) shows that the low-energy electronic transitions occur between states that belong mainly to the 2p orbitals of the sulfur atoms while Li and P provide an almost negligible contribution to the bands around the Fermi level. Nanomaterials 2022, 12, x FOR PEER REVIEW 6 of 13

β-Li3PS4 Surfaces
Although unstable surfaces are more attractive for applications such as catalysis and photocatalysis, the materials involved in lithium-ion batteries must be stable due to the many changes and stresses they undergo during operation. In fact, there are different crystallographic planes along which the bulk phase can be reasonably cut, and each surface thus exposed can give rise to different terminations. For this reason, we have developed a strategy to identify the most stable surfaces, and only subsequently do we proceed to their characterization. Finally, to complete the characterization of the bulk material, we computed its Raman spectra and compared the results with experimental literature data [5,28]. The good agreement is encouraging and confirms the reliability of the adopted model for the bulk structure. The intense peak at 427 cm −1 , Figure 2, represents the characteristic signal of the β-structure, measured at 422 cm −1 due to the collective vibrations of the (PS 4 ) units. The wide band between 150 and 300 cm −1 contains the frequencies of all Li-S modes while weak signals above 500 cm −1 are assigned to the P-S bonds.

β-Li3PS4 Surfaces
Although unstable surfaces are more attractive for applications such as catalysis and photocatalysis, the materials involved in lithium-ion batteries must be stable due to the many changes and stresses they undergo during operation. In fact, there are different crystallographic planes along which the bulk phase can be reasonably cut, and each surface thus exposed can give rise to different terminations. For this reason, we have developed a strategy to identify the most stable surfaces, and only subsequently do we proceed to their characterization.

β-Li 3 PS 4 Surfaces
Although unstable surfaces are more attractive for applications such as catalysis and photocatalysis, the materials involved in lithium-ion batteries must be stable due to the many changes and stresses they undergo during operation. In fact, there are different crystallographic planes along which the bulk phase can be reasonably cut, and each surface thus exposed can give rise to different terminations. For this reason, we have developed a strategy to identify the most stable surfaces, and only subsequently do we proceed to their characterization.
First, we noticed that the tetrahedron (PS 4 ) represents an essential structural unit to maintain surface stability, so in modeling the various slabs we tried to preserve the integrity of these units as much as possible. Then, for each crystallographic plane, we designed slabs of increasing thickness, starting from the minimum reference cell, which contains-as in the bulk-four units of Li 3 PS 4 . We explored planes corresponding to the following Miller indices: (100), (010), (001), (110), (011), (111), (210), (211), according to the distribution density of the diffraction peaks as reported in Ref. [29], and, for each surface, different termination patterns, as LiS 3 , PS 3 , PS 2 , SLi 2 , SPLi 2 , LiS 2 , SPLi, and LiS, were modeled. All the results are collected and described in the Supplementary Information.
The complete optimization of the four unit slabs, i.e., atomic positions and lattice parameters, has led to a significant reduction in candidates because some surface models had undergone a drastic rearrangement and become particularly unstable. Moving on to the eight-unit structures (64 atoms each), we had already found the most stable surface terminations for each of the surfaces, so only these have been re-optimized. The structural and energetic information of the final models, containing eight units of Li 3 PS 4 , are reported in Table 3, while the initial and fully relaxed structures of the four most stable slabs, namely the (100), (010), (011), and (210), are shown in Figure 3. Table 3. Surfaces' termination (up/down), surface energy (E surf , in J/m 2 ), formation energy (E form , in kJ/mol), and band gap (E gap , in eV) for the nine surfaces of β-Li 3 PS 4 with eight units each. rameters, has led to a significant reduction in candidates because some surface models had undergone a drastic rearrangement and become particularly unstable. Moving on to the eight-unit structures (64 atoms each), we had already found the most stable surface terminations for each of the surfaces, so only these have been re-optimized. The structural and energetic information of the final models, containing eight units of Li3PS4, are reported in Table 3, while the initial and fully relaxed structures of the four most stable slabs, namely the (100), (010), (011), and (210), are shown in Figure 3.   In some cases, as for (001) and (110), the surface geometry after the optimization undergoes a significant reconstruction. In contrast, the LiS 2 -(100) surface, which is the most stable, experiences minimal distortion.
One interesting aspect to point out is that the four more stable slabs have very different electronic structures and properties. Their surface states, hereinafter referred to as band gap for uniformity of notation, vary from 2.7 eV for (010) to 4.6 in the case of (100). The DOS profiles of (100), (011), and (210), reported in Figure 4, are very similar: the valence bands are characterized by a strong contribution of the S-2p orbitals, only slightly hybridized with Li-s and P-2p functions while the conduction bands are mainly due to the P-2p orbitals. The (010) slab shows a rather peculiar density of states with a significant S character of the lower virtual bands and a consequent decrease in the band gap.  The Hirshfeld charges [30] were evaluated and compared with those of the bulk in Table 4 (and Table S4 for the less stable two-dimensional models). As expected, the charges of the atoms in the internal layers-which coordination is preserved and did not undergo structural rearrangements during the optimization-are very close to those in The Hirshfeld charges [30] were evaluated and compared with those of the bulk in Table 4 (and Table S4 for the less stable two-dimensional models). As expected, the charges of the atoms in the internal layers-which coordination is preserved and did not undergo structural rearrangements during the optimization-are very close to those in the bulk. In the case of the (100) slab, the surface atoms also show almost the same charge as those in the inner layers. On the contrary, the lithium and surface sulfur in (010) and (011) undergo a significant depletion of their atomic charge, as a consequence of the strong structural deformation. Table 4. Hirshfield charges of β-Li 3 PS 4 Pn2 1 a bulk and eight-unit surfaces for the surface (indicated by *) and internal atoms. The electrostatic potential maps along the z-direction are drawn in Figure 5 for the (100) and (010) surfaces. In the (100) slab, the potential is symmetric with respect to the central xy plane and has the same shape on the two surfaces. The net electric dipole along the non-periodic z-direction is zero and the film can grow, preserving stability. The situation is different in the case of slab (010), which exposes two different surfaces, one rich in lithium and the other dominated by sulfur atoms, see panel (b). A charge gradient along z is then established which increases with increasing thickness. On the one hand, the dipole opposes the growth of the film, but on the other hand, it can favor the migration of lithium ions along the structure.  The electrostatic potential maps along the z-direction are drawn in Figure 5 for the (100) and (010) surfaces. In the (100) slab, the potential is symmetric with respect to the central xy plane and has the same shape on the two surfaces. The net electric dipole along the non-periodic z-direction is zero and the film can grow, preserving stability. The situation is different in the case of slab (010), which exposes two different surfaces, one rich in lithium and the other dominated by sulfur atoms, see panel (b). A charge gradient along z is then established which increases with increasing thickness. On the one hand, the dipole opposes the growth of the film, but on the other hand, it can favor the migration of lithium ions along the structure. To complete the analysis, the three-dimensional charge density of the two slabs was calculated and superimposed on the electrostatic potential to visualize the regions where the positive/negative potential is maximally concentrated, see Figure S10.
Finally, we used Wulff's theory to visualize the ideal nanocrystal, which is the one that would form based on the calculated surface energies. Using the four most stable structures, we obtained the shape shown in Figure 6a. Then, it was possible to change the relative stability between different 2D models by gradually increasing a certain amount To complete the analysis, the three-dimensional charge density of the two slabs was calculated and superimposed on the electrostatic potential to visualize the regions where the positive/negative potential is maximally concentrated, see Figure S10.
Finally, we used Wulff's theory to visualize the ideal nanocrystal, which is the one that would form based on the calculated surface energies. Using the four most stable structures, we obtained the shape shown in Figure 6a. Then, it was possible to change the relative stability between different 2D models by gradually increasing a certain amount (~0.20 eV) of the surface energy of some surfaces with respect to others. In this way, the growth of different surfaces can be favored or inhibited and the formation of different crystalline morphologies can be envisaged. In particular, in the ideal nanocrystal the (100), (011), and (210) surfaces dominate. By increasing the stability of one surface, we obtained the nanostructures as in Figure 6b. The structures in panel (c) emerge from the stabilization of two surfaces at the same time. Finally, in panel (d), the structure is represented when both the (100) and (011) surfaces are preferred: as can be seen from the inset, this morphology closely resembles that of the ultrathin nanoplates synthesized by Hood and co-authors [31]. These results show that reliable predictions are possible based on our calculations, i.e., the (100) and (011) films are among the most stable. Additionally, the different shapes for β-Li3PS4 nanocrystals are not only plausible but desirable due to the new properties they may have.

Conclusions
A DFT study of β-Li3PS4 bulk and surface stability was conducted using the CRYS-TAL program, and by applying all-electron basis-sets and the PBE0 functional. This study mainly aims to investigate the stability of the 2D thin film of β-Li3PS4 for their use in lithium batteries as a solid electrolyte.
A full and symmetry-unconstrained geometry optimization of the experimental refined Pnma bulk has led to an equilibrium structure, belonging to the Pn21a subgroup, in which the P-S lattice remains unchanged while the 12 Li atoms occupy three different Wyckoff positions: 4a′, 4a″, and 4a‴. Then, a set of slabs, whose thickness varies between 9 to 23 Å , were modeled by cutting the bulk along different crystallographic directions. Based on the calculated surface formation energies and Gibbs free energies, the (100) surface resulted as the most stable, followed by (210), (011), and (010). In addition, the (100) and (210) films showed interesting features such as good mechanical and thermodynamic stability, a high concentration and good mobility of Li ions (i.e., the Li atoms are loosely In this way, the growth of different surfaces can be favored or inhibited and the formation of different crystalline morphologies can be envisaged. In particular, in the ideal nanocrystal the (100), (011), and (210) surfaces dominate. By increasing the stability of one surface, we obtained the nanostructures as in Figure 6b. The structures in panel (c) emerge from the stabilization of two surfaces at the same time. Finally, in panel (d), the structure is represented when both the (100) and (011) surfaces are preferred: as can be seen from the inset, this morphology closely resembles that of the ultrathin nanoplates synthesized by Hood and co-authors [31]. These results show that reliable predictions are possible based on our calculations, i.e., the (100) and (011) films are among the most stable. Additionally, the different shapes for β-Li 3 PS 4 nanocrystals are not only plausible but desirable due to the new properties they may have.

Conclusions
A DFT study of β-Li 3 PS 4 bulk and surface stability was conducted using the CRYSTAL program, and by applying all-electron basis-sets and the PBE0 functional. This study mainly aims to investigate the stability of the 2D thin film of β-Li 3 PS 4 for their use in lithium batteries as a solid electrolyte.
A full and symmetry-unconstrained geometry optimization of the experimental refined Pnma bulk has led to an equilibrium structure, belonging to the Pn2 1 a subgroup, in which the P-S lattice remains unchanged while the 12 Li atoms occupy three different Wyckoff positions: 4a , 4a", and 4a . Then, a set of slabs, whose thickness varies between 9 to 23 Å, were modeled by cutting the bulk along different crystallographic directions. Based on the calculated surface formation energies and Gibbs free energies, the (100) surface resulted as the most stable, followed by (210), (011), and (010). In addition, the (100) and (210) films showed interesting features such as good mechanical and thermodynamic stability, a high concentration and good mobility of Li ions (i.e., the Li atoms are loosely bonded and the lattice is extremely porous), and a potentially small lattice mismatch with materials such as Li 2 S, widely used as a passivator in lithium technology. The (010) structure presents original characteristics, such as peculiar surface states and a net dipole moment along the non-periodic z-axis, which could eventually increase the diffusion rate of the lithium ions. Finally, Wulff's theory was applied to derive the shapes of possible nanocrystals. As proof, it has been verified that the extra-stabilization of both the (100) and (011) surfaces, done by hand, leads to nanoplates very similar to those synthesized experimentally.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/nano12162795/s1, Table S1: Cell parameters and band gap of β-Li 3 PS 4 ; Table S2: Mechanical properties of β-Li 3 PS 4 ; Table S3: Possible termination along +z of each analyzed surface. The "X" represents the termination found in the surface and the '-' the termination not found; Table S4: Hirshfield charges of β-Li 3 PS 4 bulk and eight-unit surfaces for the surface (indicated by *) and internal atoms; Figure