Amangeldi Torayev

Luke Sperrin

Maria A Gomez

John A Kattirtzi

Ceĺine Merlet

Clare P Grey

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible Local Distortions and Dynamics in Hydrated Y-Doped BaZrO 3

Y-doped BaZrO 3 is a promising proton conductor for intermediate temperature solid oxide fuel cells. In this work, a combination of static DFT calculations and DFT based molecular dynamics (DFT-MD) was used to study proton conduction in this material. Geometry optimizations of 100 structures with a 12.5% dopant concentration allowed us to identify a clear correlation between the bending of the metal-oxygen-metal angle and the energies of the simulated cells. Depending on the type of bending, two configurations, designated as inward bending and outward bending, were defined. The results demonstrate that a larger bending decreases the energy and that the lowest energies are observed for structures combining inward bending with protons being close to the dopant atoms. These lowest energy structures are the ones with the strongest hydrogen bonds. DFT-MD simulations in cells with different yttrium distributions provide complementary microscopic information on proton diffusion as they capture the dynamic distortions of the lattice caused by thermal motion. A careful analysis of the proton jumps between different environments confirmed that the inward and outward bending states are relevant for the understanding of proton diffusion. Indeed, intra-octahedral jumps were shown to only occur starting from an outward configuration while the inward configuration seems to favor rotations around the oxygen. On average, in the DFT-MD simulations, the hydrogen bond lengths are shorter for the outward configuration which facilitates the intra-octahedral jumps. Diffusion coefficients and activation energies were also determined and compared to previous theoretical and experimental data, showing a good agreement with previous data measuring local proton motion.

INTRODUCTION

Perovskites have emerged as attractive electrolytes for intermediate-temperature solid oxide fuel cells (SOFCs) owing to their good ionic transport properties as well as adequate mechanical and chemical stabilities. [START_REF] Kreuer | Proton -Conducting Oxides[END_REF] In particular, Ydoped BaZrO 3 has been widely studied as this material shows very high protonic conductivities of up to 1 × 10 -2 S cm -1 at 450 °C when exposed to water vapor. [START_REF] Yamazaki | High Total Proton Conductivity in Large-Grained Yttrium-Doped Barium Zirconate[END_REF][START_REF] Bohn | Electrical Conductivity of the High-Temperature Proton Conductor BaZr 0.9 Y 0.1 O 2.95[END_REF] When BaZrO 3 is doped with yttrium, the charge imbalance between the yttrium and zirconium ions is accommodated by the generation of oxygen vacancies leading to the composition BaZr 1-x Y x O 3-x/2 . These oxygen vacancies are responsible for the moderate oxygen-anion conductivity at high temperatures. [START_REF] Bohn | Electrical Conductivity of the High-Temperature Proton Conductor BaZr 0.9 Y 0.1 O 2.95[END_REF] Under wet conditions, protons can be incorporated in the material through absorption of water molecules that dissociate forming hydroxide ions and filling the oxygen vacancies. The material can be fully hydrated, represented by the chemical formula BaZr 1-x Y x O 3-x (OH) x , or partially hydrated represented more generally as BaZr 1-x Y x O 3-x/2-y/2 (OH) y (y < x) with oxygen vacancies remaining. The high ionic conductivities achieved in wet conditions are promising for applications of this material in intermediate-temperature SOFCs. [START_REF] Bohn | Electrical Conductivity of the High-Temperature Proton Conductor BaZr 0.9 Y 0.1 O 2.95[END_REF] Y-substituted BaZrO 3 can be prepared with high concentrations of yttrium (up to 60% [START_REF] Fabbri | Does the Increase in Y-Dopant Concentration Improve the Proton Conductivity of BaZr 1-x Y x O 3-δ Fuel Cell Electrolytes?[END_REF] ), and a number of studies [START_REF] Fabbri | Does the Increase in Y-Dopant Concentration Improve the Proton Conductivity of BaZr 1-x Y x O 3-δ Fuel Cell Electrolytes?[END_REF][START_REF] Iguchi | The Relationship between Chemical Composition Distributions and Specific Grain Boundary Conductivity in Y-Doped BaZrO 3 Proton Conductors[END_REF] have reported a dependence of the ionic conductivity on the yttrium content. As a first guess, one would expect the ionic conductivity to increase with the doping concentration following the increase in the number of charge carriers. However, it was shown that an optimum conductivity is observed for a dopant concentration of around 10%-20%. [START_REF] Fabbri | Does the Increase in Y-Dopant Concentration Improve the Proton Conductivity of BaZr 1-x Y x O 3-δ Fuel Cell Electrolytes?[END_REF][START_REF] Iguchi | The Relationship between Chemical Composition Distributions and Specific Grain Boundary Conductivity in Y-Doped BaZrO 3 Proton Conductors[END_REF] In fact, many experimental and theoretical observations suggest that the protons are trapped close to the yttrium dopant sites increasing the activation barrier for diffusion and reducing the overall ionic conductivity. [START_REF] Larring | Protons in Rare Earth Oxides[END_REF][START_REF] Kreuer | Aspects of the Formation and Mobility of Protonic Charge Carriers and the Stability of Perovskite-Type Oxides[END_REF][START_REF] Yamazaki | Proton Trapping in Yttrium-Doped Barium Zirconate[END_REF] This proposed proton trapping has attracted considerable attention in recent years both from an experimental and from a theoretical point of view. [START_REF] Yamazaki | Proton Trapping in Yttrium-Doped Barium Zirconate[END_REF][START_REF] Pornprasertsuk | Proton Conductivity Studies of Y-Doped Barium Zirconate: Theoretical and Experimental Approaches[END_REF][START_REF] Merinov | Proton Diffusion Pathways and Rates in Y-Doped BaZrO 3 Solid Oxide Electrolyte from Quantum Mechanics[END_REF][START_REF] Bjorketun | Effect of Acceptor Dopants on the Proton Mobility in BaZrO 3 : A Density Functional Investigation[END_REF][START_REF] Gomez | The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO 3 , a Cubic Perovskite[END_REF][START_REF] Draber | Nanoscale Percolation in Doped BaZrO 3 for High Proton Mobility[END_REF][START_REF] Blanc | Dynamic Nuclear Polarization NMR of Low-γ Nuclei: Structural Insights into Hydrated Yttrium-Doped BaZrO 3[END_REF][START_REF] Friman | Molecular Dynamics Simulations of Proton Diffusion in Yttrium Doped Barium Zirconate[END_REF][START_REF] Raiteri | Reactive Force Field Simulation of Proton Diffusion in BaZrO 3 Using an Empirical Valence Bond Approach[END_REF][START_REF] Van Duin | ReaxFF Reactive Force Field for the Y-Doped BaZrO 3 Proton Conductor with Applications to Diffusion Rates for Multigranular Systems[END_REF] In hydrated Y-doped BaZrO 3 , proton conduction occurs through series of transfers (between neighboring oxygens) and rotations (around an oxygen). Static density functional theory calculations (DFT) and molecular dynamics simulations (based on DFT or reactive force fields) have been proposed as powerful methods to estimate energy barriers for the intraand interoctahedral transfers as well as for rotations. Static DFT calculations have shown that local distortions arise close to the dopant atom, usually leading to stronger hydrogen bonds, which could explain the trapping effect. [START_REF] Merinov | Proton Diffusion Pathways and Rates in Y-Doped BaZrO 3 Solid Oxide Electrolyte from Quantum Mechanics[END_REF][START_REF] Bjorketun | Effect of Acceptor Dopants on the Proton Mobility in BaZrO 3 : A Density Functional Investigation[END_REF][START_REF] Gomez | The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO 3 , a Cubic Perovskite[END_REF][START_REF] Draber | Nanoscale Percolation in Doped BaZrO 3 for High Proton Mobility[END_REF] However, local distortions also induce shorter O-O distances, which would tend to decrease the activation energies for proton jumping. It is worth noting that trapping does not necessarily imply that jumps do not happen but rather that protons spend a significant amount of time in the vicinity of the dopant and that the jumps are local rather than involving long-range transport. In many studies, local activation energies are combined with a jump-diffusion model to predict long-range energy barriers and diffusion coefficients. The resulting activation energies are usually in good agreement with experiments, and trends observed for different dopant atoms could be reproduced by Bjorketun et al. 11 In this context, Blanc et al. have defined specific proton environments in the perovskite material depending on the position of the proton with respect to the dopant atom (see Figure 1). [START_REF] Blanc | Dynamic Nuclear Polarization NMR of Low-γ Nuclei: Structural Insights into Hydrated Yttrium-Doped BaZrO 3[END_REF] In their DFT calculations, the configurations with protons close to the yttrium atoms indeed correspond to the lowest energy structures.

While static DFT calculations are very useful to identify proton diffusion pathways and the effect of specific local distortions, thermal fluctuations are always present at the temperatures used in the experiments (200-700 °C). These dynamic distortions may alter the picture obtained from a static point of view. DFT-MD simulations have the advantage of being able to capture both proton hopping and thermal fluctuations that occur on short time scales as demonstrated by many studies for water in the bulk liquid phase and at interfaces. [START_REF] Marx | The Nature of the Hydrated Excess Proton in Water[END_REF][START_REF] Hassanali | On the Recombination of Hydronium and Hydroxide Ions in Water[END_REF][START_REF] Kattirtzi | Microscopic Dynamics of Charge Separation at the Aqueous Electrochemical Interface[END_REF] More recently, Fronzi et al. applied DFT-MD simulations on proton conduction in undoped BaZrO 3 and reported that proton diffusion is enhanced when compressive strain is applied on the structure. [START_REF] Fronzi | First-Principles Molecular Dynamics Simulations of Proton Diffusion in Cubic BaZrO 3 Perovskite under Strain Conditions[END_REF] The authors interpret this increase in diffusion to the shorter O-O distances observed in the compressed material. The effect of strain on proton diffusion has also been demonstrated in Y-doped BaZrO 3 for two dopant concentrations (1% and 12.5%) using MD simulations with a reactive force field based on an Empirical Valence Bond description. [START_REF] Ottochian | Influence of Isotropic and Biaxial Strain on Proton Conduction in Y-Doped BaZrO 3 : A Reactive Molecular Dynamics Study[END_REF] In these dynamical studies and other related ones, [START_REF] Friman | Molecular Dynamics Simulations of Proton Diffusion in Yttrium Doped Barium Zirconate[END_REF][START_REF] Raiteri | Reactive Force Field Simulation of Proton Diffusion in BaZrO 3 Using an Empirical Valence Bond Approach[END_REF][START_REF] Van Duin | ReaxFF Reactive Force Field for the Y-Doped BaZrO 3 Proton Conductor with Applications to Diffusion Rates for Multigranular Systems[END_REF] no strong trapping is however observed. Raiteri et al. suggest that, for low dopant amounts (less than 15%), the trapping probability is low and thus not often observed, or hard to characterize, from MD simulations. [START_REF] Raiteri | Reactive Force Field Simulation of Proton Diffusion in BaZrO 3 Using an Empirical Valence Bond Approach[END_REF] As such, there still remains a need to develop a clear understanding of the trapping effect and the effect of lattice dynamics on mobility at intermediate-to-high temperatures, the temperatures relevant to the operation of the device.

In this work, we combine static DFT calculations and ab initio molecular dynamics to obtain insights into proton diffusion and trapping in fully hydrated Y-doped BaZrO 3 . Static DFT calculations are performed on structures with two different Y-dopant concentrations (12.5% and 25%) and are used to determine the energies of a large number of configurations and correlate energies with distortions of the lattice. The effect of the relative position of yttrium atoms is explored by using a variety of structures in both static and DFT-MD simulations; we assume that no oxygen vacancies are present to simplify the analysis. Ab initio molecular dynamics simulations are then conducted on structures with a 12.5% dopant concentration at three temperatures. From the generated trajectories, the proton diffusion was analyzed through mean square displacements as well as via the numbers of proton jumps between different environments. The results show a relationship between specific local geometries and the occurrence of proton jumps. In particular, intra-octahedral jumps are only seen when short hydrogen bond lengths are formed, more common at high temperatures. The calculated activation energies are in good agreement with previously reported data corresponding to local proton motion. [START_REF] Gomez | The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO 3 , a Cubic Perovskite[END_REF][START_REF] Blanc | Dynamic Nuclear Polarization NMR of Low-γ Nuclei: Structural Insights into Hydrated Yttrium-Doped BaZrO 3[END_REF] for a dopant concentration of 12.5% were combined with results obtained in this work for a dopant concentration of 25%. All geometry optimizations were achieved using the PBE functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF][START_REF] Zhang | Generalized Gradient Approximation Made Simple[END_REF] and were done with a 2 × 2 × 2 supercell having the generic formula BaZr 1-x Y x O 3-x (OH) x . The supercells contain 41 atoms for the 12.5% dopant concentration and 42 atoms for the 25% dopant concentration; i.e., the material is considered fully protonated and there are no oxygen vacancies. Four of the structures with a 12.5% dopant concentration and all 26 structures with a 25% dopant concentration were optimized using CASTEP, [START_REF] Clark | First Principles Methods Using CASTEP[END_REF] a basis set cutoff of 680 eV, an energy convergence criterion of 10 -6 eV (as in ref 14), and a 3 × 3 × 3 Monkhorst-Pack 27 k-point mesh. The 96 other structures with a 12.5% dopant concentration have been reported in a previous work [START_REF] Gomez | The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO 3 , a Cubic Perovskite[END_REF] and were generated, using VASP, [START_REF] Kresse | Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set[END_REF][START_REF] Kresse | From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method[END_REF] as follows. The positions of all atoms in an initial 2 × 2 × 2 supercell with an yttrium dopant, without any protons, are relaxed. To compensate the charge imbalance, VASP automatically uses a uniform positive background gas. A proton is then placed in all possible accommodating sites (24 oxygens × 4 possible orientations), [START_REF] Gomez | The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO 3 , a Cubic Perovskite[END_REF] and its position is relaxed while keeping all the other atoms fixed. For these 96 proton binding structures, started from the lowest energy Glazer distortion, the geometry optimizations were done using the Generalized Gradient Approximation (GGA) PBE functional. A basis set cutoff of 600 eV, an electronic energy convergence criterion of 10 -4 eV, and a 2 × 2 × 2 Monkhorst-Pack k-point mesh were used. The geometry optimization energy convergence criterion was set to 10 -3 eV. In the end, a total of 100 structures with a 12.5% dopant concentration and 26 structures with a 25% dopant concentration were obtained. VESTA software was used for plotting all the structures shown in this article. [START_REF] Momma | VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data[END_REF] 2.2. DFT-MD Calculations. For the DFT-MD study, larger 4 × 4 × 4 supercells containing a total of 328 atoms (including 8 protons) and having a single dopant concentration of 12.5% were simulated. All the calculations were done with the CP2K software [START_REF] Hutter | Atomistic Simulations of Condensed Matter Systems[END_REF] using the PBE functional and run in the NPT ensemble under atmospheric pressure. The cell is cubic with an initial lattice size of 17.04 Å. The Gaussian Plane Wave (GPW) implementation uses a Gaussian DZVP basis set with GTH pseudopotentials, [START_REF] Krack | Pseudopotentials for H to Kr Optimized for Gradient-Corrected Exchange-Correlation Functionals[END_REF][START_REF] Goedecker | Separable Dual-Space Gaussian Pseudopotentials[END_REF][START_REF] Hartwigsen | Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn[END_REF] a 400 eV plane wave cutoff, with only the Gamma point, and a 6 × 10 -7 eV electronic energy convergence criteria. Three different distributions of the yttrium dopant atoms were considered: a uniform distribution (noted YU), a clustered configuration where two pairs of yttrium atoms are nearest neighbors (noted YC), and a planar configuration where two pairs of yttrium atoms are next nearest neighbors (noted YP). This is illustrated in Figure 2, which shows the initial configurations for the DFT-MD simulations. To obtain the starting configurations, eight protons were distributed in the system so that there are two protons in each H-Y configuration as defined in Figure 1 (near, planar, perpendicular, distant). [START_REF] Blanc | Dynamic Nuclear Polarization NMR of Low-γ Nuclei: Structural Insights into Hydrated Yttrium-Doped BaZrO 3[END_REF] Three temperatures (500, 750, and 1000 K) were investigated in order to estimate the activation energies related to proton diffusion. All the DFT-MD simulations were run for 20 ps with a 0.5 fs time step.

RESULTS AND DISCUSSION

3.1. Static DFT Calculations. As discussed in the Introduction, the relative positions of the protons and yttrium atoms, i.e., the H-Y configuration, have a large influence on the energy of the 12.5% dopant concentration cells. To investigate the impact of the relative positions of the yttrium atoms, i.e., the Y-Y configuration, the configurational energies and cell geometries of 26 structures with a dopant concentration of 25% are first analyzed. Figure 3 shows the relative energies of these structures with respect to the lowest energy structure and illustrates the Y-Y configurations.

From Figure 3, it is clear that a larger range of energies is observed for first and second shell than for third shell configurations. We also observe that the lowest energies are obtained for protons nearby Y-Y in first and second shell configurations. We note, however, that the Y-Y arrangements are determined during the high temperature synthesis of the material to form BaZr 1-x Y x O 3-x : hydration is unlikely to change the Y-Y configurations at least at moderate temperatures. Due to the high temperatures used to prepare these materials, yttrium atoms are likely to be randomly distributed in the lattice as suggested by previous researchers. [START_REF] Kreuer | Proton -Conducting Oxides[END_REF][START_REF] Bjorketun | A Kinetic Monte Carlo Study of Proton Diffusion in Disordered Perovskite Structured Lattices Based on First-Principles Calculations[END_REF] Different syntheses routes may lead to different materials properties which could perhaps be related to different dopant spatial distributions. [START_REF] Bi | Synthesis Strategies for Improving the Performance of Doped-BaZrO 3 Materials in Solid Oxide Fuel Cell Applications[END_REF][START_REF] Pergolesi | High Proton Conduction in Grain-Boundary-Free Yttrium-Doped Barium Zirconate Films Grown by Pulsed Laser Deposition[END_REF] To characterize the lattice distortions due to the presence of an yttrium atom, the metal-oxygen (Zr-O, Y-O) and metalmetal (Zr-B, Y-B) interatomic distances for different dopant concentrations (B = Zr, Y) are examined. Figure 4a andb shows a comparison of the average interatomic distances calculated from the static DFT calculations and the experimental values. [START_REF] Goedecker | Separable Dual-Space Gaussian Pseudopotentials[END_REF] The metal-oxygen and metal-metal distances obtained from DFT are in good agreement with experiments although slightly overestimated (by around 1-5%). These results are consistent with other simulation results published on similar materials [START_REF] Giannici | Long-Range and Short-Range Structure of Proton-Conducting Y:BaZrO 3[END_REF] and show that the PBE functional is suitable for such a study. From Figure 4a andb, it is clear that while the average metal-oxygen distance is larger for yttrium-oxygen than for zirconium-oxygen, the average metal-metal interatomic distances are relatively unaffected by the presence of yttrium. This suggests that there is a bending of the B-O-B units, consistent with the octahedral tilting reported previously. [START_REF] Giannici | Long-Range and Short-Range Structure of Proton-Conducting Y:BaZrO 3[END_REF][START_REF] Gomez | The Effect of Octahedral Tilting on Proton Binding Sites and Transition States in Pseudo-Cubic Perovskite Oxides[END_REF] Figure 4c shows the interatomic distances for different environments around the oxygen in the presence or absence of protons. While the presence of a single yttrium atom only leads to a limited variation in the metalmetal distance, the presence of two yttrium atoms results in a noticeable decrease in distance. We also note that the metalmetal distances for atoms surrounding an OH unit are systematically larger than in the absence of proton. This could be due to steric or electrostatic repulsions between the proton and the adjacent Zr/Y atoms.

To investigate the effect of local distortions close to the H and Y atoms on the energy of the cells, we focus on the 100 structures with a 12.5% dopant concentration that contain only one Y and proton per supercell. (In the 25% doped structures, there are two protons and two yttrium atoms, and it is thus difficult to separate the effects of the two H-Y configurations; furthermore, the relative energies of the different Y-Y configurations would need to be accounted for.) To characterize the local distortions shown in Figure 5a, we used a signed "B-OH-B bending distance" which is the nonzero component of the vector from the center of the solid black line connecting the two transition metals (Zr or Y) to the oxygen covalently bonded to the proton. The 100 structures are classified into two categories depending on the sign of the B-OH-B bending distance. If the B-OH-B bending distance is negative, the structure is assigned as "bent outward" since the hydroxyl is outside of the interior BOB angle; the structure is otherwise assigned as "bent inward" when the hydroxyl is on the inside of the interior BOB angle.

Figure 5b shows a plot of the relative cell energies as a function of the H-Y configuration. The different structures are colored according to their "bent inward" or "bent outward" bending status. It is very clear from this figure that, for a given H-Y configuration, all the structures which are "bent inward" have lower energies than the "bent outward" configurations and are thus more favorable. To confirm this observation, we plot the relative cell energies as a function of the B-OH-B bending distance. This is shown in Figure 5c where the different structures are now colored according to their H-Y configuration. While the trend is not perfect and while the energies of some structures overlap, there is a clear correlation between the cell energy and the type (inward/outward) and amplitude of the bending. Independent of the bending type, larger absolute distances (i.e., larger bendings) are correlated with lower energies. We note that all the points corresponding to specific H-Y configurations are grouped and that the lowest energies are observed for combined near/planar and bent inward configurations. As far as we know, this correlation between specific distortions and cell energies has not been highlighted before and confirms that the presence of yttrium dopant close to the protons can lead to low energy configurations not accessible otherwise.

We now investigate what features makes the lowest energy structures more stable. While the proximity of protons, oxygen, and yttrium atoms can lead to complex hydrogen bonding situations, which can influence both the local structure and proton dynamics, [START_REF] Gomez | The Effect of Yttrium Dopant on the Proton Conduction Pathways of BaZrO 3 , a Cubic Perovskite[END_REF][START_REF] Sahraoui | Structural Properties of Y-Doped BaZrO 3 as a Function of Dopant Concentration and Position: A Density Functional Study[END_REF][START_REF] Davies | Dopant and Proton Incorporation in Perovskite-Type Zirconates[END_REF] here we focus on the characterization of these hydrogen bonds to determine how important this is in determining the calculated configurational energies. To describe the different environments, a number of angles and distances, all defined in Figure 6a, are used. dHO 1 and dHO 2 are the distances between the proton and either the second or third nearest oxygen, respectively. dHO 1 is thus the traditional hydrogen bond length, while dHO is the length of the covalent bond (distance to the nearest oxygen). OHO 1 corresponds to the hydrogen bond angle, while OBO 1 and OBO 2 characterize the distortions around the metal atoms. Figure 6b-d clearly shows that the lowest energy configurations correspond to specific local environments where the hydrogen bonds are the strongest. Indeed, Figure 6b shows that all the lowest energy structures correspond to hydrogen bond lengths smaller than https://dx.doi.org/10.1021/acs.jpcc.0c04594 2.15 Å and angles larger than 130°. Figure 6c illustrates the fact that all the lowest energy structures correspond to situations where one of the hydrogen bond (HO 1 ) is noticeably stronger than the other one (HO 2 ). In these low energy environments, dHO 1 is clearly smaller than dHO 2 . Interestingly, these stronger hydrogen bonds are correlated with less distorted OBO 1 angles, i.e., with angles close to a value of 90°which would be seen in an ideal cubic structure (see Figure 6d). It is worth noting that, for the remaining configurations, a number of data points are close to the "y = x" line. In these cases, the proton is engaged in two very similar hydrogen bonds. In the case of inward bending, the formation of hydrogen bonds is intuitively favorable as it is consistent with the octahedral tilting while in the case of outward bending, the formation of hydrogen bonds goes against the expected tilting. This phenomenon is illustrated in Figure S1.

Overall, the static DFT calculations indicate that, in agreement with previous studies, the presence of yttrium as a dopant can lead to the existence of low energy configurations which could act as proton traps. These low energy structures are apparently stabilized by a strong hydrogen bond in a locally distorted environment. We now turn to DFT-MD simulations to investigate the existence and fate of such configurations in a dynamic system.

3.2. DFT-MD Calculations. In order to obtain insight into the proton dynamics in Y-doped BaZrO 3 materials, a number of DFT-MD simulations were performed on 4 × 4 × 4 supercells having a single dopant concentration of 12.5% at three temperatures (500, 750, and 1000 K). In the smaller 2 × 2 × 2 supercells used above for the DFT analysis of the static systems, there was only one Y-Y arrangement, with the different configurations arising from the different positions of the proton in the cell (Figure 1); however, the larger 4 × 4 × 4 cells require that we now consider the relative Y positions (as was performed for the smaller cells containing 25% Y, Figure 3). Three distributions of yttrium atoms were investigated (Figure 2). In one of the systems, the Y atoms are distributed uniformly and quite far apart from each other (YU); while in the clustered system (YC), some of the Y atoms are nearest neighbors; and in the last configuration, some of the Y atoms are next nearest neighbors, i.e., in the planar configuration (YP). The initial structures are generated so that the eight protons are distributed equally in all H-Y configurations (two protons in near, two in planar, two in perpendicular, and two in https://dx.doi.org/10.1021/acs.jpcc.0c04594 distant). The molecular dynamics simulations are then conducted, and the protons trajectories are recorded.

Figure 7a provides an illustration of the environments visited by one of the protons along a trajectory for the YU system at 1000 K. While the proton is originally in a planar configuration at t = 0, it jumps between planar and near configurations, ultimately visiting all possible H-Y configurations along the 20 ps trajectory. The actual trajectories of the eight protons present in the same YU system at 1000 K are shown in Figure 7b. The total residence time fractions of the protons in the different H-Y configurations, i.e., the total time protons spend in a given configuration divided by the total simulation time, for all systems are shown in Figure S2 (these time fractions are normalized by the available sites in each system in Figure S3). They confirm that all H-Y configurations are explored by protons. While the perpendicular and distant configurations are generally slightly less favorable than the near and planar ones (Figure S2), the differences between residences times between different sites are not as large as might be expected based on the static DFT calculations (Figure 1). However, direct comparison with predictions from static DFT energies is out of reach, since these total residence time fractions are not fully converged for the 20 ps trajectories conducted here (Figure S4). Note that the initial configurations of the protons were generated to explore as many microstates as possible, and since all of the configurations are sampled on the time scales used here, the short times employed in the MD runs should still allow comparisons between structural properties obtained in the static calculations and by DFT-MD.

Since it is possible to distinguish the different proton environments with respect to the yttrium atoms (near, planar, perpendicular, distant), we can calculate various properties according to each H-Y configuration. Figure 8 shows the average residence time for each H-Y configuration, i.e., the average time that a proton stays in a given site before jumping to another site. For most of the configurations, the average residence time decreases as expected when the temperature increases. For the perpendicular configuration, however, the average residence time stays short even at 500 K. In fact, the average residence time in the perpendicular configuration is the shortest for all the systems considered. In the static DFT calculations, the configurational energies for near and planar are low and close to each other, while the perpendicular and distant configurations have similar and high energies. Following these results, longer residence times could be expected for near and planar configurations, but this is not what is observed here. A low energy barrier for proton transfers between near and planar might explain the relatively small residence times determined for these H-Y configurations. In addition, fluctuations of the oxygen-oxygen distances are probably an important component in defining the time between jumps, as discussed below.

Proton transfers and rotations between environments are examined by counting how many proton jumps occur from a given type of site to another. The total number of rotations and intra-octahedral transfers are given in Table S1. It is worth noting that no interoctahedral jumps are observed in any of the simulated systems. This is consistent with an activation barrier being higher for this event compared to rotations and intra- octahedral jumps. [START_REF] Merinov | Proton Diffusion Pathways and Rates in Y-Doped BaZrO 3 Solid Oxide Electrolyte from Quantum Mechanics[END_REF][START_REF] Gomez | The Effect of Octahedral Tilting on Proton Binding Sites and Transition States in Pseudo-Cubic Perovskite Oxides[END_REF][START_REF] Munch | Proton Diffusion in Perovskites: Comparison between BaCeO 3 , BaZrO 3 , SrTiO 3 , and CaTiO 3 Using Quantum Molecular Dynamics[END_REF] There are approximately 20 times more rotations than intra-octahedral transfers, which is consistent with activation barriers reported previously. [START_REF] Merinov | Proton Diffusion Pathways and Rates in Y-Doped BaZrO 3 Solid Oxide Electrolyte from Quantum Mechanics[END_REF] When two yttrium atoms are nearest neighbors (YC), there are many more near-to-near proton jumps. This could be because the local distortions resulting from the proximity of To relate these jump events to local distortions, the jumps are now divided into jumps from inward and outward configurations. The number of rotations classified in this way are given in Table S2 and plotted in Figure 9 along with the fraction of inward and outward states (note that while the trajectories are short, these fractions are quite well converged as can be seen in Figure S5). This analysis shows that all intraoctahedral transfers occur when the protons are in the outward configuration. Figure 9 also shows that, in most cases, even though on average there are more protons in the outward states, more rotations originate from the bent inward configuration. This indicates that inward configurations are more favorable for rotations. These considerations, combined with the lowest energies observed for near and planar H-Y configurations with inward bending indicate that the presence of yttrium atoms favor inward bending, and thus rotations, thereby reducing the occurrence of intra-octahedral jumps and thus long-range diffusion. This could, at least partly, explain the proton trapping observed in Y-doped hydrated BaZrO 3 . It is worth nothing that while there are fewer protons in inward configurations than in outward configurations at the temperatures simulated here, this proportion decreases with increasing temperature, in agreement with a lower energy for inward configurations. A rough estimate of the energy difference between inward and outward from their populations at 500 and 1000 K indicates that the inward configurations energies are around 8-10 kJ mol -1 lower than the outward (corresponding to 0.10 eV cell -1 in units comparable to the static DFT calculations reported in this work) which is consistent with the range of values observed in Figure 5c.

We now try to correlate these results with a description of the hydrogen bonding in the material. To begin with, the same geometric descriptors as in Figure 6 are used to see if major differences between inward and outward configurations can be identified. The plots obtained are shown in Figure S6. As expected, due to thermal motion, the clouds of points span a much wider area than that in the case of the static DFT calculations. Moreover, the hydrogen bond lengths and the OBO 1 and OBO 2 angles are on average a bit larger for the inward configurations. To obtain further insights into the difference between inward and outward configurations, histograms of several descriptors of the hydrogen bonding environment are plotted; the histograms were all checked to ensure that they are well converged over the simulated trajectory for all systems. The histograms for the YU system simulated at 750 K are given in Figure 10. Similar results are observed for the YC and YP systems at 750 K (Figures S7 andS8). Figure 10 shows that the covalent bond lengths are similar for inward and outward configurations and as such cannot be used to interpret the proton jumps results. In contrast, inward and outward configurations show some difference with regard to the hydrogen bond length and HO-O 1 distance. On average the outward configurations show shorter hydrogen bond lengths, in contrast to what is observed in the static calculations, and shorter HO-O 1 distances which would tend c,d) Histograms of hydrogen bond lengths and HO-O 1 distances for the YU system at 1000 K. Black lines correspond to all the hydrogen bonds existing in the simulation box over time. Red and blue lines correspond to lengths at the time of a rotation or intra-octahedral jump, respectively. The number of distances in each category is very different (320,000 for the entire set of lengths, up to a few thousands for rotations and up to a few tens for intra-octahedral jumps), so, to facilitate the comparison, histograms are given as a percentage in a given category.

to reduce the activation barriers for the jump to the H-bonded O 1 site and thus increase the jump rate. As such, intraoctahedral transfers are possible in the case of protons in outward configurations, but only rotations, in which the covalent OH bond is not broken, are observed for protons in the inward configurations. It is worth highlighting here that local distortions, motion of oxygen that is covalently bonded to the proton, and rotation of protons can lead to inward/ outward exchange and so the conclusions above do not imply that protons which are in inward configurations are blocked for the entire trajectory.

As proton jumps are highly correlated with the nature of the hydrogen bonds in the structure, we then examined the effect of temperature on the hydrogen bond length more closely. Figure 11a shows the histograms for all the hydrogen bond lengths existing in the YU system at the three temperatures studied here. Similar plots for YC and YP are given in Figure S9 and show once again that the Y-Y distribution does not affect significantly the results (as confirmed by the direct comparison of histograms in Figure 11b). At higher temperatures, wider distributions of lengths are found, corresponding to larger fluctuations in the hydrogen bond lengths. Both shorter and larger lengths are explored, with the former favoring intra-octahedral jumps, which occur preferentially when the hydrogen bonds, and HO-O 1 distances, are the shortest. To confirm this last statement, histograms of hydrogen bond lengths and HO-O 1 distances when jumps occur were plotted. Such histograms are shown in Figure 11c and d for the YU system and in Figure S9 for YC and YP (other temperatures give similar results). Here, histograms are normalized by the total number of distances found in each category (total, rotations, intra-octahedral jumps) to ease the comparison. From these histograms, it is very clear that intraoctahedral jumps only occur for short hydrogen bonds, a situation that happens more often when the temperature increases. For rotations, the maximum of the hydrogen bond length distribution is only slightly shifted to larger distances compared to the average bond length. The results clearly indicate that low temperature static calculations do not necessarily capture relevant motional processes at higher temperatures. It is worth pointing out that while the possible existence of oxygen vacancies is omitted in this study, such vacancies could be present, especially at high temperatures, and introduce compositional dependent chemical expansion effects which would potentially impact the proton dynamics. [START_REF] Bjørheim | On the relationship between chemical expansion and hydration thermodynamics of proton conducting perovskites[END_REF][START_REF] Kjølseth | Determination of the enthalpy of hydration of oxygen vacancies in Y-doped BaZrO 3 and BaCeO 3 by TG-DSC[END_REF] We now examine the proton diffusion in a more global fashion to see if the local dynamics observed lead to notable difference on long-range diffusion. From the DFT-MD trajectories, the proton diffusion coefficients are determined using the Einstein relation:

r r Dt MSD ( ) 6 i t i 0 2 = ⟨ - ⟩ =
where the MSD is the mean square displacement, r i t is the position of proton i at time t, D is the diffusion coefficient, and ⟨ ⟩ denotes an average over all protons. The MSDs calculated for the different trajectories are given in Figure S10, and the diffusion coefficients obtained from the MSDs and the corresponding activation energies, calculated using Arrhenius law and the values at different temperatures, are gathered in Table 1 along with values from the literature obtained experimentally or from classical molecular dynamics. Our calculations seem to overestimate the diffusion coefficients compared to previously reported values at similar temper-atures. It is important to note that, due to the computational limitations of DFT-MD, these simulations are done on a small number of protons and the time scales and system sizes may be too short to assess diffusion coefficients, which could explain this overestimation, as discussed in more detail below.

The activation energies calculated (20.4, 16.0, and 21.3 kJ mol -1 ) are of the same order of magnitude for the three Y distributions studied. This is consistent with the relatively similar properties observed for these systems so far and with the fact that all structures contain a similar distribution of H-Y configurations (Figure S3) so that the global motion is always an average over protons in different configurations. The values calculated here are in the same range as the activation energy of 24.1 kJ mol -1 reported in the computational work of Gomez et al. in which no yttrium is present and thus no trapping can occur. [START_REF] Gomez | The Effect of Octahedral Tilting on Proton Binding Sites and Transition States in Pseudo-Cubic Perovskite Oxides[END_REF] Kinetic Monte Carlo (KMC) simulations, which simulated proton dynamics using 96 binding sites with a 12.5% dopant concentration, have predicted activation energies of 38 kJ mol -1 for a single proton [START_REF] Krueger | Centrality Measures Highlight Proton Traps and Access Points to Proton Highways in Kinetic Monte Carlo Trajectories[END_REF] and 43 kJ mol -1 when eight protons were included in the simulation. [START_REF] Gomez | Effects on the Proton Conduction Limiting Barriers and Trajectories in BaZr 0.875 Y 0.125 O 3 Due to the Presence of Other Protons[END_REF] In these simulations, the protons spend a significant amount of time in Y-OH-Z sites but travel between these sites through perpendicular and planar sites. These activation energies are closer to the values of 40.5 and 45.0 kJ mol -1 determined experimentally. [START_REF] Kreuer | Proton -Conducting Oxides[END_REF][START_REF] Yamazaki | Proton Trapping in Yttrium-Doped Barium Zirconate[END_REF] KMC simulations are more suitable to characterize long-range diffusion because of their ability to sample rare events more efficiently. It is worth noting though that, in both KMC and DFT-MD simulations, the distant sites were sampled only occasionally. In addition, the KMC and the classical molecular simulations considered in Table 1 [START_REF] Raiteri | Reactive Force Field Simulation of Proton Diffusion in BaZrO 3 Using an Empirical Valence Bond Approach[END_REF][START_REF] Van Duin | ReaxFF Reactive Force Field for the Y-Doped BaZrO 3 Proton Conductor with Applications to Diffusion Rates for Multigranular Systems[END_REF] simulate proton motion on a potential energy surface derived from static DFT calculations, typically through geometry optimizations and Nudge Elastic Band calculations. This could also explain partly the discrepancies described here.

Figure 7 shows that, in the DFT-MD simulations, in some cases, a proton can jump between perpendicular and planar sites at the beginning of a trajectory before spending a longer time in Y-OH-Z sites with occasional excursions through planar sites. While not all protons show a similar trajectory and there are also occurrences of protons escaping from near or planar to distant sites, such a behavior, consistent with KMC simulations, has been observed for several protons in the DFT- MD trajectories. This is again a hint that the current DFT-MD trajectories may be too short to look at long-range diffusion but that trapping could be observed more clearly in longer simulations, unfortunately hard to achieve due to computational cost.

Looking more closely at the experimental data, the calculated activation energies are intermediate between the activation energies determined at high (16.0 kJ mol -1 ) and low temperatures (45.0 kJ mol -1 ) by Yamazaki et al. [START_REF] Yamazaki | Proton Trapping in Yttrium-Doped Barium Zirconate[END_REF] using impedance spectroscopy and H/D exchange experiments. The difference between the high and low temperature impedance data of 29 kJ mol -1 was then ascribed to an association energy the protons need to overcome in order to experience longrange diffusion. In the same work, the authors also characterized proton motion by T 1 relaxation measurements using NMR experiments and reported a single activation energy of 13 kJ mol -1 across the whole temperature range monitored (room temperature to 873 K). The T 1 measurements are a probe of local motion while the impedance spectroscopy gives information on long-range diffusion. This picture is consistent with the fact that, due to short time and length scales, our simulations probe mainly the local motion and only rare events of "trapped" to "trap-free" positions are sampled, which are responsible for the long-range transport. As a consequence, while some trapping is captured in the simulations, the rare events corresponding to protons going from "trapped" to "trap-free" states are not sampled with sufficient frequency to affect the activation energy in a significant way.

CONCLUSION

In this work, a combination of static DFT calculations and DFT-MD was used to study the structural properties and proton dynamics in Y-doped BaZrO 3 with two dopant concentrations . The presence of yttrium leads to significant lattice distortions, which, as quantified by, for example, the calculated metal-metal and metal-oxygen interatomic distances, are in agreement with reported experimental results. Going beyond a simple description of the average interatomic distances, the local distortions of the lattice have been characterized. This allowed us to identify two types of configurations related to the bending of the metal-oxygenmetal angle. The nature of the two configurations, labeled inward bending and outward bending, were shown to affect the energies of the systems (Figure 12). In particular, (i) the larger the bending, the lower the energy and (ii) the lowest energies were observed for structures combining near and planar configurations with an inward bending. The lowest energy structures were the ones with the strongest hydrogen bonds, which is only possible in the inward configuration in the static calculations.

DFT-MD studies allowed us to explore the effect of temperature on the local structures and proton trajectories. A detailed analysis of the proton jumps according to the different H-Y configurations showed that inward configurations favor rotations while intra-octahedral jumps are only observed from outward configurations. In fact, in contrast to static DFT calculations, the DFT-MD simulations indicate that outward bending leads to shorter hydrogen bond lengths and HO-O 1 distances on average. This is consistent with favoring intra-octahedral jumps, as a decrease in the HO-O 1 distance should decrease the energy barrier for the jump. The DFT-MD simulations confirm very clearly that intra-octahedral jumps only occur when the hydrogen bond lengths and H-O 1 distances are the shortest. It is important to note that local distortions and rotations of protons can lead to inward/ outward exchange. As a consequence, the statements above do not imply that protons in inward configurations at a given time are blocked for the entire trajectory. The proton diffusion coefficients determined at three different temperatures allowed for the extraction of activation energies, giving values in the same range as some experimental results that probe shortrange motion but in general lower than previously reported simulation results using KMC or classical molecular simulations based on potential energy surfaces determined using static DFT methods. While trapping is believed to happen, it only affects the long-range diffusion and the local motion of protons is largely unaffected. As a consequence, the simulations reported in this work, probably too short to study long-range diffusion, do not capture the larger activation energies for longrange proton diffusion. Overall, our results demonstrate the importance of the bending for proton dynamics, and as local distortions are related with the yttrium content, this work is a new step toward understanding the existence of an optimal dopant concentration for the proton conduction. Finally, the results highlight the role of temperature and lattice motion on controlling the mechanisms for proton motion.

Illustration of the octahedral tilting, total residence time fractions, number of H-Y configuration sites in the systems, numbers of proton jumps, fractions of inward and outward bending as a function of time, characterizations of the interatomic distances related to the hydrogen bonding, mean square displacements (PDF)

Figure 1 .

 1 Figure 1. (a) Four possible proton-yttrium (H-Y) configurations in 2 × 2 × 2 supercells with a 12.5% dopant concentration as defined by Blanc et al. 14 Barium atoms in the A sites of the perovskite structure (at the center of the Zr/Y-O cubes) are not shown for clarity. (b) Relative energies of these cells with respect to the lowest energy structure after geometry optimization, extracted from the data of Blanc et al.[START_REF] Blanc | Dynamic Nuclear Polarization NMR of Low-γ Nuclei: Structural Insights into Hydrated Yttrium-Doped BaZrO 3[END_REF]

Figure 2 .

 2 Figure 2. Starting configurations used for the DFT-MD simulations performed on cells containing 12.5% yttrium. (a) YU: Yttrium dopant atoms are distributed uniformly in the 4 × 4 × 4 supercell. (b) YC: Some of the yttrium dopant atoms are nearest neighbors (clustered). (c) YP: Some of the yttrium dopant atoms are in the planar configuration. In all cases, there are eight protons with two of them arranged in each H-Y configuration (near, planar, perpendicular, distant).

Figure 3 .

 3 Figure 3. (a) Relative configurational energies with respect to the lowest energy structure for 26 cells having a dopant concentration of 25% and different Y-Y configurations. (b) Schematic representation of the different Y-Y coordination shells (these structures are not optimized and simply illustrate the relative positions of the yttrium atoms).

Figure 4 .

 4 Figure 4. (a, b) Metal-oxygen and metal-metal interatomic distances as a function of the Y-dopant concentration calculated in this work from static DFT (12.5 and 25%) and determined experimentally from X-ray diffraction and extended X-ray fine structure (EXAFS) analysis (0, 6, and 15%) by Giannici et al. 37 (c) Metal-metal interatomic distances for specific metal-oxygen-metal configurations obtained in this work.

Figure 5 .

 5 Figure 5. (a) Illustration of inward and outward bending, showing how the bending distance is quantified, a B-OH-B bending distance of 0 corresponding to colinear bonds (i.e., no bend). (b) Relative energies of the 100 structures with a 12.5% dopant concentration classified following their H-Y configuration and B-OH-B bending type. (c) Relative energies as a function of the B-OH-B bending distance.

Figure 6 .

 6 Figure 6. (a) Schematic representation of the local geometry around a proton and labels adopted to assess hydrogen bonding related geometries. (b) Hydrogen bond angle as a function of hydrogen bond length (i.e., dHO 1 ). (c) Distance between the proton and the 3rd nearest oxygen (i.e., dHO 2 ) as a function of the hydrogen bond length. (d) Metal-oxygen-metal angles around the hydrogen bond. For (b)-(d) the structures are colored according to their relative energy with respect to the lowest energy structure in black (energies are in eV cell -1 ).

Figure 7 .

 7 Figure 7. (a) Illustration of the various environments visited by a single proton along a 20 ps trajectory for the YU system simulated at 1000 K. The near configuration is split into two environments depending on the proton having one or two yttrium atoms as nearest neighbors ("Near Y-OH-Z" and "Near Y-OH-Y"). Note that, in YU and YP, there are no Y-OH-Y configurations so this environment cannot be visited in the trajectory shown here. (b) DFT-MD trajectories of the eight protons present, visualized using the OVITO program.[START_REF] Stukowski | Visualization and Analysis of Atomistic Simulation Data with OVITO-the Open Visualization Tool[END_REF] The proton whose trajectory is analyzed in (a) is indicated.

Figure 8 .

 8 Figure 8. Average residence time of a proton in the various H-Y configurations.

Figure 9 .

 9 Figure 9. (a) Average fraction of inward bending states observed in the different systems. (b) Total number of rotations originating from inward/ outward configurations.

Figure 10 .

 10 Figure 10. Several distances for different proton positions analyzed according to inward/outward bending states in the YU system at 750 K: (a) distance between the proton and its closest oxygen atom, (b) hydrogen bond length, (c) distance between two oxygen atoms forming a covalent and a hydrogen bond to the same proton, and (d) bending of B-O 1 -B angles.

Figure 11 .

 11 Figure 11. (a) Effect of temperature on the hydrogen bond lengths for the YU system. (b) Comparison of the hydrogen bond lengths for all the simulated systems at 1000 K. (c, d) Histograms of hydrogen bond lengths and HO-O 1 distances for the YU system at 1000 K. Black lines correspond to all the hydrogen bonds existing in the simulation box over time.Red and blue lines correspond to lengths at the time of a rotation or intra-octahedral jump, respectively. The number of distances in each category is very different (320,000 for the entire set of lengths, up to a few thousands for rotations and up to a few tens for intra-octahedral jumps), so, to facilitate the comparison, histograms are given as a percentage in a given category.

Figure 12 .

 12 Figure 12. Scheme summarizing the main conclusions of this work regarding proton transfers. Static DFT calculations have permitted the identification of two configurations, namely, inward and outward bending, which favor different proton transfers. The inward bending is usually lower in energy, especially close to dopant atoms, and favors rotations with no intra-octahedral jumps observed in the DFT-MD simulations conducted here. The outward bending on the contrary allows for intra-octahedral proton transfers although still less frequent than rotations. These observations are consistent with proton trapping close to the yttrium dopant.

  14 https://dx.doi.org/10.1021/acs.jpcc.0c04594 2. COMPUTATIONAL METHODS 2.1. Static DFT Calculations. For the static DFT study, results reported in the literature

Table 1 .

 1 Diffusion Coefficients (/10 -10 m 2 s -1 ) and Activation Energies Calculated for the Nine Systems Studied in This Work Having a Dopant Concentration of 12.5% and from the Literature Corresponding to

	Experimental (Dopant Concentrations of 10% and 20%) 1,8
	and Simulated Values (Dopant Concentrations of 0.46% and 12.5%) 16,17
	system	diffusion @500 K	diffusion @750 K	diffusion @1000 K	activation energies (kJ mol -1 )
	YU	2.60	6.51	35.96	20.4
	YC	3.79	24.10	22.67	16.0
	YP	2.19	16.26	26.29	21.3
	ref 1	0.11	3.19	16.81	41.3
	ref 8	0.07	1.16	3.30	16 (> approximately 1000 K)
					45 (< 1000 K)
					13.0 (T 1 measurements)
	ref 16	0.05	1.31	6.63	40.5
	ref 17	0.01	0.26	1.49	43.4

https://dx.doi.org/10.1021/acs.jpcc.0c04594

https://dx.doi.org/10.1021/acs.jpcc.0c04594

■ ACKNOWLEDGMENTS We thank Dr. Bartomeu Monserrat, Dr. Rachel Kerber, Dr. Fredeŕic Blanc and Profs Yoshihiro Yamazaki and Shu Yamazaki for many fruitful discussions concerning theoretical methods and proton dynamics. C.M. acknowledges an Oppenheimer Research Fellowship from the School of Physical Sciences from the University of Cambridge. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreements No. 714581 and No. 835073). Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/ L000202), this work used the ARCHER UK National Supercomputing Service. M.A.G.'s work was supported by the National Science Foundation under Grant DMR-1709975 and the Mount Holyoke College Department of Chemistry. Computational resources were provided in part by the MERCURY consortium under NSF Grant CHE-1626238.