High temperature equilibrium of 3D and 2D chalcogenide perovskites

Chalcogenide perovskites have been recently under the researchers spotlight as novel absorber materials for photovoltaic applications. BaZrS$_3$, the most investigated compound of this family, shows a high absorption coefficient, a bandgap of around 1.8 eV, and excellent environmental and thermal stability. In addition to the 3D perovskite BaZrS$_3$, the Ba-Zr-S compositional space contains various 2-D Ruddlesden-Popper phases Ba$_{x+1}$Zr$_x$S$_{3x+1}$ (with $x=$ 1, 2, 3) which have recently been reported. In this work it will be shown that at high temperature the Gibbs free energies of 3D and 2D perovskites are very close, suggesting that 2D phases can be easily formed at high temperatures. The analysis of the product of the BaS and ZrS$_2$ solid-state reaction, in different stoichiometric conditions, present a mixture of BaZrS$_3$ and Ba$_4$Zr$_3$S$_{10}$. To carefully resolve the composition, XRD, SEM and EDS analysis were complemented with Raman spectroscopy. For this purpose, the phonon modes, and the consequent Raman spectra, were calculated for the 3D and 2D chalcogenide perovskites, as well as for the binary precursors. This thorough characterization demonstrates the thermodynamic limitations and experimental difficulties in forming phase-pure chalcogenide perovskites through solid state synthesis, and the importance of using multiple techniques to soundly resolve the composition of these chalcogenide materials.


Introduction
Chalcogenide perovskites (CPVK) have recently been proposed as possible non-toxic alternatives to lead-based perovskites for photovoltaic applications thanks to their promising optoelectronic characteristics, including a bandgap that is suited for tandem solar cell applications. [1][2][3][4] Chalcogenide perovskites follow the formula ABX3 with A, B and X representing an alkaline earth cation (2+), a transition metal cation (4+) and a chalcogenide anion (2-) respectively. [5][6][7] The most studied compound in this family is BaZrS 3 , but other compositions have been reported. [8,9] Accordingly, chalcogenide perovskites have been under recent scrutiny to evaluate and validate a variety of optoelectronic parameters. Reports suggest the possibility of tunable bandgaps through compositional engineering with high absorption coefficients in the visible range, [10] low effective masses [11] and, resultantly, high carrier mobilities. [12] Moreover, these materials should present improved thermal and chemical stability compared with the hybrid halide perovskites, being resistant towards high temperatures and atmospheric conditions. [13] Notwithstanding these promising features, the development of CPVK-based devices is still hindered by the synthetic procedure necessary to prepare these perovskites. In practice, these materials need very high temperatures to be crystallized in the desired phase. For example, BaZrS3 is produced by the solid-state reaction of the elemental or binary precursors at 800-1100°C for several hours or days. This requirement for high temperatures is limiting as it does not easily allow for thin film processing and device integration, and can lead to chalcogen loss via volatile precursors.
Understanding reaction mechanisms can help design fabrication procedures potentially having lower temperatures, and recent research has started to explore which processes may allow for this.
For example, it has been proposed that favoring the formation of BaS 3 , which has a lower melting point than BaS, can trigger the formation of BaZrS3 at temperatures as low as 600°C. [14] Similarly, it has been reported that an excess of sulphur in the reaction mixture can significantly reduce the reaction time and temperature. [15] The tuning of the precursors' stoichiometry can provide a path to find milder synthetic conditions for BaZrS3 and other chalcogenide perovskites, but it can also lead to the formation of other unwanted perovskite phases. The ideal perovskite structure is formed from a three-dimensional network of corner sharing octahedra. However, there are other closely related perovskite-like phases that can be formed. For example, low dimensional Ruddlesden-Popper (RP) phases are known for chalcogenide perovskites, similarly to the oxide and halide perovskites. Ban+1ZrnS3n+1 for n = 1, 2, 3, (with the presence of low-and high-temperature polymorphs for Ba3Zr2S7) have been reported in addition to the 3D structure BaZrS3 (in which n = ∞). [6,13,[16][17][18] For the corresponding perovskites with Hf (Ban+1HfnS3n+1), additional Ruddlesden Popper phases with n = 4, 5 have been reported, [19] while no evidence is present for n > 3 in the Zr series. In these low dimensional phases, layers of 3D CPVK are alternated with layers of BaS, as represented in Figure   1. Interestingly, the bandgaps of these chalcogenide RP perovskites decrease as n increases, contrary to the oxide and halide counterparts, where the bandgap widens as n increases. [20][21][22] The different crystalline structure not only affects the bandgap of the perovskite, but can also affect other functional features, such as carrier transport and thermal or chemical stability. A deep understanding of the synthetic reaction, especially when non-stoichiometric conditions are explored, will be essential to control the formation of competing phases with distinct properties.
The most widely used technique for the assessment of different phases is X-ray diffraction (XRD).
However, due to the structural similarity of the 3D and the 2D perovskites, the diffractogram peaks overlap, which make differentiating between species challenging. For example, in the case of samples formed by the mixture of BaZrS3 and Ba3Zr2S7, long data collection times and rigorous refinement (supported by compositional techniques such as Energy Dispersive X-ray Spectroscopy, EDS) are necessary to give a quantitative estimation of their compositions. As such, even with diffractograms presenting good angular resolution, intensity, and angular range, ambiguous assignation can still occur. It follows that good practice, especially for this family of materials, would be to combine crystallographic analysis with other complementary techniques that may also provide a clearer distinction between species. Vibrational spectroscopy is a good candidate technique to carry out this role, as it probes the local structure of a material, in contrast to XRD which probes the bulk response. In this letter it will be shown that the 3D and 2D perovskite structures each have a unique vibrational fingerprint that better distinguishes between materials in the Ba-Zr-S system and that Raman spectroscopy is therefore well suited for checking the phase purity of the compounds.
In this study XRD and Raman spectroscopy are combined to assess the main composition of the product resulting from the solid-state reaction of BaS and ZrS2 at various ratios. To aid in our analysis, a first-principles thermodynamic model is used to demonstrate that at high temperatures the Gibbs free energy of the 3D and 2D CPVK materials are only a few kJ/mol apart, indicating that both can be formed during high temperature synthesis (>1000K). In addition, the Raman spectra for all known binary and ternary compositions in the Ba-Zr-S phase systems are calculated from first-principles. Despite the use of excess ZrS2 in the reaction mixtures, all the resulting powders show deficiency of Zr and S, and present mixtures of BaZrS3 and the Ruddlesden-Popper phase Ba4Zr3S10. It will also be shown that peaks' assignation can be done more confidently when Raman and XRD are used simultaneously. All calculated Raman spectra and the thermodynamic analysis code are published in open-access repositories alongside this work, allowing the adoption of our approach to other studies of BaZrS3 synthesis.

Methods:
Experimental procedure: Synthesis of powders: BaS (99.7%. Alfa Aesar) and ZrS2 (99.99%, Alfa Chemistry) were used without further purification. 300 mg of powder formed by 1:1, 1:1.05, 1:1.1 and 1:1.2 molar mixtures of BaS and ZrS2 (called Zr_0, Zr_5, Zr_10, Zr_20 respectively) were finely ground with an agate pestle and mortar, loaded into carbon-coated quartz ampules, purged with argon 3 times and sealed under a vacuum of 10 -5 mbar. Each mixture was then placed in a single zone furnace at 500°C, and the temperature was increased to 900°C at a rate of 200°C/h. The powders were kept at 900°C for 5 days, and finally quenched in water.

Computational details:
Competing phases of BaZrS3 were identified using the Materials Project database, [23] with all Ba-Zr-S compounds within 0.5 eV above the convex hull considered. This energy range has been shown to cover the 90 th percentile of all metastable materials reported within Materials Project. [24,25] First principles density functional theory (DFT) calculations were carried out with the all-electron numeric atom-centered code FHI-aims. [26] Minimum-energy crystal structures were found using parametrically constrained geometry relaxation. [27] The self-consistent field criteria was set to 10 -7 e/Å 3 and 10 -6 eV/Å for electron density and force respectively. The structures were relaxed until the maximum force component was below 5x10 -3 eV/Å. All other inputs were set to the default value within FHI-aims. All relaxations and phonon calculations were performed with the PBEsol [28] functional with a tight basis set. Electronic band structures and total energies were calculated using the HSE06 [29] functional alongside inclusion of spin-orbit coupling.
The resulting formation energies and band gaps are reported in Table S1. Phonon band structures were evaluated using the finite difference method with a 0.01 Å step size, as implemented in Phonopy. [30] The supercell size and k-point spacing used for phonon calculations are presented in Table S2. Raman intensities and peak positions were generated through Phonopy-Spectroscopy with a two-point finite difference along displacements, and are reported in the Supporting Information. [31] A Lorentzian peak width of 1 cm -1 was set on the peak positions obtained from the analysis. Macroscopic dielectric tensors were evaluated using the real space density functional perturbation theory method implemented in FHI-aims [32] with the PBE functional. [33] Gibbs free energies were calculated using Phonopy and the ThermoPot package. [34] Imaginary modes were omitted from our calculation of the thermodynamic partition function and so do not contribute to the Gibbs free energy. An online repository containing i) analysis code for generating Figure 1a and Figure S2; ii) the raw data from electronic structure calculations; and iii) the data used to generate Raman spectra is available at https://github.com/NU-CEM/2022_BaZrS3_High-T_equilibrium.

Results and discussion
During BaZrS3 synthesis competing phases can form from both the constituent elements (Ba, Zr, S) and/or from external impurities (e.g., O2). In order to focus our analysis efforts on the thermodynamically accessible competing phases , ab-initio thermodynamic calculations were performed. As in this study synthesis was carried out within a closed ampule under vacuum, phases in the Ba-Zr-S system only were considered. Competing phases were initially identified using the Materials Project database, [35] followed by re-calculation of the total energies using a higher level of theory (HSE06 functional with spin-orbit coupling) to reproduce lattice constants and bandgaps in agreement with published results from experimental studies: the formation energy, lattice constants, electronic band gap and a comparison to experimental data for all systems considered are provided in the Table S1. Five competing ternary phases were identified: Ba4Zr3S10 (I4/mmm), Ba3Zr2S7 (I4/mmm), Ba3Zr2S7 (P42/mnm), Ba3Zr2S7 (Cmmm) and Ba2ZrS4 (I4/mmm). All ternary compounds are in the Ruddlesden Popper perovskite analogue series, Ban+1ZrnS3n+1.
It follows that there are three ternary-to-ternary decomposition mechanisms to consider: To predict the relative phase stabilities a previously published methodology was followed to calculate the change in Gibbs free energy (ΔG) of each process, [36,37] as this is the potential that is minimised in equilibrium. It is important to emphasize that this calculation does not consider the effects of lattice expansion or anharmonic vibrations. However, despite these limitations, this methodology has been used to successfully predict the temperature-pressure stability window for Cu2ZnSnS4. [36] The G for all processes, reported in Figure 1a, is endothermic (positive valued) across the full temperature range considered (100K -1300K  [6] Crystal, electronic and vibrational structure information for BaZrS3 and for the lowest energy competing phase, Ba4Zr3S10, are shown in Figure 1b-g. The electronic band structure is calculated with the HSE06 functional and includes spin-orbit effects, leading to an accurate predicted value of 1.72eV [38] for the 3D perovskite direct band gap and a predicted value of 1.13 eV for the indirect band gap of the RP phase. This is the smallest predicted band gap for RP Ban+1ZrnS3n+1 materials considered, and is in line with previous reports of a decreasing band gap with increasing n. [20][21][22] Whilst an indirect bandgap can lead to a decreased absorption coefficient near the band edge, we expect this to be counter-balanced by the flatter band dispersion for Ba4Zr3S10, resulting in a larger density of states.
At room temperature BaZrS3 is reported to form in in space group Pnma, [9] which is a distortion of the idealised cubic perovskite. RP phases are reported to form in the space group I4/mmm. [20][21] It is found that BaZrS3 in the Pnma phase is dynamically stable with positive phonon modes across the Brillouin zone. In contrast, Ba4Zr3S10 is dynamically unstable with imaginary phonon modes at the zone boundaries. This indicates the presence of a symmetry lowering transition to a more stable phase at low-temperature, which is a common feature of halide and oxide perovskite materials. [39,40] Ba4Zr3S10 in the space group Fmmm has also been previously reported in the literature. [18] However this corresponds to a small lattice expansion and increase in the c/a lattice parameter ratio, so is unlikely to result from distortions along zone boundary phonon modes of the I4/mmm phase. In Figures S1 we confirm that the high temperatures structures for Ba4Zr3S10 (Fmmm), Ba3Zr2S7 (I4/mmm) and Ba2ZrS4 (I4/mmm) also produce imaginary phonon modes in the harmonic approximation. Observing the diffractograms it can be noted that the samples present high crystallinity and, except the Zr_20, present complete conversion to the ternary phase, with no evidence of unreacted BaS ( Figure S3a). The presence of unreacted ZrS2 is more challenging to assess, as its characteristics peak is located at 32.2° ( Figure S3a), where the samples present numerous small peaks which could derive also from the 2D or 3D perovskite. As discussed later in the text, Raman analysis simplified the assignation, confirming the presence of ZrS2 in the Zr_0 sample. The Zr_20 sample presents an additional peak at low angles which origin has not been identified. Given the poor quality of the XRD pattern of this sample, both in terms of crystallinity and phase purity, it has been excluded from further characterization. It has to be noted that the presence of other Ba-S and Zr-S binary compositions have been excluded observing their XRD patterns and Raman spectra ( Figure S3b). Similarly, the presence of oxides, sulfates and carbonates (all possible unwanted products in presence of air) have been excluded by comparing reported experimental XRD patterns and Raman spectra (Table S3).
The Zr_0 sample, which was expected to provide stoichiometric conversion to BaZrS 3 , [41] shows the presence of multiple phases that were assigned to the 3D BaZrS3 and to the RP phase Ba4Zr3S10.
It is important to highlight the similarity of the angular position and intensity of the BaZrS3, Ba3Zr2S7 and Ba4Zr3S10 XRD patterns, which complicates the experimental diffractograms resolution ( Figure S4). The Rietveld refinement presented in Figure S5 and Table S4-6 shows that the RP phase is the majority one, representing almost the 70% by mass. As the molar amount of ZrS2 in the starting mixture is increased by 5%, the diffractogram still reveals a mixture of phases, but with an overall predominance of the BaZrS3 over the RP phase.

b) Calculated (bottom) and measured (top) Raman spectra collected with 785 nm excitation wavelength
This can be concluded by the change in shape of the peak at 25.1° and by observing the high angle peaks ( Figure S4), and is confirmed by the Rietveld refinement. For the sample with a nominal 10% molar excess of ZrS2, Ba4Zr3Z10 returns to be the majority phase, surprisingly reaching 86% of the sample weight. In Figure 2b and Table 2  (83 and 133 cm -1 ) and the peak at 380 cm -1 are mainly assigned to the BaZrS3 phase (A 2 g, A 4 g and B 5 1g respectively), in agreement with previous works (Table 2). [5,40,41] It is worth stressing that the expected experimental spectra measured at room temperature should be red-shifted to lower frequencies compared to computational predictions, due to thermal expansion of the lattice which leads to a softening of the phonon modes. For the A 3 g and A 4 g modes in BaZrS3 the frequency shift at room temperature is reported to be approximately 8 cm -1 , [42] in agreement with what observed (a shift of 10.3 cm -1 and 6.6 cm -1 respectively). In contrast, the experimental spectra are blueshifted for the A 6 g and B 6 2g modes, as has been reported previously. [42,44] Importantly, these modes relate to displacements of the sulphur species only. In the sample Zr_0 a shoulder can be observed around 308 cm -1 , which is assigned to both the A 7 g mode of the BaZrS3 as well as the A 1 1g of the ZrS2. Since previously reported Raman spectra of BaZrS3 show that the A 7 g mode is visible only at low temperature, [42] the 308 cm -1 peak in the Zr_0 spectrum is likely to derive from unreacted ZrS2. This demonstrates that the presence of ZrS2 is easier to observe through Raman spectroscopy than through XRD.
Importantly, all the spectra present a peak at around 339 cm -1 which does not derive from the 3D perovskite. Comparing the measured spectrum with the calculated one for Ba4Zr3S10 it is possible to assign it to the secondary RP phase, which is also corroborated by XRD. The Zr_0 sample presents an additional peak at 416 cm -1 , which can be assigned to the Ag mode of the Ba4Zr3S10. This phonon mode involves displacement of the sulphur species only, as in the A 6 g and B 6 2g modes of the BaZrS3 phase, which can be used to rationalise the unexpected blue shift to higher frequencies for the experimental spectra.
Importantly, Raman spectra cannot be used for quantitative analysis unless calibration with phase pure materials is performed. The higher intensity of the BaZrS3 Raman peaks compared to the RP ones does not indicate a higher content in the sample, as demonstrated by the Rietveld quantifications. However, it is interesting to note how the ratio between the areas of the 339 cm -1 and 214 cm -1 peak decreases as 5% molar excess of ZrS2 is used in the starting mixture, suggesting that in this sample the orthorhombic 3dimensional phase is favoured over the 2D. This hypothesis is supported by the disappearance of the peak at 99 cm -1 (which does not derive from BaZrS3 and is attributed to the Ba4Zr3S10 mode located at 92 cm -1 ) and to the increased intensity of the 380 cm -1 peak. In the Zr_10 sample, instead, the ratio between the two peaks returns in favour to the Ba4Zr3S10, with a reduction of the BaZrS3 peaks, as confirmed by Rietveld analysis. It should be noted that the measured Raman spectra were compared against the spectra of the binary precursors ( Figure S3) as well as against the oxide counterparts (Table S3), [45][46][47][48][49] confirming their absence in the synthesised mixture.
So far it has been shown that, in these synthetic conditions, the reaction between binary precursors leads to the formation of a secondary 2D ternary phase. Zr-or S-poor conditions can trigger the formation of RP phases, and a careful compositional analysis is necessary to obtain sensible conclusions. For this reason EDS compositional analysis has been performed on all the synthesized powders and on the Zr_0 precursor mixture before heat treatment, and the results are reported in Table 3. The Zr_0 precursor mixture confirms that the Ba/Zr atomic ratio before heat treatment was 1:1.
Surprisingly, all the resulting powders after treatment present Zr deficient compositions, with the Ba/Zr atomic ratio decreasing in the Zr_0, Zr_5 and Zr_10 series. However it is worth stressing that the confidence range of these values is enlarged by the limitations associated with powders analysis with EDS (see Supporting Information). The starting material and the synthesised powders also present sub-stoichiometric amounts of sulphur, which become more unbalanced as the starting excess of ZrS2 is higher. To exclude accidental losses during the loading of the samples in the quartz ampules, a second batch repeating the Zr_0 conditions was prepared and characterized showing very similar XRD, Raman spectra and EDS compositions to Zr_0 ( Figure S6). The systematic loss of Zr in these synthesis hints to possible unwanted reactions with the quartz ampule, even if carbon-coated. More investigation is needed to address this phenomenon, which is out of the scope of this work. However this evidence hints that other precursors, rather than ZrS2, should be used, especially in high-temperature solid-state synthesis.
The observation of sub-stoichiometric amounts of Zr and S gives an additional explanation for the formation of Ba4Zr3S10. As the RP Ban+1ZrnS3n+1 family of materials are a ZrS2-deficient analogue of BaZrS3, the RP phases are expected to be more readily formed in Zr-or S-poor environments.
To explore this further, a ternary phase diagram for the Ba-Zr-S system was constructed using the first-principles thermodynamic model introduced earlier ( Figure S2). This allows a prediction of which products are formed for the sub-stoichiometric amounts of Zr and S as measured by EDS.
At high temperature all RP phases and most binary phases (all except BaS3) lie on the convex hull, so that the predicted products are very sensitive to composition. Using Ba-Zr-S composition values within the range of values reported for each sample in Table 3, the model predicts that Ba4Zr3S10 and BaZrS3 will be formed at at 900°C, alongside a smaller proportion of ZrS2 (Table S7). This suggests that there is an additional driving force for phase separation into BaZrS3 and Ba4Zr3S10 resulting from the under stoichiometric amount of Zr and S. Our modelling also suggests that the S concentration in the samples limits the formation of BaZrS3. This is why the Zr_10 sample, with the lowest proportion of sulphur measured by EDS, produces the highest proportion of Ba4Zr3S10.
On the other hand, this sample has the Ba/Zr ratio closer to unity, suggesting that the cation ratio is not the driving force for the preferential formation of 3D over 2D perovskites.
The compositional mapping presented in Figure 3 offers an additional insight in the nature of the synthesised powders. In the image, brighter zones correspond to areas of high atomic concentration, although it is recognised that shadowing effects give large dark areas (observable in the dark large areas in the mappings but not in the secondary electron images). Notwithstanding this, it may be observed that in all the samples Zr is less uniformly distributed than Ba and S, as  hinders the uniformity of the reaction environment possibly creating local regions even more Zrdeficient than the average composition, favouring the formation of 2D perovskites. Furthermore, the morphology of the synthesised crystals has been analysed through scanning electron microscopy (SEM) in both secondary electron and the more compositionally sensitive backscattered electron (BSE) modes. Looking at the BSE image in Figure S10 it can be noted that, in the Zr_0 sample, the crystals have a different appearance accordingly with the grain composition, as it could also be appreciated in Figure S7. Finally, this work has demonstrated the complexity of the Ba-Zr-S phase diagram, stressing the importance of using multi-experimental techniques to soundly resolve the reaction products of synthesis, and the necessity to find alternative synthetic routes involving lower temperatures and different precursors. Raman spectroscopy, here suggested as a technique complimentary to XRD, has been proven to provide a relatively easy and quick differentiator between CPVK phases, and even for several binary compositions of the Ba-Zr-S elements. In addition, the use of Raman spectroscopy in other CPVK synthesis environments (such as nanoparticle synthesis) may help to identify the reaction mechanisms where traditional techniques, such as XRD, cannot be used, or where the presence of organic materials complicates the analysis. Hopefully, the use of the database created in this work will support further progress in developing a low-temperature synthetic procedure for chalcogenide perovskites, promote increased control of the 2D-3D phase equilibrium during synthesis, and ultimately enable their thin film deposition and integration into optoelectronic devices.

Supporting information
High temperature equilibrium of 3D and 2D chalcogenide perovskites.

X-Ray Diffraction
The powder X-ray diffraction patterns were acquired using Rigaku SmartLab SE, in Bragg Brentano geometry and with Cu Kα1 (1.54056 Å) radiation. The diffractograms were collected from 7 to 80°, with long acquisition time (0.2°/min) to have sufficient angular resolution for Rietveld refinement.

Rietveld Analysis
Full-profile refinement and quantitative phase analysis was performed using Rietveld method as implemented in FullProf suite. [1,2]

Raman Spectroscopy
Raman spectra were collected using Horiba LabRAM Soleil system, equipped with a 785 nm laser. The power at the sample was measured to be 4.7 mW. The used objective was a 50X with a 11mm working distance, the confocal hole was 200 mm and the grating 1200 g/mm.

Scanning Electron Microscopy and Energy Dispersive X-Ray (SEM-EDS)
Scanning electron microscopy images were acquired using Tescan Mira 3 FEG-SEM, employing both secondary electron and back-scattered electron detectors. Energy dispersive Xray and elemental mapping were collected using a X-Max detector and Aztec software (Oxford Instruments), using 20KV electron acceleration and ZAF correction algorithm. The elemental compositions presented in Table 3 were obtained by averaging the values obtained over 10 different areas (600*600 µm) for each sample. These values are reported with a 5% uncertainty due to the nature of the measurement: EDS can be a very accurate technique to measure sample composition, provided that the sample has a flat and smooth surface, uniform thickness and density, and that a calibration standard for each element of interest is measured together with the samples. The powders analysed in this work presented very rough surfaces and strong variation in the sample thickness. The morphological characteristics of the samples can strongly affect the accuracy of the measurements, due to electron scattering and X-ray re-absorption. Additionally, analysis was done without standards for Ba, Zr, S (even if quantification were been verified on other standard samples). For this reason, even if precise measurements were carried out (ie, low standard deviation between the 10 points analysed for each sample, ranging around 0.5%), the overall accuracy is reduced by the standardless analysis, and for this reason it was assigned to 5%.

First-principles predictions of Raman peak position and intensity
Calculation details are included in the main text.   [5] 2.397 (direct) ZrS2 (P-3m1) -1.734 a = 3.691, c = 6.611, a = 3.663, c = 5.821 [6] 1.595 (indirect) ZrS3 (P21/m) -1.378 a = 3.6572, b = 5.2023, c = 9.504, a = 3.624, b = 5.124, c = 8.98 [7] 1.876 (indirect) ZrS (Fm-3m) -1.459 a = 5.24814, a = 5.240 [8] 0 ZrS (P4/nmm) -1.592 a = 3.61935, c = 5.54735, a = 3.55, c = 6.31 [9] 0 Zr3S4 (Fd-3m) -1.371 a = 10.37, a = 10.25 [8] Figure S4. Phonon dispersions showing imaginary modes for a) Ba4Zr3S10 in the Fmmm phase b) Ba3Zr2S7 in the Immm phase c) Ba2ZrS4 in the Immm phase. As is common in the literature, we omitted the imaginary modes from our calculation of the thermodynamic partition function and so they do not contribute to the Gibbs free energy. This approximation can be justified when considering that the thermodynamic properties are calculated statistically across all phonon branches, and that there is a relatively small number of imaginary modes across the whole Brillouin Zone. In addition, it has been shown for a number of related systems that renormalising the imaginary phonon modes to an effective real harmonic frequency, one strategy for inclusion within the Gibbs free energy term, has little impact on the calculated free energies. [15,16] Figure S2. Phase diagram for the Ba-Zr-S system at 1200K. The phase diagram has been calculated using ThermoPot [17] and Pymatgen. [18] The blue points indicate stable phases on the convex hull, whilst the red square indicates an unstable phase.