Cooperative origin of proton pair diffusivity in yttrium substituted barium zirconate

Proton conduction is an important property for fuel cell electrolytes. The search for molecular details on proton transport is an ongoing quest. Here, we show that in hydrated yttrium doped barium zirconate using X-ray and neutron diffraction that protons tend to localize near the dopant yttrium as a conjugated superstructure. The proton jump time measured using quasi-elastic neutron scattering follows the Holstein-Samgin polaron model, revealing that proton hopping is weakly coupled to the high-frequency O-H stretching motion, but strongly coupled to low-frequency lattice phonons. The ratio of the proton polaron effective mass, m*, and the proton mass is m*/m = 2, when coupled to the Zr-O stretching mode, giving experimental evidence of proton pairing in perovskites, as a result of proton-phonon coupling. Possible pathways of a proton pair are provided through Nudge Elastic Band calculations. The pairing of protons, when jumping, is discussed in context of a cooperative protonic charge transport process. Proton conduction through ceramic electrolyte membranes involves the ionic transport of H+ ions through a crystalline lattice. Here, the authors investigate the underlying molecular transport mechanism and demonstrate that proton conduction in yttrium-doped barium zirconate occurs through cooperative proton motion.

P roton diffusivity and conduction is a well-established phenomenon observed in rare earth transition metal oxides with the ABO 3 perovskite structure. The conduction is an ionic transport process of technological interest for protonconducting electrolyte membranes used in intermediatetemperature electrochemical energy converters, such as solid oxide fuel cells and electrolyzers [1][2][3][4][5][6][7] . The traditional, accepted mechanism for proton incorporation in perovskites requires doping with electron-acceptor (p-dopant) ions 8,9 . When an ABO 3 perovskite is doped at the tetravalent B sites with trivalent ions, one oxygen vacancy results from every two doped sites. This vacancy dissociates water from the ambient water vapor into a hydroxyl ion, which fills the vacancy. The second proton binds to another nearby oxygen in the perovskite lattice 10,11 . Upon thermal activation, the O-H bonds dissociate 12 . The freed protons conduct with polaron-type temperature characteristics 13 through the oxygen ion network of the perovskite. One extensively studied perovskite proton conductor system is yttrium-doped barium zirconate (BZY) 3 .
Quasi-elastic neutron scattering (QENS) reveals that the protons assume a rotational mode around the oxygen with a very low, barely detectable, activation energy 14,15 . This qualitatively agrees with computational results, where binding sites are separated by a very small rotational barrier (<0.02 eV) within a plane and a 0.14 eV barrier to rotate between planes as shown in Fig. 2 of ref. 16 . Simulations indicate that there are other, higher energy rotations in the system and more in systems with dopant 17,18 . At higher temperatures, the activation energy for these higher energy rotations as well as intra-and inter-octahedral barriers are overcome and the O-H dissociate and the protons may jump to an adjacent oxygen. With sufficient available thermal energy, net proton transport takes place.
Investigations on tuning the ion mobility by lattice dynamics have received emerging attention in recent years [19][20][21][22] . For proton conductors, experiments and calculations show that proton motion is coupled to local, transitory distortions in the perovskite crystal lattice. The experimental observations show proton conductivity following exponential growth. The temperature dependence reveals that proton conductivity is a thermally activated process with an identifiable activation energy, which is a function of material parameters 23 . Moreover, this is reminiscent of a transport process that has been recently postulated as polaron type 24,25 . In particular, the variation of the proton jump rates with temperature follows, with surprising accuracy, a Holsteintype polaron model that was developed by Samgin 24,26 . Consequently, the proton conductivity can be interpreted as a cooperative process governed by the proton-phonon coupling. It has been shown that the protons are propelled by the lattice vibration 13 . The polaron mass calculated for our case turns out to be almost exactly the mass of two protons 13 . However, the role of each specific phonon mode remains a question. For instance, whether there exist one or several specific phonon modes that govern the proton transport? In our earlier calculations and in calculations in this paper, we always see motion of ions near the proton that aid in the proton motion [16][17][18][27][28][29] . In Samgin's work 24 , this assistance is through low-frequency lattice vibrations. The lattice assistance we see in our earlier computational work [16][17][18][27][28][29] is not directly tied to low-frequency lattice modes in the absence of protons but rather oxygen and other ion motions seen in the transition pathways, i.e., regional motions.
The polaron model was initially created to describe general charge carriers dressed by an elastic field. The development of polaron theory has evolved from electrons in hard condensed matter to charged particles in similar environments. In perovskite-type high-temperature superconductors, charge carriers follow polaronic behavior and condense into pairs 30,31 . The strong electron-phonon coupling plays an important role in pairing of charge carriers 32 . One may therefore ask whether there are further processes in proton diffusion which take place in cooperation. In our experimental and theoretical studies on BZY, we found that pairing of protons can explain some empirical observations. The proton jumping distance and rates can be determined from the quasi-elastic broadening of the neutron peak during thermal activation. We find that the jumping lengths for yttrium-doped barium cerate are approximately 5 Å, which is twice the distance of two adjacent oxygen positions in the yttrium-doped barium cerate lattice (3 Å) 33 . We can illustrate the proton jump scenario by choosing one oxygen lattice position as occupied by one proton, and only the next oxygen site is available for the jumping proton. This inter-octahedra transfer is indicated in Fig. 9 in ref. 33 . A similar pathway is seen in our previous work, for example, in Fig. 4(c) in ref. 34 where in the upper half of the trajectory a proton does a sequence of two intra-octahedral transfers to hop over another proton. The pathway also includes inter-octahedral transfers and rotations. Note that, from thermodynamic arguments, the pairing of protons are energetically favorable for proton conductivity, providing additional support for conclusions based on QENS.
To form conduction pathways, simulations only show movements between binding sites in close proximity 34 . This principle is satisfied if the oxygen vacancies spontaneously arrange to form a crystallographic oxygen vacancy superstructure as has been suggested 35 . The proximity of binding sites is based on the Rietveld refinement of neutron diffraction data. Neutrons are excellent probes for light elements such as oxygen and hydrogen or protons. Also, systemic deficit of oxygen ions (oxygen deficiency) in a crystal lattice in neutron diffractograms can be quantitatively analyzed with neutron diffraction 36 . We found in element-specific anomalous X-ray diffraction studies on BZY Bragg reflections, which we could assign to the yttrium positions in the BZY. This indicated that the yttrium is arranged as a superstructure in BZY at 300 K 37 . Neutron diffractograms on deuterated BZY (after exposure to vapor from heavy water) showed superstructure reflections, which we could assign to the deuteron positions in the crystal lattice. Moreover, these Bragg reflections coincide with the Bragg reflections from the yttrium superstructure. We interpret this systemic coincidence of Bragg reflections as that the deuterons settle-at least at 77 K-preferentially and with statistical significance in correlation with the yttrium. It appears that protons could organize and possibly also move in a landscape set by the yttrium 37 . In addition, nuclear magnetic resonance spectra have shown a low energy structure of two close by protons in Sc-doped BaZrO 3 38 .
Computational studies also show that the dopant identity can have a considerable impact on the proton motion 17,18,28,29,39,40 . Proton motion through BZY tends to have long-range proton conduction pathways right through the dopant regions. This is in contrast to proton motion through aluminum (Al)-doped BaZrO 3 where long-range pathways go through regions without the dopant. In these cases, the dopant still serves as a trap 41,42 and the proton periodically escapes to areas that contain no dopant ion. These regions serve virtually as proton conduction "highways." Kinetic Monte Carlo (kMC) simulation and other mixed electronic structure and graph theory approaches have been used to find proton pathways and conduction possibilities on the landscape set by the dopant 17,18,28,29,34 .
Many computational works have proposed concerted motion 43,44 or other cooperative motion 45,46 of multiple ions during transport. While exploring proton conduction pathways of multiple protons on the landscape set by the dopant, one computation study 34 finds that proton conduction pathways and limiting barriers are affected by the presence of other protons. The lowest energy binding site for a second proton is quite close to the first proton as the first proton distorts the landscape, which facilitates hydrogen bonding. Further, a nearest neighbor peak in the radial distribution function in a kMC simulation at about 4.0 Å shows that conduction often keeps two protons in close proximity. Finally, the long-range limiting conduction barrier rises with proton-proton interaction.
Despite computational simulations showing the cooperative motion between multiple ions during transfer, direct experimental evidence of cooperative ion motion in solids has been rare. In this work, we use a combination of X-ray, neutron diffraction, vibrational spectroscopy and QENS experiments, and electronic structure and nudged elastic band (NEB) simulations to understand what features control the proton conduction in BZY at the molecular level and to explore what consequences occur when multiple protons are involved. We begin with presenting the superstructure of yttrium and protons in BZY, followed by the analysis of proton diffusion mechanism. To find out the specific phonon modes that govern the proton transport, we investigated in detail the influence of various phonon vibration modes in the BZY phonon spectrum through the proton polaron model, from which we show evidence of cooperative proton motion. Then we provide possible moving pathways when two protons are in close interaction. The new two-proton pathways give direct evidence of how the lattice distortions caused by one proton can assist the motion of a second proton.

Results and discussion
Correlation between proton and yttrium superstructure. The crystal structure of BZY determines the landscape for the proton transport. The structure refinement of the neutron diffraction data for BaZr 0.9 Y 0.1 O 2.95 (BZY10) showed a space group of Pm3m with unit cell length of 4.208 Å at room temperature 37 . The neutron diffractogram of deuterated BZY10 reveals additional Bragg reflections at Q~2.78 Å −1 , 3.48 Å −1 , 4.38 Å −1 , and possibly at 2.44 Å −1 (Fig. 2 in ref. 37 , middle panel, shown for the deuterated sample), which match very well with a superstructure k-vector [0.0729824, 0.0729824, 0.0729824]. Absence of these reflections in the dry BZY10 suggests that protons and deuterons occupy specific lattice positions. Naively, the proton concentration should scale with the 5% oxygen vacancy concentration, which means 1 vacancy per 20 unit cells. This translates to a linear dimension of 20 1/3 , or~2.7, unit cell lengths for the periodicity of the proton position, which corresponds to a distance of 11.36 Å. The proton induces a local elastic or electrostatic distortion and propagates throughout the lattice and interacts with the next protons, giving rise to long-range order.
It has been shown that the yttrium participates chemically in the protonation process. It is believed that ionic transport is impeded at higher temperatures when the oxygen vacancies are more ordered. The proximity of the Y 3+ dopant and the oxygen vacancy suggests that oxygen vacancy ordering is a response to Y 3+ ordering. The Y L 1 -edge K-shell X-ray absorption edge of dry and protonated BZY10 37 shows evidence that the yttrium becomes partially reduced when the crystal is protonated.
The neutron diffractograms and Yttrium L 1 -edge X-ray absorption spectra (NEXAFS) suggest that the interaction of yttrium and protons can have a significant effect on order and conduction pathways. Element-specific X-ray diffractograms 37 confirm that, at X-ray absorption threshold of Y at 16 keV, which are particularly sensitive to Y, we see characteristic and systematic extra Bragg reflections ( Fig. 2 in ref. 37 , right panel). The reflection at Q = 4.49 Å −1 for 16 keV, the energy before the Y absorption edge, is split into two peaks at Q = 4.49 and 4.51 The extra reflection is not present at the higher photon energies.
Hence, Y atoms likely form an ordered superstructure. Since the superstructure reflections for protons and for yttrium are not at the same Q positions, the superstructures are not identical. If both superstructures were identical, every Y would have a neighboring proton, which requires saturation with protons. The proton superstructure reflections are at lower Q values than the Y reflection, which indicates that the proton unit mesh is larger than the Y unit mesh and accounts for the under-saturated proton loading in the samples. The proton concentration is 3 mol% as determined by measuring the weight change upon hydration 47 , much less than the Y concentration of 10 mol%.
The Y positions are the anchor points that make the landscape on which the protons move. The Y locations offer additional structure information that may help to solve the underdetermined sets of equations, which describe the proton trajectories. The positions of the Y limit the degrees of freedom for the proton movement.
Phonon-assisted proton transport. When a proton is incorporated in an ionic crystal, it interacts with the lattice vibration. The interaction between the charged proton and its surrounding crystal lattice field can be treated as a quasiparticle: a proton polaron 13,24 . Transport properties in the ionic crystal are dominated by the proton-phonon interactions. The Holstein-Samgin polaron model 24,26 describes the proton jump as due to the interactions with lattice vibration. Here we apply this model to elucidate the correlation between the lattice vibration and proton hopping. The interaction is used to calculate the spectrum of oscillations of the crystal lattice. In this section, we use the QENS data for BZY10 47 to derive parameters related to the proton polaron transport, such as the proton-phonon coupling constants, proton polaron hopping activation energies, and proton polaron effective mass.
In the neutron scattering function for BZY10 47 , the quasielastic broadening becomes wider with increasing temperatures above 600 K, indicating activated proton diffusive motion. The proton jump time is extracted through a Chudley-Elliot (CE) model 48 fit to the QENS half-width at half-maximum (HWHM). Being widely used to study proton diffusion 13,49,50 , the CE model considers the proton transport as successive jumps in random directions from site to site with jump length d and jump time τ. The model QENS HWHM (Γ) is the Q-dependent function: Figure 1a shows the Q-dependent HWHM of QENS scattering function for BZY between 600 and 900 K. The solid lines are the CE model fits to Eq. (1). To obtain reasonable values for the jump time τ, the proton jump length d is set to 1.6 Å for all temperatures. Such a distance is adopted because it is equal to an intra-octahedral transfer, which occurs with the highest probability of 87% among long-range path limiting steps 51 . The thermal expansion coefficient for protonated BZY10 is 4.56 × 10 −6 K −1 37 so that the change of d over the temperature range is <0.01 Å. The obtained proton jump time decreases from 3.55 to 1.79 ps, as the temperature increases from 600 to 900 K, shown as solid symbols in Fig. 1b. In addition, we have also applied the Hall and Ross model for the fit, which considers the proton jump length as a Gaussian distribution 52 . We further considered a wide range of jump lengths, including a rotation with d = 1.4 Å, an intra-octahedral transfer with d = 1.6 Å, an inter-octahedral transfer with d = 1.9 Å, and two successive intra-octahedral transfer with d = 3.2 Å. Computationally, a minimum in between two end points signals that the move considered is composed of two elementary steps. Two intra-octahedral transfers in succession cover a distance of 1.6 × 2 = 3.2 Å. The fit results for both models with d fixed at different lengths are shown in Supplementary Fig. 1. The corresponding jump times τ and jump distances d are listed in Supplementary Table 1. With d fixed at different lengths, the difference between jump times obtained from the two models are within 15%.
In the Holstein-Samgin polaron model 24,26 , the jump time τ is expressed as: where A and u are both constants for the interaction of protons with the lattice related to the specific vibration mode, ω is the phonon frequency that assists the proton transport, k is the Boltzmann constant, and T is the temperature. By inserting a particular excitation phonon frequency into Eq. (2), the jump time becomes a temperature-dependent function. The proton jump times obtained by QENS experiment are in good agreement with the Holstein-Samgin polaron model, as presented in Fig. 1b, where the fitting function is shown as the solid curve. In this way, both constants u and A are determined from the model fit. We will discuss below the situation for different frequency regions.
Interaction of charge carrier with the lattice vibration increases the polaron effective mass. A fundamental task in polaron theory is the determination of the ratio between polaron effective mass m* and the mass m of the charge carrier. For this, we invoke the protocol of electron polarons, where the polaron effective mass is positively correlated to the coupling constant. For weak and intermediate coupling with u ≤ 4.5 53 : and for strong coupling with large u 54 : The polaron hopping activation energy E a can be obtained from the Holstein-Samgin model 24 , which is: By combining QENS, Raman and Fourier transform infrared (FTIR) spectrum, it is possible to map out a full phonon spectrum for hydrated BZY10, as presented in Fig. 2. The Raman and FTIR spectra are in good agreement with literature reports on BZY 55,56 . The phonon bands <110 meV belong to lattice vibration motions. For the generalized density of states extracted from QENS, in the energy range <30 meV, two peaks at 12 and 20 meV becomes prominent as the temperature increases from 500 to 800 K. The phonon band at 12 meV is also visible in the Raman spectrum, corresponding to the acoustic phonon T 2u 57 , a Ba-ZrO 6 stretching motion 58 . Following the density functional theory study on BaZrO 3 Table 1. Supplementary Table 2 compares the coupling constant u and polaron mass m Ã m employing both CE and models and the chosen jump distances between 1.4 and 3.2 Å. The variation in coupling constants and polaron mass are within about 2%. The proton-phonon coupling constant u and effective mass determine the dynamics of proton transport in general. The O-H stretching mode with high phonon energy of 430 meV is weakly coupled to protons with u = 1.8 and a light polaron effective mass ratio to the proton mass of 1.38, suggesting that the proton can be easily knocked out by the O-H stretching modes, which leads to fast translation to a neighboring proton. In simulations, we also observe significant proton displacement caused by such highfrequency modes. The effective mass ratio of 1.6 is found for the O-H bending mode at 200 meV with u = 2.5. On the other hand, the low-energy acoustic phonon at 12 meV yields strong coupling with constant u = 10.2 and a polaron effective mass ratio of around 112 54 , which reflects the known difficulty of proton transport driven by low-energy phonons at low temperatures. It is worth noting that, for the lattice phonon Zr-O stretching at 90 meV with u = 3.7, which fall in the range of intermediate coupling, the polaron effective mass ratio are around 2, with an error <2%. A similar effective mass ratio of 2.08 was obtained for Y-doped BaCeO 3 , when the proton jump is excited by the Ce-O stretching vibration at 355 cm −1 13 , which means: The effective mass m * of a polaron describes the extent of interaction between the charge carrier and the ionic crystal. With electrons as charge carriers, the polaron mass is considerably larger than a factor 2: the ratio m * /m was reported to be 5-10 or even more for high-temperature superconductors 60 . A small polaron mass corresponds to weak interaction with the crystal field, hence a fast transport. Nevertheless, a superprotonic conductivity does not exist in BZY, raising the question whether the underlying phonon modes ν(Zr-O), and ν(B-O) in general, which provides the fortuitously low ratio m Ã polaron =m proton of 2, is the actual driver of proton conductivity in perovskites or maybe the driver of a different state, such as proton pairing, as the factor 2 might suggest. Such state could then be considered a degenerate proton polaron or even a proton bi-polaron. In polaron physics, pairing of charge carriers is energetically favored. When charge carriers pair, the energy of the induced distortion of the lattice is shared and their effective mass is reduced 61 . These results suggest the possibility of proton pairing when jumping originating from proton-lattice interaction, which leads to a reduced polaron mass in perovskite-type proton conductors.
The polaron radius r p determined from activation energies E a and the distances between oxygen ions, d OO , is 24 : where k p is a proportional constant. Equation (7) indicates that the activation barrier is strongly dependent on the lattice parameters, in line with our experimental data measured by impedance spectroscopy when the sample was kept under high compressive strain 62 . Inserting the activation energy 62 in Eq. (7) and the intra-octahedron OO distance at respective pressure 63 yields a polaron radius of 3.34 Å, which is in the order of the lattice constant. This is a small number when compared to other polaron radii and thus may reflect an essential limitation of proton transport in perovskites. This small radius implies that the proton is easily trapped by the local deformation of ions, preventing it from moving in a continuum polarization cloud. In Table 1, the activation energy of around 0.10 eV applies to the temperature range between 600 and 900 K. It is worth to mention that the activation energy for proton diffusion measured by QENS is 0.13 eV at above 700 K 37 . We note that, in the experimental characterization of materials, the activation energy may depend on the length scale that the analytical tool can cover. QENS, a neutron method, probes microscopic or atomic motion such as ion jump from site to site 50 , yielding activation energy of around 0.1 eV, whereas impedance spectroscopy, an electroanalytical technique, probes at the macroscopic scale, and usually gives the activation energy above 0.4 eV 64 . Two-proton configurations and possible pathways. Inspired by earlier studies 17,34 which show that a proton significantly alters the local environment in BaZr 0.875 Y 0.125 O 3 leading to low-energy configurations and trajectories with two protons in very close proximity and the experimental results in this study, we find twoproton configurations and possible two-proton pathways. To be sure that the earlier studies did not have additional effects due to small size of the simulation box, this study considers longer simulation boxes while maintaining the 12.5% doping level. Figure 3 shows the configuration of a 2 × 2 × 4 unit cell system with 12.5% doping without a proton in Fig. 3a and the same structure with a single proton optimized starting from the lowest energy binding site of the smaller system in Fig. 3b. In our earlier study 34 , we added a second proton to the lowest energy first proton configuration at all possible second proton sites. Figure 3 shows the larger system configurations for the lowest and second lowest configurations from our earlier two-proton study 34 in Fig. 3c, d, respectively. Incorporation of protons enhances inter-octahedral OO distance contractions in the horizontal direction. In this work, this inter-octahedral OO distance decreases from 3.8 Å without proton (Fig. 3a) to 3.4 Å (Fig. 3b) with one proton and further to 3.2 (Fig. 3c) and 3.1 Å (Fig. 3d) with two protons. The difference in energy between the lowest energy two-proton structure Fig. 3c and the second lowest energy two-proton structure Fig. 3d is 0.04 eV. Comparison with our earlier work 34 shows that these contractions are effectively the same in the 2 × 2 × 2 and 2 × 2 × 4 unit cell simulation boxes. Since lattice changes are very important to finding lower energy sites of two nearby protons, we did an NEB calculation in the larger simulation box, as shown in Fig. 4, with two protons starting on the same octahedron in Fig. 4a and one of the protons ending in a different octahedron in Fig. 4f. The sequence of steps clearly shows how one proton influences the moves of another nearby proton. Each proton distorts its surrounding lattice to increase hydrogen bonding creating energetically favorable distorted regions for a moving proton to gravitate towards. After the first intra-octahedral transfer to Fig. 4b, the hydroxyl bond rotates to the distorted plane prepared by the second proton Fig. 4c and triggering off subsequent motions. The two protons deform the lattice to Fig. 4d further enhancing hydrogen bonding by decreasing the inter-octahedral distance and lowering the energy. The proton originally on the plane transfers intraoctahedrally leading to two protons in extremely close proximity.
At this point, one of the protons is pushed away to a neighboring oxygen site (Fig. 4e) and rotates out of the plane (Fig. 4f).
In another NEB calculation, the tilting of one hydroxyl within a plane assists another proton's two intra-octahedral transfers as seen in Fig. 5. Since this path contains only two elementary steps, it was divided into two paths and the climbing NEB method was used to find transition states that are shown in between the three minima. The pathway starts at a 0.04 eV minimum in Fig. 5a with both protons on oxygens adjacent to the dopant. The 0.22 eV transition state in Fig. 5b shows a transfer of one proton to a neighboring site as the left hydroxyl tilts downward. Subsequent hydroxyl tilting leads to the 0.14 eV minimum structure (Fig. 5c). The first intra-octahedral transfer has an activation barrier of 0.19 eV. Hydroxyl tilting facilitates a second intra-octahedral transfer with a 0.44 eV transition state at Fig. 5d, which shares the proton between the transferring oxygen ions. After the transfer occurs, there is further hydroxyl ion tilting to get to the 0.19 eV minimum in Fig. 5e. This second move has an activation barrier of 0.30 eV. We observe large displacement of oxygen ions when comparing the images in each step, reflecting the interactions between the backbone lattice and the positively charged protons, and how the proton transport is coupled to the lattice vibration as discussed before.
These are just two possible pathways including all interactions while specifying initial and final configurations. Others could have different sequences of events. While our earlier work on 12.5% Y-doped BaZrO 3 in reference 34 does not include the additional backbone distortions induced by close presence of two protons, it does include many more pathways than this work. For example, the path shown in Fig. 5 where one hydroxyl tilts while the other proton does two intra-octahedral transfers was seen in our earlier work as a part of a probable path from a graph theory search shown in Fig. 4(c) in ref. 34 . Reference 34 uses graph theory to find probable long-range conduction pathways for one moving proton in the presence of a stationary one as well as kMC to find long-range trajectories with all protons moving. This work showed that the average long-range limiting conduction barrier increases from 0.39 to 0.45 eV at 1000 K as the number of protons increases from 1 to 8 in the absence of external field and is most  Fig. 3 The effect of protons on the structure of a 2 × 2 × 4 Y-doped BaZrO 3 unit cell system. The green, blue, red, and white spheres represent barium, yttrium, oxygen, and hydrogen ions, respectively. A variety of distances are highlighted in the lowest energy structure without any protons (a). The lowest energy structure with one proton (b) shows an inter-octahedral OO distance near the proton of 3.4 Å in contrast to 3.8 Å inter-octahedral OO distances further away from the proton. The lowest (c) and second lowest (d) two-proton configurations show inter-octahedral OO distances of 3.1 and 3.2 Å, respectively. The second lowest energy twoproton configuration (d) is only 0.04 eV higher than the lowest energy twoproton configuration (c). Proton-enhanced inter-octahedral OO distance contractions are in agreement with earlier work 34 on a smaller system with very slight differences in the distances shown. commonly an intra-octahedral transfer barrier. This range is in excellent agreement with experiments 64 , which show long-range conduction barrier of 0.43 eV for a system at 10% doping. One result of the kMC simulation not reported in ref. 34 is that the average barrier is only about 0.15 eV, which is significantly smaller than the average limiting conduction barrier and much closer to the local barrier found by the proton polaron model.

Conclusion
In summary, by combining the experimental and computational data, it is possible to depict a generalized scheme for the interaction of protons with the BZY lattice. X-ray diffractograms at the Y K-shell absorption edge suggest that Y atoms form an ordered superstructure in BZY. Y participates chemically in the proton incorporation process, as indicated by the Y L 1 NEXAFS spectroscopy. Neutron diffractograms provide evidence that the protons are organized in a superstructure as well, with larger separation than Y. The density functional theory conjugate gradient minimization calculations reveal that the incorporation of protons causes enhancement in the inter-octahedral OO distance contraction in the crystal lattice. The Holstein-Samgin polaron model fits on proton jump time extracted from QENS data suggest that the proton jump is weakly coupled to high-frequency O-H stretching motion with light proton polaron mass but strongly coupled to the low-frequency lattice vibration yielding heavy polaron mass. The proton polaron effective mass is 2 times of the proton mass, when the proton jump is coupled to the Zr-O stretching vibration at 90 meV, pointing out the possibility of proton pairing during transport. Sample NEB calculations show how one proton can prepare the way for another and how the close proximity can lead to proton transfer. These are only possible paths observed when setting up a specific initial and final configuration but they do illustrate how some cooperative motions could occur. We note that the proton polaron model involves quantum proton treatment to extract parameters related to polaron transport from experiments-a classical model would not have successfully matched the model. However, the computational section only treats electron quantum mechanically and describes proton motion on the Born-Oppenheimer ground state potential. The kMC calculations still treat the proton classically as the temperature of the calculation is 1000 K and therefore the thermal activation is more significant than the tunneling.
Extracting proton conductivity mechanisms in ceramic materials from experiments and simulations has historically used classical rationals, such as diffusion. Only tunneling describes proton transport across high kinetic energy barriers within very short distances at low temperatures 65 . The lack of application of tunneling principles in the proton conductivity framework is surprising, given that there is no fundamental limit that discriminates electrons from protons, except for their large mass difference. There is also no reason to exempt protons from quantum mechanics. It remains to be shown that the polaron density of states can be computed from the phonon density of states. It remains also to be seen whether proton pairing has in addition to its cooperative characteristic also a true synergistic effect, which, if verified, can be optimized by materials tailoring or by optimized device operation. The degree of analogy between electron-and spin-based processes with rest mass zero and proton (and ion in general) based process with finite mass needs to be investigated in future experimental and computational studies.

Methods
Experimental methods. BZY10 was first prepared as powders by solid state synthesis, then compacted into pellets or slabs and sintered at 1700°C in an airvented furnace and subsequently cooled and stored in ambient air. For further analyses, all samples in the form of powder, pellets, or slabs were annealed in dry nitrogen to form the dry sample. Protonation was achieved by annealing a dry sample in a nitrogen gas humidified with water (protonated sample) or heavy water (deuterated sample) 47,66 .
Neutron powder diffractograms with λ = 1.1545 Å were recorded over the range of temperatures from 300 to 773 K in steps of around 100 K, with the HRPT diffractometer 67 at SINQ, the Swiss Neutron Spallation Source. Anomalous X-ray diffraction data on the dry sample for the yttrium K-edge were collected at the Swiss Light Source 68 at 16, 16.9, 17, and 18 keV photon energy from 300 to 773 K in steps of around 100 K. NEXAFS were collected at beamline 9.3.1, Advanced Light Source, Berkeley, CA in the energy range from 2370 to 2400 eV in steps of 0.1 eV and a resolution of ΔE/E~6000.
The FOCUS time-of-flight instrument at the Swiss Spallation Neutron Source was used for QENS. Protonated slabs were arranged azimuthally in a platinum cylinder so that the neutron beam passed through one absorption length of the sample material. The platinum cylinder was equipped with a pressure relief valve to maintain constant pressure and minimize water escape while heated within the QENS beamline. The Q-range of 0.2-0.8 Å −1 was investigated with 6 Å wavelength neutrons, corresponding to an incident energy of 2.273 meV. The instrument resolution was obtained by a separate measurement of a vanadium sample. Reduction of the raw data was performed with the DAVE software package 69 . The QENS data obtained at 300 K served as a reference state for background subtraction, assuming that protons are sufficiently stationary at this temperature. Data from 500 K and the data points at Q = 0.2 Å −1 and 0.3 Å −1 at 600 K was omitted from the data analysis due to poor signal-to-noise ratio. Raman spectrum was measured at room temperature on hydrated BZY10 by using a Senterra R200-L spectrometer with a 514-nm laser excitation line. FTIR was measured at room temperature by Nicolet 6700 Spectrometer. To prepare the specimen for FTIR measurement, a homogeneous mixture of 0.001 g hydrated BZY10 powder and 0.1 g KBr was pressed into a pellet using a load of 7 tons.
Computational methods. Density functional theory using the Vienna ab initio simulation package (VASP) 70 was used. As in our earlier work 17,18 , VASP [71][72][73][74] with the Perdew-Burke-Emzerhof density functional using a generalized gradient approximation was used. The projector augment wave method 75 was used with the valence states of Ba(5s, 5p, 6s, 5d (only 0.01)), Zr(4s, 4p, 5s, 4d), O(2s, 2p), and Y (4s, 4p, 5s, 4d). The labels on these PAW method files within the VASP library are Ba_sv, Zr_sv, Y_sv, Al, O, and H. Further, a 2 × 2 × 2 Monkhorst-Pack k-point mesh was used. The conjugate gradient minimization was used to find minimum energy structures. The NEB method 76 with a force-based conjugate gradient optimizer 77 was used to find two-proton trajectories. As pathways near the Y dopant tend to be lower in energy, we chose to explore trajectories traversing binding sites on oxygen ions adjacent to the Y dopant. The total force was Fig. 5 Illustration of a nudged elastic band pathway where a second proton moves around a first proton via two intra-octrahedral transfers. The green, blue, red, and white spheres represent barium, yttrium, oxygen, and hydrogen ions, respectively. The first intra-octahedral transfer starts at a 0.04 eV minimum structure (a), goes through a 0.22 eV transition state (b), and finishes at a 0.14 eV minimum structure (c). The second intra-octahedral transfer starts at minimum (c) and goes through the 0.44 eV transition state (d) to the 0.19 eV minimum (e). All energies are relative to the lowest energy twoproton configuration found in a 2 × 2 × 4 Y-doped BaZrO 3 unit cell system. converged to <0.1 eV Å −1 . A spring constant of 5.0 eV Å −2 was used between the images. Four images were initially used. When an intermediate low-energy point is found, that image was optimized and new four image NEB pathways were made from the initial minimum to the intermediate minimum and from the intermediate minimum to the final minimum. After convergence, the climbing NEB method 76 was used to find the transition state in each pathway. For cases of multiple minima in between initial and final points, only NEB is done. Since our earlier kMC 34 did not show clustering of more than two protons, only two-proton trajectories were considered.