Structural Aspects of the Superionic Transition in AX2 Compounds With the Fluorite Structure

Some AX 2 binary compounds with the fluorite structure (space group Fm3¯m ) are well-known examples of materials exhibiting transitions to ionic superconducting phases at high temperatures below their melting points. Such superionic states have been described as either highly defective crystals or part-crystal, part-liquid states where the A ions retain their crystalline order whilst the X ions undergo partial melting. However, no detailed description of the structure of these phases exists. We present here the results of our investigation of the structural changes that occur during these transitions and the structural characteristics of the resulting superionic materials. This work is based on atomic-scale molecular dynamics modelling methods as well as computational diffraction techniques. We employed a set of empirical potentials representing several compounds with the fluorite structure to investigate any potential-dependent effect. We show the importance of small-scale structure changes, with some local environments showing a hexagonal symmetry similar to what is seen in the scrutinyite structure that has been documented for example in UO2.

Besides the increased ionic conductivity, a known feature of the superionic transition is a peak in the constant-volume heat capacity C P (T) with a characteristic Λ shape (Goff et al., 1991). The cause for this is an increase in the enthalpy H(T) beyond what would be expected from the trend measured from the same materials at low temperatures (Naylor, 1945;Bredig (1963, 1968)). This has been commonly referred to as anomalous heat (Dworkin and Bredig, 1968;Szwarc, 1969). Coincidently, an anomalous increase in the lattice parameter of the materials can also be measured, or equivalently a Λ peak in the linear thermal expansion coefficient α(T) (Goff et al., 1991). Both peaks have their highest point at the same temperature for any given compound, which is conventionally taken as the superionic transition temperature T S .
The fluorite structure (space group Fm3m) is common for many compounds with the AX 2 form (Hayes (1974)). In the classic fluorite structure, following the prototype CaF 2 , A is a cation-such as Ca 2+ , Sr 2+ , U 4+whereas X is an anion-commonly F − , Cl − , or O 2− . This structure can be viewed as the combination of a face-centred cubic (FCC) sublattice containing the A chemical species and a simple cubic (SC) sublattice with the X species, as shown in Figure 1. The anti-fluorite structure is very closely related, the only difference being that in this case A is an anion, and X is a cation. Both Na 2 O and Li 2 O are common examples of the antifluorite structure. However, this difference in ionic character does not change the geometry of the crystals. In the rest of the article, both structures will be referred to as either fluorite or Fm3m, but the conclusions apply equally to anti-fluorites.
Compounds with the fluorite structure have been very well studied for industrial applications (e.g. UO 2 , CeO 2 , yttria-stabilised zirconia), or because they are a simple model for ionic conductors (Hayes, 1974). Their tendency to behave as superionic conductors at high temperatures has been a particular area of interest since the 1970s. The nature of the hightemperature superionic phase in AX 2 compounds with the fluorite structure is not fully explained yet, and conflicting interpretations exist. One of them is that high conductivity is the result of the availability of many crystalline sites for diffusion resulting from extensive Frenkel disorder and defects clustering on the X sublattice (Hutchings et al., 1984;Hull et al., 2011).
Another view is that the superionic phase is the combination of a melted X sublattice and an A sublattice that retains its crystalline character, and thus the integrity of the material (Boyce et al., 1977;Castiglione and Madden, 2001). Recent works on the dynamical behaviour of the superionic phase of UO 2 and Li 2 O seem to validate this description, showing string-like collective diffusion mechanisms otherwise common in supercooled liquids (Gray-Weale and Madden, 2004;Annamareddy and Eapen, 2017;Zhang et al., 2019). This interpretation is also strengthened by measurements of entropy changes during superionic transitions, which produced results consistent with some partial melting (O'Keeffe, 1976;Boyce et al., 1977).
In this work, we set out to investigate the structural aspects of the superionic transition in compounds with the fluorite structure. We aim to provide some insight on key structure changes happening at that transition and on the structural characteristics of the superionic phase. Molecular dynamics (MD) techniques have been recognised as a powerful tool for a very long time for the study of ionic conducting oxides (Catlow et al., 1978). It is particularly well-suited for this type of investigations that bridge thermodynamical properties and structural features. Indeed, its atomic-scale resolution makes accessible structure details such as point defects, local environments, or structure changes. This has resulted in a large collection of empirical potential that can be used to simulate a variety of compounds. Instead of focusing on one specific material, we considered a list of them for which successful potentials have been published. This way, we intend to describe phenomena that are characteristic of ionic compounds in the Fm3m structure rather than specific to one composition or one potential.
Structures from MD simulations of superionic phases can be very difficult to analyse because of temperature-induced effects. This includes not only structural defects that form normally at high temperatures, such as Frenkel pairs, Schottky trios, or larger clusters, but also the resulting strain fields that distort the structure, and large thermal vibrations. All of this obscures the underlying structure of the materials and makes any structural analysis challenging. To overcome this, we employed computational diffraction techniques to generate virtual X-ray diffraction (XRD) patterns. These provide another way of characterising structural features that are difficult to describe otherwise, by showing structures in the reciprocal space that are not subject to the same perturbations as the atomic positions at high temperatures.

Empirical Potentials
We selected a broad range of empirical potentials that have been designed to simulate crystals with the fluorite structure. Thanks to the sustained interest for AX 2 materials with the Fm3m structure, a large number of such potentials are available and can be used for our purposes. The final set of compounds that can be simulated using the selected potentials included seven different chemical compositions, with some fluorides, oxides and one chloride. Most of these compositions were represented by several different potentials listed in Table 1 with their main characteristics.
The seven selected compounds include some variety in the properties of both A and X ions. The X ions were Li, F, Cl, or O, having electric charges of respectively +1, −1,−1, and−2. The ions on the A sublattice were O, Ca, Sr, Ba, Pb, or U with masses ranging from 16 u to 238 u, and with A X mass ratio from 2.31 (Li 2 O) to 14.88 (UO 2 ). The potentials we used take different analytical forms. Most of the potentials, particularly older ones, use simple Buckingham interactions. However, they are not identical to each other and differ in the numerical values of their parameters and their fitting procedure. Indeed, they have been initially fitted to reproduce different properties such as defect energies at room temperatures, which are not necessarily related to the superionic transition. In some more recent potentials, the electric charges of the different species are a parameter that was adjusted during the potentials fitting in the same way as the parameters of the non-electrostatic pair interactions. This reflects the imperfect ionic nature of these materials. Amongst the potentials used in this study, it is the case for the potentials representing Li 2 O and UO 2 . Some potentials use analytical forms other than Buckingham, such as Morse (Pedone potential for Li 2 O), 4-terms Buckingham (Morelon potential for UO 2 ), or additional embedded atom model (EAM) many-body contributions (CRG potential for UO 2 ). The Pedone potential was not even designed to simulate the Fm3m structure, but instead complex minerals and glasses. The analytical form of the electrostatic interactions was the Wolf summation (Wolf et al., 1999).
In some instances, several potentials were proposed in the original references. This is the case for the Catlow potential for BaF 2 , CaF 2 , and SrF 2 , and the Oda potential for Li 2 O. For the Catlow potential, we chose the first model of the three proposed. The difference between the first and second model is only a different set of parameter for the core-shell interactions, and the third model showed less accurate lattice parameters at higher temperatures in preliminary simulations. For the Oda potential for Li 2 O, we selected the FIT-GGA parameters, which offered the most accurate thermal expansion and melting point.
All the potentials we used predicted the Fm3m structure as the ground state at room pressure, and showed reasonable thermodynamical properties. However, no assessment regarding the quality of the potentials, or their ability to reproduce quantitatively given materials properties, was made otherwise.
Some of the potentials were designed with polarisable coreshell models to improve dielectric properties and elastic constants. However, the shell models tested tended to show instabilities and large fluctuations of the thermodynamical properties at high temperatures. We therefore ignored coreshell interactions when using these potentials, and considered only rigid ions. This should not affect the structure at low temperatures because of the way the interactions were fitted. In fact, the first step in the potential design process was to adjust the pair interactions to reproduce the structure correctly (Catlow and Norgett, 1973). Adjusting the parameters for the shell models was done in a following step, to improve dielectric constants. In any case, these potentials were not intended to reproduce very disordered structures or melts, therefore they should not be less accurate in principle than other, rigid-ions only potentials at high temperatures.
In addition, the Catlow potentials for BaF 2 , CaF 2 , and SrF 2 , as well as the Evangelakis potential for CaF 2 showed some instability at high temperatures, due to the relatively low energy barrier in their F-F Buckingham interactions. To avoid issues during the simulations with ions falling in the Buckingham singularity at short separation distances, a Ziegler-Biersack-Littmark (ZBL) contribution was added (Ziegler et al., 1985). We used the switching scheme implemented in LAMMPS with an outer radius of 2.0 Å and inner radii of 1.5 Å and 1.6 Å for the Catlow and Evangelakis potentials respectively. The Catlow potential for α-PbF 2 was derived independently and did not suffer from the same issue, and was therefore not modified.
Considering this diversity of compounds and potentials, it is expected to see some variance in the properties and structures of each specific material depending on the potential used. However, this diversity is also expected to limit any bias in our results by separating effects shown with only some of the potentials from common features independent of potentials details. Our intent was to identify features and qualitative behaviour universal in crystals with the Fm3m structure, rather than finding the absolute most accurate potentials. For this reason, we privileged a larger number of potentials rather than a more limited selection. This means that our results should be interpreted in general in terms of trends across compounds and potentials, and not quantitative predictions of materials properties. We will point out when a given potential deviates significantly from experimental results, or from the behaviour of other potentials.

Molecular Dynamics Simulations
The input structures for the MD simulations were 10 × 10 × 10 supercells of the conventional Fm3m unit cell. This resulted in cubic simulation boxes with 12 ,000 atoms and sides around 50 Å long, depending on potential and temperature. With all the potentials, electrostatic interactions were calculated using the Wolf sum with a damping parameter of 0.3 Å -1 (Wolf et al., 1999). All pair interactions, including electrostatics, were truncated and shifted using a cutoff radius of 11 Å. Equilibrium MD calculations were carried out from room temperature to the temperature at which only liquid phases were observed. This upper bound for the temperature range depends on both composition and potential. To sample these ranges similarly for each potentials, 50 evenly-spaced temperatures were selected in each case. The code used was LAMMPS (Plimpton, 1995), with both temperature and pressure controlled by the time-reversible Tuckerman integrator (Tuckerman et al., 2006) based on the Nosé-Hoover (Nosé, 1984) method and the Parrinello-Rahman strain energy (Parrinello and Rahman, 1981). When structure optimisation was useful, it was done by minimising the potential energy of the structure whilst keeping the simulation box fixed, using the implementation of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm from the L-BFGS library (Byrd et al., 1995).
Each simulation consisted of a several successive relaxation steps. The first one was a constant-temperature, fixed-volume (NVT) equilibration for 10 ps in order to thermalise the simulation box. This was followed by a fixed-temperature and room-pressure (NPT) equilibration for 240 ps, during which the average box lengths were calculated. Both cell size and shape were allowed to fluctuate during this step in order to avoid constraining the relaxation process. Although the box shape could change during this step, it remained very close to cubic on average except in simulations of liquids. Then, the box shape was fixed to match the average shape calculated during the NPT step. Following this, each supercell underwent a last NVT relaxation for 10 ps. Finally, data was accumulated over a 20 ps simulation at constant volume and energy (NVE). The time step for all these simulations was 1 fs, and the damping time for the thermostat and barostat were 0.01 and 0.5 ps respectively. The enthalpy H(T) was determined for each temperature T in K by taking its average over the NVE runs. The constant-pressure heat capacity C P (T) was then determined by differentiating the enthalpy This was done numerically using a second-order centred finite difference scheme. The lattice parameter a(T) was calculated by taking the average of the box size along the x, y, and z axes and dividing it by 10, the number of conventional cells along each direction. The linear thermal expansion coefficient was calculated using the measured lattice parameter, with the same numerical scheme employed to calculate C P (T). Eq. 2 was used instead of the equivalent and more common expression α 1 a za zT for the increased numerical precision of its finite differences discrete form. Statistical sampling of the phase space is important, particularly when considering highly disordered structures. For this reason, all simulations were run three times independently, with identical initial atomic positions but different initial velocities. The figures presented here are averages over these three simulations, but figures showing the complete data are available in the Supplementary Material for this article.
Diffusion being an important aspect of the superionic transition, we also calculated diffusion coefficients from the MD simulations. This was done by using the Einstein relation where t is the simulation time, and u is the atomic displacement. The mean squared displacements 〈|u(t)| 2 〉 were calculated for both A and X species separately, to produce specific diffusion coefficients. Point defects are an important structural feature at high temperatures. Since they play an important role in diffusion and have been linked to the onset of the superionic transition (Hiernaut et al., 1993), it is useful to be able to determine how common they are. Given that the structure considered here are stoichiometric to ensure electric charge neutrality and that the number of atoms is fixed and does not change over the course of a simulation, the only possible defects are Frenkel pairs, where an ion leaves its initial lattice site to form an interstitial elsewhere in the crystal. Counting the number of Frenkel pairs requires counting either the number of vacancies or the number of interstitials. We did so by mapping each ion in structures of interest taken from MD simulations to a reference crystal, which was a Fm3m structure with the same composition and lattice parameter. Thus, there are three possibilities for each lattice site. It can be: (i) empty, indicating the presence of a vacancy; (ii) occupied by one ion; or (iii) occupied by two ions, which indicates the presence of an interstitial. This analysis has been performed separately for the A (FCC) and X (SC) sublattices. It can result in principle in under-counting defects compared to experimental methods, if some ions are displaced enough to leave their site but not enough to be mapped onto another one. We also did not consider the interstitial octahedral site as a separate site for the purpose of the defects analysis.
We performed additional structure analysis using the polyhedral template matching (PTM) (Larsen et al., 2016) technique implemented in Ovito (Stukowski, 2010). This method assigns a kind of symmetry environment to each atom by comparing the relative positions of its neighbours to reference polyhedra for prototype structures such as simple cubic, hexagonal close packed, body-centred cubic, or face-centred cubic. This was used to show the structural disorder across the superionic transition. The RMSD cutoff for the PTM analysis was set to 0.5 for all structures.

Diffraction Patterns
Powder diffraction patterns were generated using the Debyer code, which implements the Debye scattering equations (Farrow and Billinge, 2009; Wojdyr, (https://github.com/wojdyr/debyer)). We use the copper Kα 1 wavelength λ 1.54056 Å for these calculations. In order to increase the signal-to-noise ratio, the 12 ,000-atom fluorite cells were first duplicated 7 times in all three directions, resulting in 4,116 000-atoms supercells. To investigate separately the structures of the A and X sublattices, additional supercells were set up by copying the full supercells and then removing either the cations or the anions. Thus, for each simulation we calculated a full pattern accounting for all the atoms, plus partial patterns representing respectively the structures of the A and X sublattices. A sinc damping function was applied to the radial distribution functions. This, in combination with the supercells, is helpful to limit artefacts from the Fourier transform, which would otherwise be large if the 12 ,000-atoms boxes had been used instead. We set the cutoffs for the Fourier transforms to less than one half of the lengths of the supercells. The intensity was finally plotted in powder XRD patterns as a function of the scattering angle 2θ.
Two-dimensional single crystal diffraction patterns were generated using the approach described in (Jin et al., 2020). So called reciprocal space maps (RSMs) were computed in the HHL planes, where H and L are non-integer multiples of the reciprocal lattice unit cell of a given fluorite structure 1/a AX 2 (T), where a AX 2 (T) is the lattice parameter of the AX 2 crystal at the considered temperature T. In order to increase the signal-tonoise ratio the HHL planes have a "thickness" of 1 4 reciprocal lattice unit cell over which the intensity was integrated. The RSMs were computed with the 12 ,000-atoms cells which results in a visible broadening of the reciprocal lattice points, and the presence of finite box size interference fringes, especially visible at low disorders. In order to quantify the intensities of the different reflections observed in the RSMs, line scans (with a 0.2 width in HKL units) were first extracted from the 2D data. These scans were further modelled with pseudo -Voigt functionsi.e. a linear superposition of a Gaussian and a Lorentzian function-to represent the Bragg peaks, an asymmetrical Gaussian function to represent the isotropic diffuse scattering and an additional linear background. The intensities were obtained after a conventional least-square fitting of the model to the 1D scans.

Thermodynamical Properties
As a first step, the general behaviour of the potentials was assessed by calculating their thermodynamical properties as a function of temperature: enthalpy H, constant-volume heat capacity C P , and thermal expansion coefficient α. These were used to verify that all potentials predicted a superionic transition for their compound, and to estimate the transition temperature T S . Some room-temperature properties calculated using the empirical potentials are shown in table 2. The constantpressure heat capacities are very similar across the potentials, and close to the value expected from the Dulong-Petit law of 2 | Thermophysical properties at room temperature for the different potentials, compared to experimental references: (1) Wicks and Block (1963); (2) Roberts and White (1986); (3) Dworkin and Bredig (1968); (4) Johnston and Bauer (1951); (5) Hull et al. (1988); (6) Popov et al. (2017); (7) Schröter and Nölting (1980); (8) Smith et al. (1963); (9) Pavlov et al. (2017); (10) Martin (1988).

Compound
Potential C P /10 -4 eV/(atom K) α/10 -5 K -1 2.59 · 10 -4 eV/(atom K). This indicates a strong harmonic behaviour from all the potentials, consistently with experimental observations. The main discrepancy with experiments is with Li 2 O, which has a significantly smaller heat capacity. The linear thermal expansion coefficient is more potential-dependent than the heat capacity. In general, the models tend to over-estimate α, sometimes significantly. For example, the Cazorla and Sayle potentials for BaF 2 , as well as the Bingham and Sayle potentials for CaF 2 , the Gillan potential for SrCl 2 , and the Bingham potential for SrF 2 have errors between 40 and 50%. This is not surprising, considering that α is very sensitive to the details of the potentials, and that the parameters are almost never fitted to it. Overall, the values for both C P and α show that all the potentials behave well at room temperature, even though some are more accurate than others for these properties. The simulation were performed in a temperature range from 300 K to the temperature T M at which the initial crystal turned into a liquid. This transition does not correspond to the experimental melting point T m , and is in general higher as shown in Table 2. By definition, the melting point measured experimentally is the temperature at which the liquid becomes thermodynamically more favourable than the crystal. However, this does not mean that the crystal is strictly unstable: it can be metastable instead and require climbing a non-negligible energy barrier to melt. The small length-scale and short times considered in our simulations inhibits nucleation of a liquid phase, which is the first step of a usual melting transition, by preventing the structure from climbing that energy barrier. The solid to liquid transition we observe is more analogous to a mechanical melting mechanism, and indicate the temperature at which the crystal becomes unstable rather than metastable (Wang et al., 1997). Other methods, such as the so-called moving interface technique (Morris et al., 1994), would be required to measure the thermodynamical melting points T m accurately. For this reason and to avoid confusion, the thermodynamical (experimental) and mechanical melting points will be noted T m and T M respectively throughout this article.
The most visible manifestation of the superionic transition from an experimental point of view is the Λ-shaped peak in C P , or equivalently the inflection point in the enthalpy H(T). Such peaks have been reported for several fluoride (Hull, 2004), chloride (Schröter and Nölting, 1980;Hull et al., 2011), and oxide (Ralph, 1987) compounds with the Fm3m structure. Amongst the potentials considered here, some, such as the Pedone potential for Li 2 O shown in Figure 2, showed a well-defined, sharp peak as generally expected. It was not as clear for other potentials such as the Gillan potential for SrCl 2 , the Cazorla potential for BaF 2 , or the Sayle potential for BaF 2 , which is also shown in Figure 2.
These potentials actually showed a plateau instead of a peak, and the superionic transition is not very visible when looking at thermodynamical properties for such potentials. The other potentials show intermediate behaviours, with C P peaks with various degrees of sharpness. Two potentials for the same compound can have different C P characteristics. An example of this is BaF 2 , where the Catlow potential has a very sharp peak, whilst no peak is visible with the Cazorla potential. Figures equivalent to Figure 2 for the other potentials can be found in the Supplementary Material for this article. A similar behaviour, also with a sharp Λ-shaped peak, was observed in the linear thermal expansion coefficient α(T). Interestingly, the potentials showing a strong C P peak similarly tend to show a clear α peak at the same temperature. From the C P (T) and α(T) curves, we can define four qualitatively different temperature ranges: i) a low-temperature regime, in which the heat capacity varies slowly as a polynomial; ii) an exponential increase leading to the top of the Λ peak and the superionic transition, attributed to the formation of Frenkel defects (Szwarc, 1969); iii) a decrease forming the right-side of the Λ peak; iv) a high-temperature regime in which C p varies slowly, leading to T M .
Some potentials show only some of these features. For instance, regime (i) can be very short for potentials with a low superionic transition temperature. The potential in which the C P peak is not visible would not have a clear stage (iii), instead going directly from (ii) to (iv). For the potentials that show all four stages, the C P peak is the simplest way of measuring the superionic transition temperature T S , which is by convention the temperature corresponding to the top of the peak. This is naturally more difficult to measure for potentials that show a plateau or no peak at all.
The potentials show significant variation in their predictions for T S , with differences of ∼300 K between the lowest and highest values for each compound. This highlights the variety of purposes for which the potentials were designed. More recent potentials tend to be more accurate, presumably thanks to more sophisticated potentials fitting techniques. Potentials that were intended to reproduce the superionic transition, such as Oda for Li 2 O, are unsurprisingly the most accurate. Overall, the trend across materials is reproduced, with SrCl 2 and UO 2 having respectively a lower and a higher T S on average compared to the other materials. The main outlier is the Catlow potential for β-PbF 2 , which overestimates T S by 700 K. However, few other potentials are available for this compound, which makes the discussion of any potential-dependent effect difficult. It should also be noted that Pb 2+ ions in PbF 2 are quite polarisable, which  (9) 3,126 (10 Frontiers in Chemistry | www.frontiersin.org October 2021 | Volume 9 | Article 723507 could point to a limitation of the rigid ions models to simulate this material. This point has been made in the literature (Castiglione and Madden, 2001), where calculations were done using shell models, with results closer to experimental references. The mechanical melting temperatures T M are not available experimentally, therefore those obtained from our simulations cannot be compared directly with references. We could still verify that they were greater than the reported melting temperatures T m for each compound.
The transition temperature T S has been linked to the melting point T m . The values for both T S and T M from our simulations are summarised in Table 3 and compared with experimental values. There is some variation across potentials, but overall there is a qualitative agreement for T S , i.e. compounds with a high transition temperature experimentally also tend to have higher values using empirical potentials. It has been reported that the TS Tm ratio is generally around 0.8 (Hiernaut et al., 1993). However, this is not very accurate, as experimental references summarised in Table 3 show T S T m ratios between 0.63 (for β-PbF 2 ) and 0.85 (for SrF 2 and CaF 2 ). Whilst we did not measure the thermodynamical melting points T m that would be directly comparable with experiments, we did measure the mechanical melting points T M . The T S T M ratios from our simulations are similarly scattered. They are potential-dependent, and range from 0.57 (for CaF 2 with the Evangelakis potential) to 0.79 (for SrCl 2 with the Bendall potential).

Point Defects and Diffusion
We calculated the fraction of point defects as a function of temperature for each potential. In general, we found as expected that the fraction of A defects was negligible in most of the considered temperature ranges (Catlow and Norgett (1973)). In all cases, the number of A Frenkel pairs just below T M was still two orders of magnitude smaller than that of X defects. No defect at all could be seen on the A sublattice for temperatures lower than about 100 K below T S , therefore they will be ignored and the following discussion is focused on the X defects.
Populations of X defects were significant at higher temperatures with all the potentials we considered. Figure 2D shows this for both BaF 2 and Li 2 O with the Sayle and Pedone potentials respectively. The simulations showed a pattern with a transition from a low-temperature Arrhenius regime, in which the X defect fractions takes the form where H f is the formation enthalpy of an X Frenkel pair, to a hightemperature regime in which the number of defects grows with a smaller effective formation energy. The effective formation enthalpy below the superionic transition H f is shown for all potentials in Table 4. The difference between both regimes and the inflection point can be seen in Figure 2D. The T S temperature marks the start of the deviation of the number of defects calculated from MD form its low-temperature fit. This links the change in defect properties to the onset of the superionic transition. The same pattern was observed for all potentials, with some variation in the magnitude of the transition. In general, potentials showing a sharp C P peak also had a large difference between the low-temperature and high-temperature behaviours.
Compound Potential H f /eV

This work Exp
BaF 2 Catlow and Norgett (1973) 3.67 1.6 (1) , 1.81 (2)  Beyond the superionic transition, the X sublattice is expected to be highly disordered. However, defects are detected using the Fm3m structure as a reference, which may be inadequate in presence of large distortions and possible extended defects. For this reason, the number of defects is much less reliable in the superionic phase, and no attempt was made at obtaining a quantitative measurement of the number of point defects in this case. Crystals in the Fm3m structure are known to exhibit three different diffusion behaviours (Voronin and Volkov, 2001). The diffusion is extrinsic when it is limited by the availability of extrinsic defects to provide diffusion pathways. It becomes intrinsic, with a higher activation energy, when the temperature is high enough to provide a significant population of Frenkel pairs. Diffusion is then mediated by point defects hopping, and depends on the availability of X Frenkel pairs. Finally, in the superionic phase, diffusion is faster and is thought to involve collective displacements (Annamareddy and Eapen, 2017), with a lower activation energy than the intrinsic Frenkel pair mechanism (Voronin, 1995;Voronin and Volkov, 2001). Whilst our simulations could not reproduce extrinsic diffusion mechanisms because of the lack of extrinsic defects, we found intrinsic and superionic behaviours consistent with what could be expected from the literature. Both behaviours are visible in Figure 2C. We fitted exponential functions of the form D(T) D 0 e −Ea/kBT to the calculated diffusion coefficients on both sides of the superionic transition to obtain diffusion activation energies E c a and E s a for the perfect crystal and the superionic phase respectively. The resulting activation energies are shown in Table 5. All the potentials considered here showed an activation energy change during the superionic transition, with E c a > E s a , albeit to different degrees. The crossover between the low-and high-temperature behaviours does not occur at T S for the potentials for which T S could be calculated from the Λ peak in the C P (T) curves. Instead, as was the case for the number of defects, T S marks the point at which the diffusion coefficients from MD start deviating from the low-temperature fit. The switch to the high-temperature behaviour itself occurs over a temperature range of a few hundred Kelvins. The upper bound of this range corresponds to the transition between regimes (iii) and (iv) as described in Section 3.1, i.e., to the high-temperature end of the Λ peak. For potentials that do not show a clear peak in C P (T), the change in diffusion coefficients can be used to estimate T S . This method is less accurate, however, because the deviation of D(T) to the fit is more difficult to characterise than a clear peak of C P . All the transitions temperatures are summarised in Table 3.  (2001) Frontiers in Chemistry | www.frontiersin.org October 2021 | Volume 9 | Article 723507 The activation energies in the Fm3m phases E c a from the MD simulation agree with the experimental references in general.
The fact that all potentials showed two different regimes for diffusion of the X chemical species regardless of the presence of a peak in C P and α confirms that the superionic transition can hapen without resulting in Λ peaks. This means that it is not inconceivable that such a transition could also happen without Λ peak in real materials. In fact, there are reports of this in SrCl 2 (Dent et al., 2004).

Simulated Powder XRD Patterns
As mentioned previously, the analysis of the structural evolution of AX 2 fluorites with temperature by a direct analysis of the MD simulations proved very challenging, particularly in the superionic phases. An approach to avoid this issue is to visualise the information contained in MD cells in a more compact representation. This is for instance the case when considering the radial distribution functions (RDFs). In the present study we choose to compute powder XRD patterns, i.e., Fourier transforms of RDFs. They are more sensitive to long-range atomic order and its alteration (Debelle et al., 2014), and are therefore a very useful complement to realspace structure analysis. We illustrate our findings with the example of BaF 2 with the Sayle empirical potential, bearing in mind that all conclusions drawn below are qualitatively valid for all the other potentials. This particular case was selected firstly because it exhibits XRD features common to all the composition-potential combinations investigated, with clearly visible temperature-dependent evolutions. It also has a plateau in C P and α instead of a Λ peak. It is hence both a representative example of all configurations investigated in this work, and a demonstration of a superionic transition without a Λ peak.
The evolution of powder XRD patterns as a function of temperature is illustrated in Figure 3 (2θ-T figure), considering all atoms in the MD cell ( Figure 3A) and only the X sublattice ( Figure 3B). We will refer to them as full-crystal XRD and X-only XRD patterns, respectively. Additionally, we report selected full XRD and X-only XRD patterns at two temperatures on Figure 4.
Let us first consider the full XRD patterns shown in Figures 3,  4 A. At 300 K, the characteristic peaks of the fluorite unit cell are clearly visible and can be indexed within the Fm3m symmetry group. Unsurprisingly, these peaks gradually shift towards smaller angles with increasing temperature, following the thermal expansion of the material. This is more pronounced at higher angles in agreement with the derivative of Bragg's law FIGURE 4 | Simulated powder XRD patterns from MD simulations of BaF 2 using the Sayle potential (Sayle et al., 2003), at temperatures below and above the superionic transition: (A) full crystal; (B) F sublattice only. Peak indices refer to the full Fm3m structure. The superionic transition temperature for this potential is 1,300 K.
ΔT being the temperature change.
The intensity of the peaks also decreases with increasing temperature, as a result of increasing local disorder, i.e. increasing thermal motion. The decrease in intensity is quantified by the Debye-Waller factor (Trueblood et al., 1996), which, in the case of harmonic and isotropic vibrations, is expressed by the simple expression T(θ) exp −8π 2 〈u 2 〉sin 2 θ/λ 2 (6) where u is the displacement of the atom from its regular lattice site. The effect of intensity damping is more pronounced at higher angles, as was the case for the peak shift. Finally, it can be noticed that a background of diffuse scattering develops around the most intense peaks. The intensity of this background increases up to the melting point where only diffuse scattering remains, but with a significantly redistributed intensity. The origin of this diffuse scattering background lies in the presence of thermal motion also responsible for the above-mentioned intensity damping, but it can also reveal the presence of structural defects.
In all the cases we investigated, the A-specific pattern was virtually identical to the full-crystal pattern. For this reason, the A patterns will not be discussed here, but 2θ − T maps can be found int the Supplementary Material. The similarity between full-crystal and A-only patterns can be simply understood by considering the difference in atomic numbers Z between light elements constituting the X sublattice (F, Cl, O), as compared to the heavier elements occupying the A sublattice (Ca, Ba, Pb, U). The latter dominate the overall XRD signal the intensity of which is proportional to Z 2 , thus making the contribution from lighter elements fainter.
Let us now consider the X-only XRD patterns shown in Figures 3, 4 B. Besides the similar general features related to thermal expansion (peak shift) and thermal motion (decrease in intensity) there are two noticeable differences with the full-crystal patterns. The first difference is the positions and intensities of the peaks. This is easily explained by the fact that the X sublattice has a SC symmetry, as opposed to the FCC symmetry of the full fluorite structure. The lattice parameter of the SC sublattice is half of that of the overall Fm3m structure. To avoid any ambiguity, we index the peaks of the X-only patterns using the same Miller indices as the full fluorite structure. This means that, for instance, the peak labeled 200 would correspond to the 100 peak of a SC unit cell of the X sublattice.
The second difference between X and full-crystal powder XRD patterns is the diffuse scattering background. Its intensity relative to the Bragg peaks is much higher, which indicates a much higher level of disorder on the X sublattice than on the A sublattice. The diffuse scattering is also more broadly distributed across the 2θ range, roughly centred around the 200 peak. This is a sign of the presence of amorphous or highly disordered local environments with a first neighbour distance approximately given by the X-X distance in the corresponding sublattice.
An additional 111 peak (or 1 2 1 2 1 2 in a base SC cell) emerges at temperatures around T S in the X-only patterns. This peak could not be detected in the full XRD pattern because its position exactly coincides with the 111 peak of the fluorite structure. The appearance of a superstructure peak with halved Miller indices points to a doubling of the periodicity of the unit cell of the X sublattice, which becomes equal to that of the overall Fm3m crystal. However, the lack of other superstructure peaks such as 113 and the weak intensity of the 111 peak makes it difficult to draw firm conclusions from those XRD patterns. For this reason, we now consider single-crystal diffraction.

Simulated Single-Crystal XRD Patterns
Single-crystal XRD provides another way of understanding the structural changes related to the superionic transition. Similarly to powder XRD, because it operates in the reciprocal space, it can reveal long-range features that are difficult to observe directly in simulation boxes due to thermal noise and increasing disorder at high temperatures. However, in contrast with powder XRD, there is no orientation averaging involved in the calculation. Hence, it is possible to precisely select definite regions of the reciprocal space to analyse. In order to be able to detect the appearance of superstructure peaks of the type 111, 113, etc., we computed the intensity distribution in (HHL) reciprocal lattice planes (that is, with a [110] zone axis). As for the powder diffraction study, we focus on BaF 2 with the Sayle potential and we use Miller indices relative to the Fm3m structure even when discussing the X sublattice. Data showing the evolution of the RSMs of all investigated cases are available in the Supplementary Material. We first briefly discuss the RSMs obtained from the A sublattice shown in Figures 5A-D. From room temperature up to T S , besides diffuse scattering concentrated around the Bragg peaks, no significant changes were observed. The peak coordinates have here been corrected for thermal expansion to allow for an easier comparison between different temperatures. At T M , a diffuse scattering halo characteristic of an amorphous structure is observed pointing to the melting of the material ( Figure 5D). These observations lead to two important conclusions. Firstly, apart from the observed thermal diffuse scattering, the fact that the signal from the A sublattice remains constant throughout the process -with no visible peak shift, peak broadening or splitting, or additional diffuse scattering-seems to indicate that the A sublattice is free of structural defects. This is in stark contrast with the X sublattice, which we will discuss shortly. We also know from direct analysis of the MD simulations shown in Section 3.2 that these defects are rare up until T M . Secondly, the thermal diffuse scattering forms streaks along the main 〈111〉 directions, indicating that atomic displacements are more pronounced in the close-packed {111} planes. This is also consistent with the weak and inhomogeneous diffuse scattering observed in the powder XRD data.
As far as thermal diffuse scattering is concerned, similar conclusions can be drawn for the X sublattice shown in Figures 5E-H. Diffuse scattering streaks appear along the {001} and {110} planes of the SC structure, pointing to enhanced thermal motion in the corresponding lattice planes. However, as soon as the temperature increases, major changes are observed on the X sublattice. The first difference is the appearance of weak intensity high-order 33L superstructure peaks (for instance at 722 K, indicated by the arrows on Figure 5F). The second noticeable difference is the increase of the diffuse scattering connecting the Bragg peaks, together with the appearance of a weak diffuse scattering halo with [002] radius, indicated on the figure by a dotted circle. The halo is barely visible in the RSM in this case, but it can be definitely detected by analysing line scans performed along the [111] direction (see also Figure 6). The RSMs of all cases investigated are given in the Supplementary Material. All of them exhibit the diffuse halo appearing almost concomitantly with the 33L peaks with varying intensity. In the superionic phase ( Figure 5G, representing a FIGURE 6 | Evolution of the intensity of the additional single-crystal XRD peaks in the superionic phase of BaF 2 using the Sayle potential: (A) evolution of the [111] intensity distribution as a function of temperature; (B) evolution of the 111 (blue circles), 222 (black circles), 333 (red circles) and diffuse scattering (green circles) peaks as a function of temperature. All intensities are on a logarithmic scale and normalised with respect to the room temperature of the 222 peak. Lines represent the evolutions of FP (thick blue), isolated HCP environments (red dotted) and clustered HCP environments (dotted blue). The different coloured areas correspond to the different regimes (i)-(iv) mentioned in the text.
Frontiers in Chemistry | www.frontiersin.org October 2021 | Volume 9 | Article 723507 structure at 1,472 K), additional 11L superstructure peaks are formed and the diffuse halo is now clearly visible. It should be noted that the apparent ellipsoidal shape of the halo is due to the partial superposition of the halo with the diffuse streaks along [110] and [001]. The 33L peaks are still visible, but their intensity decreases because thermal motion starts affecting high angle reflections. Finally, all Bragg peaks disappear at T M . Only the diffuse halo remains then, characteristic of an amorphous structure, in this case of the molten phase. All the potentials showed both the additional peaks and the diffuse halo. However, the relative intensity of these features is potential-or materialdependent. For instance, all the potential for Li 2 O showed a particularly intense diffuse halo, which is explained by the fact that the light element Li is significantly affected by thermal motion. Some potentials also produced fainter 33L and 11L peaks than others. All these features are visible in the animations, which can be found in the Supplementary Material. This evolution is summarised in Figure 6, which depicts the changes of the [111] line scans and the different intensities with temperature. There, it can be clearly observed that the diffuse halo starts to appear at a temperature slightly higher than room temperature (∼ 400-500 K, green circles), hence pointing to thermal vibrations. It then steadily increases with temperature. The 333 peak appears almost concomitantly (∼ 500 K, red circles) and reaches a maximum close to T S . Around 850 K, the 111 peak starts to increase significantly (blue circles). This corresponds to the end of stage (i) and the beginning of stage (ii). The end of stage (ii), marking the superionic transition correspond to a temperature where the intensity of the 111 peak exhibit a change in slope and where the intensity of the 333 peak starts decreasing. Finally at T M , all intensities drop to 0 except the diffuse scattering, which remains constant.
This evolution is consistent with the powder XRD data, in particular the appearance of a 111 peak at T S and a broad diffuse scattering background in X-only patterns. However, because of the orientation averaging inherent in powder XRD and the consecutive peak overlap, it is not possible to detect low intensity details as those discussed here. Finally, the fact that the diffuse scattering steadily increases until the melting temperature without any discontinuity, reinforces the description of the X substructure as a liquid made in the literature. In the next section, we discuss the origin of the superstructure peaks. Visual examination of the structures during MD simulations at different temperatures reveals some structural changes. To highlight them, we extracted slices from supercells from MD simulations of BaF 2 using the Sayle potential. The slices are 3 Å thick and centred on {100} atomic planes containing only F ions. At room temperature, the potential predicts a perfect Fm3m crystal. The slice at 300 K ( Figure 7A) indeed shows a square structure, corresponding to the {100} faces of the cubes forming the X sublattice. At intermediate temperatures, such as 722 K as in Figure 7B, some local distortions are visible, showing isolated triangular features. Visualisation of these features is helped by the use of the polyhedral template matching technique, which provides a characterisation of the local environment around each particle. At this temperature, almost all of the F ions are still in their ideal (SC) local environment. The closest matching structure for the ions with non-cubic local environments is HCP. The apparition of HCP local environments has been verified for all the potentials. With most of them, it coincides with the end of regime i and beginning of regime (ii). Most of the time, it also seemed to happen before Frenkel pairs were detected. However, this should be treated with caution. Indeed, the methods we used to count Frenkel pairs, though qualitatively accurate, can result in an under-estimation of the number of defects. On the other hand, the PTM analysis could also be too sensitive, resulting in false positives.
This changes around the superionic transition temperature T S . Indeed, in the superionic phase of regimes (iii) and (iv), we could see some clustering of the atoms with HCP environments, which were much more likely to be first neighbours. This is shown in Figure 7C, where we can see clusters of HCP environments mixed with cubic regions. This pattern holds until the mechanical melting point T M , after which the structure becomes liquid ( Figure 7D).
The ions with unknown local environments are cases where the PTM algorithm could not decide how to categorise the particle. There are very few of them below T S . For example, in the case of BaF 2 with the Sayle potential, they represent 5% of all F ions at the superionic transition. They are more prevalent in the superionic phase, but still amount to 13% at the mechanical melting point. They are much more present in the liquid phase, where the structure is not well-suited to the PTM analysis. These ions with disordered local environments are expected to contribute to the diffuse scattering observed in both powder and single-crystal diffraction patterns. The MD simulation boxes in the superionic phases also show that there can be a high degree of misorientation between the different types of local environments. This also contributes to the diffuse scattering.
This evolution of the number of ions with HCP local environments is plotted in Figure 6, together with the XRD data. This figure clearly confirms the findings described previously. Isolated HCP environments appear at relatively low temperatures (∼ 600 K), well before Frenkel pairs. The concentration of those isolated HCP environments, noted hcp i , are neatly correlated with the apparition of the 33L superstructure peaks, which indicate that these defects reduce the unit-cell symmetry to a period of ∼ 1 3 of the (111) fluorite lattice spacing. In terms of temperature range, the fraction of clustered HCP environments, noted hcp c , is almost perfectly correlated both with the fraction of Frenkel pairs, and the intensity of the 111 peak. This indicates that the Frenkel pairs are at the origin of (i) the clustering of HCP environments, and (ii) the formation of ordered local environments on the X sublattice, with a periodicity equal to that of a fluorite unit-cell.
The PTM analysis in itself does not give a complete picture of the local environments it detects. For example, the actual local environments also depend on the position of neighbours on the A sublattice. Those were completely ignored in our PTM analysis, which was done separately for each sublattice. Determination of the local environments directly from MD simulations proved all but impossible due to distortions and large thermal fluctuations in the superionic phases. To work around this, we performed constant-volume geometry optimisations from the final structure of the MD simulations for each potential and temperature. During this operation, the ions were moved towards more energetically favorable positions, thus removing thermal noise from vibrations and unstable defects. The results were structures that were easier to interpret and characterise. It should be noted, however, that the energy-minimised and high-temperature structures are not the same. Indeed, some temperatureactivated phenomena disappear during geometry optimisation. Thus, the energy-minimised structure can only provide some clues as to how to interpret high-temperature MD simulations.
The structures optimised from high-temperature simulations showed both cubic and non-cubic local environments after PTM analysis. The trend was also similar to what was observed in hightemperature simulations. In regime i, all potentials produced a perfect Fm3m crystal after optimisation. This changed in regime (ii), where we could observe some isolated ions with non-cubic environments, as well as some trapped Frenkel pairs, which could not recombine during energy minimisation. Above the superionic transition temperature, most potentials showed a separation into two phases, which contained the cubic and non-cubic environments respectively. Visual examination of the noncubic phase showed the same triangular pattern observed in MD simulations, except that the pattern covered extensive regions of the simulation box instead of small, elongated clusters. This is visible in the inset in Figure 7C. In addition, because all thermal noise was removed, we could isolate the noncubic environment for a more thorough characterisation. We found that the structure we isolated had the symmetry elements of the Pbcn space group. This space group corresponds to the scrutinyite structure (α-PbO 2 ), and is well-known in some compounds of the AX 2 form. In particular, there are some instances of such compounds undergoing a pressure-induced phase transition to distorted fluorite structures Léger and Haines (1997); López-Moreno et al. (2016). A transition from Fm3m to Pbcn hasn't been documented in materials with the fluorite structure, however, except for UO 2 , where it was seen during fracture (Zhang et al. (2012)) or under tensile mechanical loading (Fossati et al. (2013)). Even in that case, the Pbcn phase was not thermodynamically stable.
In the Pbcn structure, the X sublattice has a distorted HCP structure, whilst the A sublattice remains close to FCC, as shown in Figure 8B. The pattern shown in the (0001) basal plane of the X sublattice in that structure, shown in Figure 8, is the triangular pattern visible in both MD simulations and energy-minimised structures. Thus, we can conclude that the non-cubic local environments detected by the PTM analysis are actually X ions which have a local environment consistent with that of the Pbcn structure. The overall structure of the superionic phase, however, is quite different from what was observed during the Fm3m → Pbcn transition (Fossati et al., 2013). Indeed, it does not consist in a clear separation between a Fm3m and Pbcn phases. It is formed instead of a mixture of local environments with different structures, along with some ions with a disordered environment. Moreover, we can see on Figure 7C that there can be a large degree of disorientation between different local environments.
We can now propose a mechanism for the superionic transition in materials with the fluorite structure. Upon heating, the X sublattice is subject to thermal vibrations compatible the with symmetry of the crystal. This gives rise to the diffuse scattering in the {001} planes observed in the RSMs of the X sublattice. There is a known vibration mode that fits this description. It is the B 1u mode with a [0 0 1 2 ] wavevector, which has been linked to the superionic transition in CeO 2 (Klarbring et al., 2018). The polarisation vector of this mode corresponds exactly to the displacement of the X ions at the beginning of the Fm3m → Pbcn transition shown in Figure 8C. These distortions result in the appearance of local environments with HCP symmetry elements, initially isolated and distributed in the Fm3m structure. Upon further heating, the crystal starts forming X Frenkel pairs. This coincides with the clustering of isolated HCP local environments, with symmetry elements characteristic of the Pbcn structure, throughout the crystal. These local environments are responsible for the anomalous enthalpy, and the exponential increase in the heat capacity leading to the Λ peak. The structure change could also be the reason of the different diffusion mechanism in the superionic phase. In contrast, the volume fraction of the disordered (liquid-or amorphous-like structure), responsible for the diffuse scattering halo, steadily increases from room temperature up to melting without any noticeable change. It is hence unlikely to be involved in the drastic changes observed in the thermodynamic properties.

CONCLUSION
We used a wide selection of empirical potentials to simulate several different compounds with the fluorite structure and investigate structural changes between the low-temperature perfect structure and that of the superionic phase. Although a full validation of each potential was not done here, we verified that the properties they predict are realistic by comparing them to available experimental values. In particular, the superionic transition temperatures, the enthalpies of formation for the X Frenkel pairs, and the activation energies for X diffusion in both the Fm3m and superionic phases were found to be adequate. Therefore, whilst the potentials are quantitatively imperfect when taken separately, we are confident that the qualitative trends they show are not potential dependent.
From a thermodynamical point of view, not all the potentials reproduced the characteristic C P and α peaks usually associated with the superionic transition, despite showing other aspects such as structural disorder and a change in the diffusion properties. This demonstrates that the transition could occur without a clearly-visible peak. Different potentials for the same compound could also exhibit different behaviour. For example, in BaF 2 , the Catlow potential showed a sharp Λ peak in its C P (T) curve, whereas the Sayle potentials showed a broad plateau, and the transition was almost invisible from the thermodynamical properties of the Cazorla potential.
All potentials are characterised by a 2-regime diffusion behaviour, with a change in the apparent activation energy at the superionic transition consistent with available experimental data. Moreover, populations of X Frenkel pairs showed a similar trend, with a low-temperature Arrhenius regime followed by a progressive transition towards a lower apparent formation energy in the superionic phases. This is consistent with a transition from a diffusion process mediated by isolated point defects in lowtemperature crystals to a collective process involving clusters in the superionic phases.
Powder and single-crystal XRD patterns were calculated from the MD simulation boxes for investigate structural changes from the reciprocal space. Whilst the A sublattices did not change from room temperature to the melting point, additional Bragg peaks on the X sublattices showed the emergence of new structural features distinct from the simple partially molten structure sometimes assumed. This observation could not be made using experiments, which could not have separated signals from both sublattices. Subsequent analysis of the MD simulation boxes could determine that the additional peaks are the result from a the presence of local environments with a structure close to the Pbcn structure. In these local environments, the X sublattice has HCP features. This work is a demonstration of the power of combined MD and computational diffraction techniques to described structures associated with widespread disorder.
Following our observations, we can describe the evolution with temperature of the structure of the fluorite compounds, following the four stages defined from thermodynamical properties and outlined in Section 3.1. Only the X sublattice is discussed in the following: the A sublattice remains FCC from room temperature to T M , with some distortions at high temperatures.

Stage (i)
This is the low-temperature regime, in which the heat capacity varies slowly as a polynomial. Near the end of this temperature range, isolated hexagonal local environments are formed in an otherwise perfect Fm3m matrix. This stage is the extrinsic diffusion regime. No diffusion could be measured from the MD simulations because for their short time scale and lack of extrinsic defects.

Stage (ii)
In this regime, C P (T) follows an exponential increase leading to the top of the Λ peak, or the left-hand edge of the plateau depending on the potential. The end of this stage is T S and the onset of the superionic transition. The structure is characterised by the accumulation of Frenkel pairs and hexagonal environments becoming more present. The number of isolated HCP environments peaks near the end of this stage. This corresponds to the intrinsic diffusion regime, limited by Frenkel pairs, where diffusion coefficients follow an Arrhenius law.

Stage (iii)
In this transition stage, the C P (T) curve decreases sharply, forming the high-temperature side of the Λ peak. The number of isolated hexagonal environments also decreases, to become negligible by the end of this stage, at which point all the hexagonal local environments are part of clusters. The beginning of this stage is associated with a change in slope of the Frenkel pair concentration, diffusion coefficient, and the concentration of clustered hexagonal local environments, which all increase at a reduced rate. The diffusion coefficients progressively change from the stage (ii) to the stage iv Arrhenius regimes. This stage might not be visible if the potential does not predict a C P peak.

Stage (iv)
The growth of the C P (T) curve resumes slowly, leading to the mechanical melting at T M . The structure is characterised by a mixture of regions with hexagonal local environments, and regions that retain the cubic local environment of the perfect fluorite structure. The X diffusion coefficient follows another Arrhenius law, with an activation energy smaller than that of stage (ii), underlying the collective diffusion mechanism underlines elsewhere (Annamareddy et al., 2014).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.