Cooperative excitations in superionic PbF2

Chris Erik Mohn, Marcin Krynski, Walter Kob, Neil L. Allan

► To cite this version:

Chris Erik Mohn, Marcin Krynski, Walter Kob, Neil L. Allan. Cooperative excitations in superionic PbF2. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Math. or Phys. Character (1896-1934), 2021, 379 (2211), pp.20190455. 10.1098/rsta.2019.0455. hal-03440257

HAL Id: hal-03440257

https://hal.science/hal-03440257v1

Submitted on 11 Mar 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mohn, C. E., Krynski, M., Kob, W., & Allan, N. L. (2021). Cooperative excitations in superionic PbF₂. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 379(2211). https://doi.org/10.1098/rsta.2019.0455

Peer reviewed version

Link to published version (if available):
10.1098/rsta.2019.0455

Link to publication record on the Bristol Research Portal
PDF-document

This is the accepted author manuscript (AAM). The final published version (version of record) is available online
via The Royal Society at https://doi.org/10.1098/rsta.2019.0455. Please refer to any applicable terms of use of
the publisher.

**University of Bristol – Bristol Research Portal**

**General rights**

This document is made available in accordance with publisher policies. Please cite only the
published version using the reference above. Full terms of use are available:
http://www.bristol.ac.uk/red/research-policy/pure/user-guides/brp-terms/
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

CrossMark
click for updates

**Subject Areas:**

physical chemistry

**Keywords:**

β-PbF₂, superionic conductivity, DFT, collective diffusion, molecular dynamics

**Author for correspondence:**

Chris E. Mohn
e-mail: chrism@geo.uio.no

Cooperative excitations in
superionic PbF₂

Chris E. Mohn<fnref n="1" />, Marcin Krynski<fnref n="2" />, Walter

Kob<fnref n="3" />, Neil L. Allan<fnref n="4" />

Links between dynamical Frenkel defects and collective diffusion of fluorides in β-PbF₂ are explored using Born-Oppenheimer molecular dynamics. The calculated self diffusion coefficient and ionic conductivity are 3.2 × 10⁻⁵ cm<fnref n="2" /> s⁻<fnref n="1" /> and 2.4Ω⁻¹cm⁻<fnref n="1" /> at 1000 K in excellent agreement with results from pulsed field gradient and conductivity measurements. The calculated ratio of the tracer diffusion coefficient and the conductivity diffusion coefficient (the Haven ration), is slightly less than unity (about 0.85), which in previous work has been interpreted as providing little evidence for collective "multi-ion" diffusion. In contrast, our molecular dynamics simulations show that fluoride diffusion is highly collective. Analysis of different mechanisms shows a preference for direct collinear "kick out" chains where a fluoride enters an occupied tetrahedral hole/cavity and pushes the resident fluoride out of its cavity. Jumps into an occupied cavity leave behind a vacancy, thereby forming dynamic Frenkel defects which trigger a chain of migrating fluorides assisted by local relaxations of the lead ions to accommodate these chains. The calculated lifetime of the Frenkel defects and the collective chains is around ~1 ps in good agreement with that found from neutron diffraction.

<footnote n="1">¹ Centre for Earth Evolution and Dynamics, University of Oslo, N-0315 Oslo, Norway</footnote>

<footnote n="2">² Faculty of Physics, Warsaw University of Technology, ul. Koszykowa 75, 00-662 Warsaw, Poland</footnote>

<footnote n="3">³ Laboratoire Charles Coulomb, UMR 5221, University of Montpellier and CNRS, 34095 Montpellier, France</footnote>

<footnote n="4">⁴ School of Chemistry, University of Bristol, Bristol BS8 1TS, UK</footnote>
# 1. Introduction

Fast-ion conductivity in the solid state was first observed by Faraday himself in the early 18th century when the resistivity of a solid piece of lead fluoride, PbF₂, suddenly dropped when heated to high temperature [1,2]. This opened up the field of solid-state ionics and the β phase of PbF₂ is now widely considered a parent of a broad class of “type II” superionics that on heating undergo a rapid and *continuous* increase in the ionic conductivity [1,2]. These compounds include many fluorite structured halides (XF₂ with X= Ca, Sr, Ba) as well as functional oxides such as CeO₂, UO₂, ThO₂ and ZrO₂ [2]. In contrast to “type I” superionics (e.g. δ-Bi₂O₃), where the transition to the superionic phase is first order due to disordering of the inherent vacancies, β-PbF₂ contains no such vacancies, and the superionic transition is associated with transient and extended Frenkel defects [2–4]. The ionic conductivity in β-PbF₂ increases continuously by eight orders of magnitude when heated from room temperature to the superionic transition temperature, T_c = 711 K [5], levelling off to ~3 Ω⁻¹ cm⁻¹ just above T_c [5,6], where it remains nearly constant until the melting point at 1158 K. The β-phase of PbF₂, thus, possesses one of the highest known values of ionic conductivity at intermediate temperature, more than 3 orders of magnitude higher than that of stabilized zirconia frequently used as the electrolyte in Solid Oxide Fuel Cells (SOFC). The exceptionally high conductivity makes β-PbF₂ particularly attractive as a reference for designing the next generation of fuel cells and gas separation membranes to lower the current high operating temperature of these devices. There is, however, at present very little atomistic understanding of the nature of the Frenkel defects in β-PbF₂ and how they trigger diffusion processes. Such knowledge should give crucial input into the design and development of new fast ionic conductors in future devices and links between structure and conductivity therefore warrant further attention. Here we use *ab initio* Born–Oppenheimer Molecular Dynamics simulations (MD) to investigate the nature of ionic transport within β-PbF₂ which allows us to learn more about collective (multi-ion) dynamics. In particular we aim to probe the links between the dynamical Frenkel defects, nature of fluoride migration and the role played by the immobile lattice in promoting such fast ion behavior.

The crystal structure of β-PbF₂ is shown in Fig. 1. It has a cubic fluorite structure (space group *Fm*3*m*) in which the lead ions form a face centered “immobile” cubic sublattice. The fluorides are located in the tetrahedral holes/cavities, but crystallographic interpretations from different neutron diffraction (ND) experiments report large variations in the average fluoride positions within these cavities (see e.g. Ref. [2]) reflecting that the extensive thermal “activity” of the mobile sub lattice is not so easy to capture using simple average crystallographic models. From these ND experiments, different models of the dynamic and highly transient Frenkel defects have been suggested. The best fit to ND data was found using a defect model with as many as 9 vacant cavities with an equal amount of fluorides significantly displaced from their regular tetrahedral site [3]. These extended Frenkel clusters are highly transient with lifetime ~10⁻¹² s which is on the magnitude of a jump-time of the fluorides [5]. This in turn, indicates that formation of Frenkel defects is intimately linked to - or are manifestations of - collective migrating fluorides.

Atomistic insight into the nature of diffusion can be obtained from experiment by the comparison of the *single particle* self diffusion or tracer diffusion coefficient (i.e. *D*<sup>self</sup> or *D*<sup>*</sup>) and the diffusion coefficient calculated (indirectly) from conductivity experiment, *D*<sup>σ</sup>, which carries information about the *net* displacement of all mobile species. * The self diffusion coefficient can be found from different techniques such as nuclear magnetic resonance (NMR) which requires some *a priori* knowledge of the mechanism of diffusion and also more *directly* from pulsed field gradient (PFG) experiments. The ratio of tracer diffusion coefficient to the conductivity diffusion coefficient is the Haven ratio, *H* = *D*<sup>*</sup>/*D*<sup>σ</sup>. Using *D*<sup>self</sup> from Ref. [7] and *D*<sup>σ</sup> from conductivity measurement [6], a value of *H* between 0.7 and 0.9 (at around 1000 K) has been taken as evidence for little collective diffusion or for the presence of a simple *non-collinear* interstitialy

<footnote>*Here we can assume that the tracer diffusion coefficient and the self diffusion coefficient have the same value since isotope effects in tracer diffusion experiment can be ignored.</footnote>
mechanism [7] shown in Fig. 2a. Here an interstitial fluoride, F1, at the octahedral 4b site of the fluorite structure kicks one of its nearest neighbours, F2, at the regular tetrahedral 8c site, into a new octahedral site and replaces the atom that has been kicked out. Such a mechanism is consistent with *H* being slightly less than 1, but appears to be inconsistent with results from ND which shows that the fluorides do not occupy the interstitial 4b site [3]. Other possible mechanisms, not considered previously, are “kick out” cooperative events in which a fluoride, F1, at the regular tetrahedral 8c site of the cubic fluorite structure enters a neighboring occupied tetrahedral cavity and pushes the fluoride in that cavity, F2, into a new cavity. This new cavity, C, is normally occupied but can also be “instantaneously” empty if a fluoride has just left it. If C is empty, as shown in Fig. 2c, the collective chains of fluorides is terminated by the jump into cavity C. By contrast, if C is occupied, as shown in Fig. 2f, longer “caterpillar” chains may form. Short collinear strings such as that visualized in Fig. 2c, have been shown to play an important role in understanding fast ion conduction in e.g. δ-Bi₂O₃ [8] and could also contribute to collective diffusion in β-PbF₂. Indeed with longer multi-ion chains a combination of many different mechanisms can contribute to the conductivity in superionic PbF₂. In addition to the “caterpillar” mechanism mentioned above (Fig. 2f), “rings” (Fig. 2g), for example, can also add to the conductivity. Interestingly, an eight-ion ring similar to that shown in Fig. 2g, has been suggested to contributes strongly to diffusion in β-PbF₂ [9], but such closed “loops” can not explain its high value of conductivity alone since the mechanism involve no net charge transport! In principle, combinations of different mechanisms are consistent with the value of the Haven ratio which is just an *average* value of different possible collective paths. As long as the *net* mean square displacement of the cooperative diffusing chain of fluorides is longer than that of uncorrelated diffusion, the Haven ratio is less than unity.

## 2. Theory

To investigate the collective dynamics of β-PbF₂ from the atomic Born–Oppenheimer MD trajectories we first calculate the ionic conductivity by integrating the charge-current correlation function

$$ \sigma = \frac{1}{k_B T V} \int_0^\infty J(t) dt, \qquad (2.1) $$

where V is the unit-cell volume, $J(t) = e^2 \sum_{i,j} Z_i Z_j \langle \mathbf{v}_i(0) \cdot \mathbf{v}_j(t) \rangle$, $Z_i$ is the ionic charge number, $\mathbf{v}_i$ is the velocity of ion i respectively and $k_B$ is the Boltzmann constant. The conductivity diffusion coefficient, $D^\sigma$, can then be calculated from Eq. 2.1 using the Nernst-Einstein relation

$$ D^\sigma = \frac{k_B T}{\rho (e Z^-)^2} \sigma. \qquad (2.2) $$

Following Castiglione and Madden [10] we also calculate the conductivity in the form of a mean-square displacement of the centre of mass:

$$ \sigma = \frac{e^2}{k_B T V} \lim_{t \to \infty} \frac{1}{6t} \langle \left| \sum_i Z_i^+ \delta \mathbf{r}_i^+ (t) + \sum_i Z_i^- \delta \mathbf{r}_i^- (t) \right|^2 \rangle. \qquad (2.3) $$

Here the summations within the brackets are the charge weighted *net* displacements of all cations and anions respectively. The advantage of using Eq. 2.3 rather than Eq. 2.1 is that the integral over the charge current correlation function converges slowly due to the non-negligible contribution of the tails in $J(t)$. Comparison of results from Eq. 2.1 and Eq. 2.3 can thus help us to quantify the effect of correlations due to these tails in $J(t)$.

To calculate *H* we also need $D^{\text{self}}$ which can be obtained from Eq. 2.3 by ignoring correlations in the atomic positions between the fluorides. We can then write down the well known Nernst-Einstein equation

$$ \sigma^{\text{self}} = \frac{e^2}{k_B T} \lim_{t \to \infty} \frac{1}{6t} \left( \langle \rho^+ | Z^+ \delta \mathbf{r}^+ (t) |^2 \rangle + \langle \rho^- | Z^- \delta \mathbf{r}^- (t) |^2 \rangle \right), \qquad (2.4) $$
where $\rho^+$ and $\rho^-$ are the densities of the cations and anions respectively.

### (a) Single particle and collective dynamics from a jump diffusion model

We can also calculate $D^{\text{self}}$ and hence $\sigma^{\text{self}}$ from a jump diffusion model as follows: Since vibrations and jump frequencies occur at different time scales, we can decompose the ionic trajectory into vibrations and jump events. We divide the simulation cell into space-filling primitive cubes [11–14] with the cube-center located at the tetrahedral 8c site of the fluorite structure, and the octahedral 4b site and the lead atoms at alternate cube-corner positions. This division allows us to calculate the jump frequency by counting the number of times an ion leaves a primitive cube and hence a tetrahedral cavity and enter a new one. It should be noted that the fluorides vibrate with large amplitudes that sometimes extend into a neighboring fluoride tetrahedral cavity. To avoid counting such vibrations as two consecutive jumps (jump relaxations) we require that the fluoride must stay within the new cavity for at least $\tau^{\text{thresh}}$ before we label the event as a jump. The value of $\tau^{\text{thresh}}$ should be at least on the size of a low frequent vibrational mode, but not much larger either because then two consecutive rapid jumps could be identified as a single long jump. We find a reasonable choice of $\tau^{\text{thresh}} = 1$ ps which we confirm below by the comparison of the diffusion coefficient calculated from a jump diffusion model (using this value of $\tau^{\text{thresh}}$) with that calculated from the MSD.

We also record correlations between two successive jumps of a single fluoride and distinguish between back-jumps (jump-relaxations) which we define as movements in which an atom jumps back to where it came from. Sideway jumps involve two jumps in different crystallographic directions (i.e. a jump in the $\langle 100 \rangle$ direction is followed by a jump in either the $\langle 010 \rangle$ or the $\langle 001 \rangle$ directions), whereas forward jumps are two consecutive jumps in the same crystallographic direction i.e. a jump in the $\langle 100 \rangle$ direction is followed by a jump in the $\langle 100 \rangle$ direction. We then calculate the self-diffusion coefficient using

$$D^{\text{self}} = \frac{1}{6} \sum_{i} f_{i} \Gamma_{i} a_{i}^{2}, \qquad (2.5)$$

where $i$ refers to the jump directions i.e. between nearest neighbours cavities aligned in the $\langle 100 \rangle$ direction, next-nearest neighbours aligned in the $\langle 110 \rangle$ direction or third-nearest neighbours aligned in the $\langle 111 \rangle$ direction. $\Gamma_i$ and $a_i$ are the jump frequency and the jump distance in direction $i$. Correlations factors, $f_i$, can be calculated using $f_i = 1 + 2\langle\cos\theta_{l,l+1}\rangle$, where $\theta_{l,l+1}$ is the angle between two consecutive jumps $l$ and $l+1$.

We also investigate collective diffusion from a hopping model where correlation between diffusing neighbouring pairs are identified. We label two nearest neighbouring cavities A and B and assume that these cavities are occupied by fluorides F1 and F2. We identify events where F1 jumps to B and where F2 jumps to a third cavity C. If the cavities A, B and C involved in the diffusion of the two fluorides are all aligned in the same direction as shown in Fig 2c, the mechanism is collinear. If the three cavities are not aligned along the same crystallographic direction (as shown in Fig. 2e), the mechanism is non-collinear whereas if A is identical to C, the fluorides exchange positions which gives no net charge transport (see Fig. 2b). Longer chains of collective diffusing fluorides can form a range of different patterns, such as the caterpillar mechanism (Fig. 2f), arcs and also loops (Fig. 2g).

### (b) Single particle and collective dynamics from dynamic heterogeneities

Analysing the nature of dynamical heterogeneity of the highly mobile species is another possible approach to understand the nature of particle dynamics and, in contrast to the hopping model, there is no reference to the underlying crystallographic lattice. Instead, a "rearrangement indicator" $R_i$ [15,16] of a fluoride $i$, is used to measure the variance in the fluoride positions
within a time window $\delta t$ as:

$$R_i = \sqrt{\langle \langle (\mathbf{r}_i - \langle \mathbf{r}_i \rangle_{t-\delta t/2})^2 \rangle_{t+\delta t/2} \langle \langle (\mathbf{r}_i - \langle \mathbf{r}_i \rangle_{t+\delta t/2})^2 \rangle_{t-\delta t/2}} \quad (2.6)$$

where $\langle \rangle_{t+\delta t/2}$ and $\langle \rangle_{t-\delta t/2}$ denotes appropriate time averages. The value of $\delta t$ is typically of the size of the jump time (similar to $\tau^{thresh}$ used in the jump diffusion model described in the previous subsection). After calculating the “rearrangement indicators” of the fluorides along the trajectory, we next identify the fluorides that are probably hopping as those with $\sqrt{R_i} > 0.95 a_{nn}$, where $a_{nn}$ is the distance between nearest neighbouring tetrahedral 8c positions of Fm$\overline{3}$m in the simulation box ($a_{nn} = 3.05$ Å). These “rearranging fluorides” can cluster together in space possibly capturing collective chains of migrating fluorides. We define the size of a cluster as the total number of rearranging fluorides where each ion in the cluster is closer than a given cutoff value, $d_{max}$ to at least one other fluoride within the same cluster. Note that if we do not find any other rearranging fluorides within $d_{max}$, of a given rearranging fluoride, the cluster contains only this fluoride (i.e. the cluster size is 1).

### 3. Computational details

Ab initio Born-Oppenheimer MD simulations were carried out using Density Functional Theory (DFT) in the NVT ensemble with a cubic $4 \times 4 \times 4$ simulation box containing 768 atoms with $V = (24.4 \AA)^3$ and $T = 1000$ K. The time step was 2.5 fs and a Nosé-Hoover thermostat was used. During a 120 ps simulation we counted more than 10 000 fluoride jumps, which ensured sufficient statistics for a qualitative comparison with the ionic (d.c.) conductivity obtained experimentally. Such a long trajectory is required to converge the centre of mass diffusion to sufficient precision since only one point (i.e. the centre of mass of the ions) is collected at each time step. The d.c conductivity was thus calculated from the collective MSD after the ballistic regime from $t = 2$ ps to $t = 4$ ps. At each step along the trajectory the total energy was minimized using the Perdew-Burke-Ernzerhof functional [17], an energy cut-off of 300 eV and sampling the Brillouin zone at the gamma point only, as implemented in VASP [17–19].

### 4. Results and discussion

After the ballistic regime, we calculate the single particle and collective diffusion coefficients from the slopes of the MSD (shown in Fig. 3). The calculated $D^{self} = 3.2 \times 10^{-5}$ cm$^2$ s$^{-1}$ at 1000 K is in good agreement with result from a PFG experiment ($D^{PFG} = 4 \times 10^{-5}$ cm$^2$ s$^{-1}$) [7]. From the collective MSD we find that $D^{d.c.}$ is $3.6 \times 10^{-5}$ cm$^2$ s$^{-1}$ which is also in good agreement with the value extracted from conductivity measurements using Eq. 2.2. Our calculated conductivity from Eq. 2.3, $\sigma^{d.c.} = 2.4 \Omega^{-1}$ cm$^{-1}$ which is in good agreement with that measured experimentally $\sigma^{d.c.} \approx 3 \Omega^{-1}$ cm$^{-1}$ at around 950 K [5], and with that from a MD study at 1000 K using pair potentials [20].

The Haven ratio is $0.85 \pm 0.01$ when Eq. 2.3 is used to calculate the collective diffusion coefficient and $0.80 \pm 0.09$ when Eq. 2.1 is used to calculate $D^{\sigma}$. The good agreement between these two methods indicates that our MD run is sufficiently long to ensure convergence of the integral in Eq. 2.1. Our Haven ratio is markedly higher than the Haven ratio estimated from MD calculations using pair potentials [10], but in very good agreement with that reported by Carr et al [6], in which the self diffusion coefficient from the PFG experiment [7] and $D^{\sigma}$ extracted from conductivity measurement, gave a Haven ratio between 0.7 and 0.9 at around 1000 K [6]. Bearing in mind the modest precision in the PFG data and also the difficulty of measuring the conductivity above 950 K due to sublimation losses of the sample [6], the overall agreement is satisfactory.

Our calculated Haven ratio is lower than 1 but larger than those of other (type 1) superionic conductors with a high concentration of inherent vacancies such as $\delta$-Bi$_2$O$_3$ ($H \approx 0.5$ [8]), indicating that possibly collinear events are less important in type II superionic conductors, such
as β-PbF₂, than in type I. To provide some mechanistic insight into the diffusion of the fluorides we now carry out an atomistic investigation of the mechanism of diffusion in β-PbF₂.

### (a) Single particle dynamics

Analysis of the *ab initio* MD trajectories at 1000 K using a jump diffusion model shows that about 80% of the fluorides migrate to the nearest neighbouring cavity aligned in the ⟨100⟩ direction, about 17% migrate to the next nearest neighbours aligned in the ⟨110⟩ direction, whereas only about 4% jump in the ⟨111⟩ direction. This is remarkably consistent with MD results for CaF₂ where 80%, 16% and 3% of the jumps were in the ⟨100⟩, ⟨110⟩ and ⟨111⟩ directions respectively [21]. The fluorides vibrate mostly within the tetrahedral cavities centered at the 8c site, as seen in Fig. 4, where we plot the ionic density and two snapshots of the MD trajectory. The average residence time, τ<sup>residence</sup>, is about 5.1 ps and the jump times, τ<sup>jump</sup>, between the tetrahedral cavities is typically between 0.5 and 1 ps. The fluorides follow irregular and highly curved pathways as highlighted in Fig. 4 and the anions sometime resides, for a very short time, near an octahedral site or near the edge between two nearest neighbours tetrahedral cavities. The octahedral site, however, is *not* a resident site for the fluorides but a site avoided. If we attempt to quench dumped MD configurations along the trajectory by a steepest descent (structural optimisation), all fluorides in the simulation box end up at the regular tetrahedral site and none at the octahedral site. This confirms that the octahedral site is not a resident site for the fluorides. The lead ions also relax appreciable distances ~ 0.5 Å to accommodate the jumping fluoride. This highlights the important role played by the “soft” cation sublattice in promoting fast-ion conductivity.

We find only evidence for correlations between successive jumps when these are between nearest neigbours tetrahedral cavities. That is there is a slightly higher probability of back jumps (prob. ≈ 0.2) compared to that of two uncorrelated consecutive jumps (prob.≈ 0.167). Inserting *f*<sub>&lt;100&gt;</sub> = 0.167/0.200 = 0.83 and the calculated jump frequencies into Eq. 2.5, we find that *D*<sup>self</sup> = 2.6 × 10<sup>-5</sup> cm<sup>2</sup> s<sup>-1</sup> (with τ<sup>thresh</sup> = 1.0 ps) which is in very good agreement with *D*<sup>self</sup> calculated above from the single particle MSDs.

A preference for backwards jumps is consistent with previous experimental and computational results for other fast-ion conductors [22,23,26]. However, jump-relaxations seems to be less pronounced in compounds without inherent vacancies compared to other “type I” fast-ion conductors with a high concentration of vacancies such as δ-Bi<sub>2</sub>O<sub>3</sub> [23], α-CuI [13], Ba<sub>2</sub>In<sub>2</sub>O<sub>2.50</sub> [24], La<sub>2</sub>Ce<sub>2</sub>O<sub>7</sub> [25] etc.

### (b) Collective diffusion from a jump diffusion model and dynamical heterogeneities

An interstitialcy mechanism (see Fig. 2a) suggested in Ref. [7] is consistent with a Haven ratio slightly less than unity. However, since the interstitial 4b octahedral site is not a resident site for the fluoride, as is evident from the ionic density plot and the trajectories shown in Fig. 4, this mechanism cannot be used to explain the nature of collective dynamics within β-PbF₂. Instead, our MD trajectories show that the fluorides jump rapidly and directly between tetrahedral cavities. We start to analyse the nature of collective migrating fluorides from a hopping model by identifying all events where an fluoride, F1 in cavity A, jumps to cavity B (which is occupied by F2) and kicks/pushes F2 out of its resident cavity. The fate of F2 can be captured by classifying these “kick out” events as collinear (Fig. 2c), non-collinear (Fig. 2e) or exchange (Fig 2b). In this analysis we count all pair events where F1 and F2 reseides in the same cavity, B, for τ<sup>thresh</sup> = 0.1 ps, before we label the event as a pair-jump.

Compared to a random diffusion process, we find a marked preference for collinear “kick out” pairs (Fig. 2c) over non-collinear (Fig. 2e) and very few exchanges (Fig. 2b). Indeed about 40% of the “kick out” events are collinear (compared to 1/6 for uncorrelated diffusion), slightly less
than 60% are non-collinear ones (compared to 2/3 for uncorrelated diffusion) and less than 2%
exchange positions (compared to 1/6 for uncorrelated diffusion). This suggests that, together, a
preference for collinear "kick out" events and very few exchanges explain, at least in part, why the
Haven ratio is less than unity. This analysis does not include events where F1 and F2 do not share
a common cavity such as that shown in Fig 2e), i.e where a jump to cavity B can possibly "trigger"
a neighbouring fluoride, F1 in C to jump to D. Analysis of the trajectories shows that these events
are much less correlated than the "kick out" pair mechanism. The main contribution to collective
dynamics is therefore probably collinear "kick out" chains (Fig 2) but diffusion processes in PbF2
clearly involve a range of different mechanisms. Many of these are hard to capture in an jump
diffusion analysis since we cannot always distinguish jump relaxations and vibrations with large
amplitudes.

Difficulties in analysing the nature of collective chains with more than two atoms from a
hopping model due to large amplitudes of vibrations and local distortions of the sub-lattices
has motivated us to investigate the collective diffusion of fluorides from a cluster size analysis
using Eq. (2.6). The cluster sizes are plotted in Fig. 5 with different cutoff values, dmax, chosen
to capture different crystallographic alignments of fluorides involved in collective transport. If
we count clusters along the MD trajectory using a cutoff value of dmax = 1.1ann, we will only
include events where two nearest neighbouring fluorides are jumping at about the same time
i.e. the time difference between the two jumps is shorter than τthresh. Thus this includes "kick
out" events such as those shown in Fig. 2c and Fig. 2e, whereas events where the two fluorides
are aligned further apart i.e. as next-nearest neighbours before they jump, as shown in Fig. 2d, are
not counted. If dmax = 1.5ann in the cluster analysis, we also include many of the events displayed
in Fig. 2d and if we increase dmax further to 1.8ann, collective chains in which two fluorides are
third-nearest neighbours are also included in the analysis. With an even higher value of dmax
we include so many uncorrelated events that these start to outnumber the correlated events
included in the analysis. Thus the red and gray lines in Fig. 5 are therefore very similar when
dmax = 2.1ann. This provide an upper constraint for the correlation length of the collective chains
of migrating fluorides. When dmax = 1.1ann we find, compared to uncorrelated events, about 30%
less isolated jumps (cluster size = 1) and, two times more 3-atom chains. By contrast, clusters with
less than 4 or 5 atoms involving nearest neighbours (dmax = 1.1ann) are rare, which rules out the
long caterpillar-like mechanisms shown in Fig. 2c. When dmax = 1.1ann we find that the average
cluster size is 1.78 (compared to 1.38 for a random distribution of jumps) and when dmax = 1.5ann
the calculated average cluster size is about 2.45 (compared to 2.03 for random a distribution). In
Fig. 6, representative snapshots of the geometries of the collective chains of migrating fluorides
are shown. These illustrate that a range of different mechanisms are at play.

(c) Linking the nature of Frenkel defects and collective dynamics

Strings of collective migrating fluorides are intimately connected to the formation of short lived Frenkel defects seen in neutron diffraction [5]. A jump into an occupied cavity such as that in Fig 2c, leaves behind a vacant cavity and the migrating fluoride, F1, shares for a short time a cavity with F2, thereby creating a Frenkel defect. These Frenkel defects are highly dynamic. If we attempt to capture these in a local minimum configuration on the potential energy surface by quenching the MD configurations by a steepest descent (energy minimisation), the fluorides always end up in their regular tetrahedral sites. Formation of dynamic Frenkel defects distorts both the mobile and immobile lattices suggesting that these defects are highly extended in agreement with defect models suggested from neutron diffraction [3]. The average lifetime of a vacant cavity is about 0.5 ps, but about 5% of these are vacant for more than 1 ps. The lifetime of a doubly occupied cavity (where two fluorides are resident within the same tetrahedral cavity) is less than 0.4 ps. If we assume that the average lifetime of the vacancy is about the same as the lifetime of the collective fluoride chains, we roughly estimate the lifetime of the Frenkel defects to be around 1 - 2 ps. This is consistent with the lifetime of the highly dynamic Frenkel defects seen in neutron diffraction ≈ 1 ps [5]. These Frenkel defects are manifestations of collective diffusion in β-PbF2.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Crystal structure of β-PbF₂ (*F*m*3*m). The green and purple spheres are the fluorides and the lead ions located at the tetrahedral 8c and the 4a sites respectively, whereas the small black spheres are the (empty) octahedral 4b site. A tetrahedral hole, which we also refer to as a tetrahedron "cavity" in this work, is the tetrahedron formed by the four nearest neighbour lead ions at a 4a site surrounding a given tetrahedral 8cs site.</figcaption></figure>

**5. Conclusions**

Links between dynamical Frenkel defects and collective diffusion of fluorides in β-PbF₂ are investigated from a hopping model and from the nature of dynamical heterogeneities using *ab initio* molecular dynamics simulations. The ionic density plots and the MD trajectories show that the octahedral holes/cavities of the fluorite structure are empty and that the fluorides stay within the tetrahedral cavities for 5 ps, on average, before they jump rapidly within 1 ps to an occupied tetrahedral cavity. This can initiate a chain of migrating fluorides and a range of different collective mechanisms contributes to the conductivity. Analysis of correlations between pairs of simultaneously jumping fluorides shows that a collinear "kick out" mechanism (Fig 2c) occur twice as frequently as a random collinear kick out event and that direct exchange (Fig 2b) is rare. Other pair mechanisms investigated include that in Fig 2d where a jump into an occupied cavity can also trigger a nearest neighbour fluoride to jump. Such events, however, are much less correlated than the "kick out" events which also rules out correlations between two even more distant jumps that take place at about the same time. A combination of collinear kick out chains and few exchanges can thus explain why the Haven ratio is slightly lower than 1 [6]. Large local distortions of the "immobile" cation sub-lattice and rapid fluoride back-jumps make it difficult to capture all collective fluoride chains from a jump diffusion model in particular if these chains contain several fluorides.

This has motivated us to identify the nature of collective chains from the dynamical heterogeneity: We first calculate the "rearrangement indicator" of all fluorides within a short time interval using Eq. 2.6. We then identify the fluorides with a large rearrangement indicator typical
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Different possible mechanisms of collective diffusion in β-PbF₂ suggested in the literature and investigated in this work. Fig. a) is an example of a non-collinear interstitialcy mechanism suggested in Ref. [7]. Here a fluoride, F1, is assumed to sit at the interstitial octahedral site and then it migrates to the neighbouring tetrahedral site and kicks F2 into a new octahedral site. Fig b) shows a pair of fluorides exchanging positions. We also investigate different correlated events shown in c)-g). In c), e) and f) a fluoride, F1, jumps to another occupied cavity and kicks/pushes the resident fluoride, F2, out of its cavity. These “kick out” or “caterpillar” mechanisms can be collinear, as shown in c), or non-collinear as visualized in e). Fig. d) shows a fluoride, F1, that jumps to an instantaneously empty cavity B, which trigger a jump of a neighbouring atom F2 located at C to jump to D. The mechanisms in f) and g) are examples of longer collinear “caterpillar” like chains and “ring-like” pattern representing extreme examples of mechanisms with very low and very high Haven ratios respectively.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The collective (thick green line) and single-particle Nernst–Einstein (thin red line) mean-square displacements of the fluorides at 1000 K from an *ab initio* Born–Oppenheimer MD simulation.</figcaption></figure>

of a jump and calculate the distances between them. All rearranging fluorides within a given
cutoff distance that jump at about the same time (< 1 ps) are members of the same cluster. We find
that these clusters are more compact than those found by distributing the rearranging fluorides
at random over all fluoride positions, providing additional evidence for collective dynamics
within β-PbF₂. Compared to a mean of a random distribution, we find fewer isolated rearranging
atoms, twice as many 3-atom clusters and equally many clusters with 4 or more atoms when
dmax = 1.1ann.

The nature of collective diffusion in β-PbF₂ is strikingly similar to those seen in many
other fast-ion conductors [27,28,30–36]). Oxygen diffusion in δ-Bi₂O₃, fast Li⁺ diffusion in
Li₁₀GeP₂S₁₁₂ and silver ion diffusion in AgI, for example, all involve “multi-ion” chains with
H < 1. Collective diffusion in the fast-ion regime is advantageous since the nature of energetically
preferred local motifs and their connectivity patterns are kept intact as the chain of ions migrates
minimizing bond breaking and bond-formation. The “immobile” sublattice also plays important
roles as they relax appreciable distances away from their average lattices position generating
sufficient “local space” to accommodate the diffusing multi-ion chains. Thus, experimental and
atomistic studies are about to unveil the “universal” nature of diffusion in the superionic regime,
an important next step to design new devices for future energy technology.

**Data Accessibility.** The simulation output file (which also contains input parameters) can be downloaded
from https://www.doi.org/10.17635/lancaster/researchdata/

**Authors' Contributions.** CM performed calculations, analysis and wrote the paper. MK did analysis and contributed to discussion. NLA and WK contributed to analysis, discussions and writing of the paper. All authors read and approved the manuscript.

**Competing Interests.** The authors declare that they have no competing interests.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The top picture shows ionic density contours of the fluorides from an MD simulation at 1000 K. The bottom (left) picture shows trajectories of the fluorides (green lines) in a 5 ps time window and the bottom (right) picture shows two typical examples of two jump events. The jump of F1 follows a highly curved trajectory in the crystallographic <100> direction whereas F2 jumps several times including a jump in the <10¯1> direction.</figcaption></figure>

**Funding.** The Centre for Earth Evolution and Dynamics is funded by CoE-grant 223272 from the Research Council of Norway. The computing time was provided by the Norwegian metacenter for computational science (NOTUR) through a grant of computing time.

**Acknowledgements.** W.K. is a senior member of the Institut Universitaire de France.

## References

1. Faraday M 1838 Experimental researches in electricity—Twelfth series. *Phil. Trans. R. Soc.* 83-124
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The red full lines are the fraction of different cluster sizes containing fluorides that move distances of a typical jump i.e. with $\sqrt{R_i} > 0.95a_{nn}$ where $R_i$ is calculated from Eq. 2.6. The gray dashed lines are the corresponding values found from a mean of random distributions of rearranging fluorides over all fluoride positions.</figcaption></figure>

2. Hull S 2004 Superionics: crystal structures and conduction processes. *Rep. Prog. Phys.* **67** 1233-1314

3. Dickens MH, Hayes W, Hutchings MT, Smith C 1982 Investigation of anion disorder in PbF₂ at high-temperatures by neutron-diffraction. *J. Phys. C:Solid State Phys.* **15** 4043-4060

4. Goff JP, Hayes W, Hull S, Hutchings MT 1991 Neutron powder diffusion study of the fast-ion transition and specific heat anomaly in beta-lead fluoride. *J. Phys.:Condens. Matter* **3** 3677-3687

5. Derrington CE, Keeffe MO 1973 Anion Conductivity and Disorder in Lead Fluoride. *Nat. Phys.* *Sci* **246** 44-46

6. Carr VM, Chadwick AV, Saghafian R 1978 The electrical conductivity of PbF₂ and SrCl₂ crystals at high temperature. *J. Phys. C:Solid State Phys.* **11** L637-L641

7. Gordon RE, Strange JH 1978 NMR relaxation and self-diffusion in PbF₂. *J. Phys. C:Solid State Phys.* **11** 3213-3223

8. Mohn CE, Krinsky M 2020 Collective diffusion within the superionic regime of Bi₂O₃. *Phys. Rev. B* **101** 104309

9. Goodenough JB, 1984 Fast ionic conduction in solid. *Proc. R. Soc. Lond. A* **393** 215-234
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Different snapshots of fluoride positions during an MD run at 1000 K with 768 atoms in a cubic simulation box. The green spheres are fluorides with a rearrangement indicator, $\sqrt{R_i} > 0.95a_{nn}$, representing jumping fluorides (the arrow marks the jump-direction), whereas the gray are those with $\sqrt{R_i} < 0.95a_{nn}$ (i.e. the non migrating fluorides).</figcaption></figure>

10. Castiglione MJ, Madden PA 2001 Fluoride ion disorder and clustering in superionic PbF$_2$. *J. Phys. Chem.* **13** 9963-9983

11. Jacucci G, Rahman A 1978 Diffusion of F ions in CaF$_2$. *J. Chem. Phys.* **69** 4117-4125

12. Wengert S, Nesper R, Andreoni W, Parrinello M 1996 Ionic diffusion in a ternary superionic conductor: An ab initio molecular dynamics study. *Phys. Rev. Lett.* **77** 5083-5085

13. Mohn CE, Stølen S, Hull S 2009 Diffusion within alpha-CuI studied using ab initio molecular dynamics simulations. *S J. Phys.: Condens. Matter* **21** 335403

14. Mohn CE, Stølen S, Hull S 2012 Comment on "Dynamic Correlation between Superionic Coppers in alpha-CuI". *J. Phys. Soc. Japan* **81** 106001

15. Candelier R, Widmer-Cooper A, Kummerfeld JK, Dauchot O, Biroli G, Harrowell P, Reichman DR 2010 Spatiotemporal Hierarchy of Relaxation Events, Dynamical Heterogeneities, and Structural Reorganization in a Supercooled Liquid. *Phys. Rev. Lett.* **105** 135702

16. Schoenholz SS, Cubuk ED, Sussman DM, Kaxiras E, Liu AJ 2016 A structural approach to relaxation in glassy liquids. *Nat. Phys.* **12** 469-471

17. Perdew JP, Burke K, Ernzerhof M 1996 Generalized gradient approximation made simple. *Phys. Rev. Lett* **77** 3865-3868

18. Kresse G, Hafner J 1993 Ab initio molecular dynamics for open-shell transition metals. *Phys.*
Rev. B 47 558

19. Kresse G, Joubert D 1999 From ultrasoft pseudopotentials to the projector augmented-wave method. *Phys. Rev. B* **59** 1758-1775

20. López JD, Dioza JE, Correa H, 2019 Molecular dynamics simulation of ionic transport and thermodynamic properties in β-PbF₂. *ionics* **25** 5383-5390

21. Kaneko Y, Ueda A 1988. Molecular Dynamics Studies of Superionic Conductor CaF₂. *J. Phys. Soc. Japan* **57** 3064-3073

22. Mohn CE, Stølen S, Nordberg ST, Hull S 2009 Ab initio molecular dynamics simulations of oxide-ion disorder in the δ-Bi₂O₃. *Phys. Rev. B* **80** 024205

23. Mohn CE, Stølen S, Nordberg ST, Hull S 2009 Oxide-Ion Disorder Within the High Temperature delta Phase of Bi₂O₃. *Phys. Rev. Lett.* **102** 15550

24. Mohn CE, Allan, NL, Stølen, S 2006 Sr and Ga substituted Ba₂In₂O₅: Linking ionic conductivity and the potential energy surface. *Solid State Ionics* **177** 223-228

25. Kalland L, Norberg ST, Kyrklund J, Hull S, Eriksson SG, Norby T, Mohn CE, Knee CS 2016 C-type related order in the defective fluorites La₂Ce₂O₇ and Nd₂Ce₂O₇ studied by neutron scattering and ab initio MD simulations. *Phys. Chem. Chem. Phys.* **22** âÁR 24070-24080

26. Funke K 1993 Jump relaxation in solid electrolytes. *Prog. Solid State Chem.* **22** 111-195

27. Mohn CE, Allan NL, Freeman C, Ravindran P, Stølen S 2005 Collective ionic motion in oxide fast-ion-conductors. *Phys. Chem. Chem. Phys.* **6** 3052-3055

28. Shimojo F, Inoue T, Aniya M, Sugahara T, Miyata Y 2006 Ab initio molecular-dynamics study of static structure and bonding properties of molten AgI. *J. Phys. Soc. Jap.* **75** 114602

29. Kendrick E, Kendrick J, Knight KS, Islam MS 2007 Cooperative mechanisms of fast-ion conduction in gallium-based oxides with tetrahedral moieties. *Nat. Mat.* **6** 871-875

30. Xu M, Ding J, Ma E 2012 One-dimensional stringlike cooperative migration of lithium ions in an ultrafast ionic conductor. *Appl. Phys. Lett.* **101** 031901

31. Kalland LS, Mohn CE 2020 First principles calculations on order and disorder in La₂Ce₂O₇ and Nd₂Ce₂O₇. *Phys. Chem. Chem. Phys.* **22** 13930-13941

32. Tealdi C, Mustarelli P, Islam MS 2010 Layered LaSrGa₃O₇-Based Oxide-Ion Conductors: Cooperative Transport Mechanisms and Flexible Structures. *Adv. Funct. Mat* **20** 3874-3880

33. Hull S, Norberg ST, Eriksson SG, Mohn CE 2013 Neutron powder diffraction and molecular dynamics study of superionic SrBr₂. *J. Phys. Cond. Mater.* **25** 454205

34. Vasileiadis A, Wagemaker M 2017 Thermodynamics and Kinetics of Na-Ion Insertion into Hollandite-TiO₂ and O₃-Layered NaTiO₂: An Unexpected Link between Two Promising Anode Materials for Na-Ion Batteries. *Chem. Mater.* **29** 1076-1088

35. de Klerk NJJ, van der Maas E, Wagemaker M 2018 Analysis of Diffusion in Solid-State Electrolytes through MD Simulations, Improvement of the Li-Ion Conductivity in β-Li₃PS₄ as an Example. *Adv. Energy Mater.* **10** 2001486

36. Zou Z, Ma N, Wang A, Ran Y, Song T, Jiao Y, Liu J, Zhou H, Shi W, He B, Wang D, Li Y, Avdeev M, Shi S 2020 Relationships Between Na⁺ Distribution, Concerted Migration, and Diffusion Properties in Rhombohedral NASICON. *Adv. Ener. Mat.* **10** 2001486