Nuclear dynamics in BaZr0.7Ce0.2Y0.1O3−δ proton conductor as observed by neutron diffraction and Compton scattering

Concurrent neutron Compton scattering (NCS) and neutron diffraction experiments at temperatures between 70 K and 300 K have been performed on proton-conducting hydrated BaZr0.7Ce0.2Y0.1O3−δ (BZCY72) fabricated by spark plasma sintering. A combined neutron data analysis, augmented with density functional theory modelling of lattice dynamics, has enabled, for the first time, a mass-selective appraisal of the combined thermal and nuclear quantum effect on nuclear dynamics and thermodynamic stability of this technologically important proton conducting perovskite oxide. The analysis suggests that the nuclear dynamics in hydrated BZCY72 is a result of a subtle interplay of harmonic, anharmonic and thermal effects, with the increased anharmonic character of the lattice dynamics above the orthorhombic to rhombohedral phase transition at 85 K. The anharmonic effect seems to be most pronounced in the case of oxygen and cerium. The analysis of the proton momentum distribution reveals that the concentration of the hydrogen in the BZCY72 lattice is constant across the orthorhombic to rhombohedral phase transition and further down to the room temperature. Moreover, the average hydrogen concentration obtained from our analysis of the mass-resolved neutron Compton scattering data seems to be commensurate with the total vacancy concentration in the BZCY72 framework. The calculation of the vibrational enthalpy of both phases allows obtaining the value of the enthalpy of the orthorhombic to the rhombohedral phase transition of −3.1 ± 1 kJ mol−1. Finally, our analysis of the nuclear kinetic energy of the proton obtained from NCS and the oxygen-oxygen distance distributions obtained from ND allows to conclude that BZCY72 in both the orthorhombic and rhombohedral phase at 70 K and 100 K respectively falls into the category of the KDP-type crystals where proton is probably under the influence of a double-well potential and forms hydrogen bonds of moderate strength. The obtained results have important ramifications for this technological important material.


Introduction
The concept of an economy based on hydrogen as a fuel is one of the most promising guidelines for future energy vectors [1][2][3]. Proton conducting electrolyte membranes are essential key materials for electrolysis, fuel cells and hydrogen pumps [4][5][6]. Acceptor-doped perovskites offer high protonic conductivity at an intermediate temperature range form 400°C-700°C, giving them thermodynamic and economic advantages over electrochemical cells operating at low temperatures [7]. In those structures, protons are present as foreign species formed by the incorporation of water by dissociative absorption replacing the oxygen vacancies: ( ) 2. Experimental 2.1. Sample preparation BZCY72 powder was obtained from Cerpotech, Norway. To ensure homogeneity, the powder was ball milled inside a zirconia container over a period of 12 h. Subsequently, the powder was calcined at 1200°C with heating and cooling rates of 10 K min −1 for 3 h in static air. Spark Plasma Sintering was carried out in the Tycho Sinterlab at the University of Rostock using HP D5 unit from FCT Systeme GmbH Rauenstein, Germany. About 8 g of BZCY72 powder was loaded into a graphite die with an inner diameter of 20 mm. Sintering was done at 1450°C for 5 min at an axial pressure of approximately 90 MPa. The heating and cooling rate was maintained at 150 K min −1 . An optical pyrometer focused on a small hole at the top punch was used for temperature monitoring.
After the sintering, the sample was separated from the die by means of a hydraulic press and additional heat treatment at 900°C for 2 h was necessary to remove the remaining graphite from the sample. Phase purity was confirmed by XRD using a Bruker D8 x-ray diffractometer. The BZCY72 sample was loaded into a tube furnace and hydrated at 680°C for 24 h under 3% H 2 O/Ar atmosphere. For the neutron scattering measurements, the cylindrical pellet was sealed into a square-shape (90 × 90 mm 2 ) aluminium container.

Neutron compton scattering
Neutron Compton scattering (NCS) was carried out at the VESUVIO time of flight spectrometer at ISIS Rutherford Appleton Laboratory, United Kingdom [38]. Measurements were performed on the sintered BZCY72 sample at 70 K, 100 K, 150 K, 200 K and 300 K. The details of the NCS data treatment are described in supplementary material available for this article. For the sake of clarity, we only mention here that, in what follows, when referring to the widths of recoil peak widths of protons and heavier nuclear species, the values of standard deviations of purely Gaussian peaks obtained in steps 9 and 5, respectively, of the traditional NCS data treatment procedure described in subsection 1.1 of the supplementary material (SI) is available online at stacks. iop.org/JPCO/4/045004/mmedia, will be used.
It is important to note at this point, however, that, as a matter of feasibility study yet another new to NCS data treatment procedure has been tested in this work and is outlined in the SI-the parametric bootstrap technique (PBT). Despite the existence of an extensive body of work on data treatment in neutron Compton scattering [27,29,30,32,36,[39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55], the subject of the treatment of systematic errors, parameter and error correlation and error propagation in Monte Carlo methods applied during subsequent steps of data analysis, described in subsection 1.1 of the SI, has largely been left aside. Moreover, a wider discussion on error propagation in NCS is also very much needed in the light of the recent developments concerning the substantial community growth and an upsurge in the scale and variety of scientific topics covered by the method, both triggering discussion about a substantial upgrade of the VESUVIO instrument in the future [28-31, 33, 56-61].
For the discussion on the performance of both methods of NCS data treatment, the reader is referred to subsections 1.1 and 1.2 of the SI. Here, we limit ourselves to noting that both techniques seem to give mutually consistent results on the test sample of NCS data on BZCY72 recorded at 70 K. Thus for the sake of brevity, the further discussion about NCS result will be held within the mainstream data analysis technique and more detailed aspects of the NCS data analysis, captured by the PBT method, will be outsourced to the SI.

Neutron powder diffraction
For details of the procedure of the VESUVIO neutron powder diffraction data treatment, the reader is again referred to the supplementary material. Here we will only limit ourselves to mentioning that, for the sake of comparison with the XRD data and benchmarking against the results of the ab initio simulations, only the ND data recorded at 70 K and 300 K were used.

Modelling
Modelling of BZCY72 was performed by means of the first-principles calculation based on density functional theory and the pseudo-potential method implemented in the CASTEP code version 19.1 [62]. The starting point of the calculation was the crystal structure of the BZCY72 loaded with hydrogen as specified in the CIF file obtained as a result of the refinement procedure of the powder neutron diffraction data recorded by the highresolution VESUVIO backscattering diffraction bank, a configuration that has already proven successful in the refinement of diffraction data [32,57,63,64]. The assignment of the positions of hydrogen atoms in the unit cell of BZCY72 was not attempted and the obtained CIF file contained only the information about the positions of the framework atoms. Thus, the effect of the presence of hydrogen in the structure was accounted for indirectly in the refinement, whereby the final, refined positions of the heavy atoms in the BZCY72 framework reflect the effect of the chemical lattice expansion and distortion due to the presence of hydrogen. This effect was very small and present chiefly in the case of oxygen atoms. In consequence, we have performed the lattice dynamics calculations only for the heavy BZCY72 framework atoms, leaving the task of simulating hydrogen dynamics in BZCY72 for future work and discussing the proton momentum distributions and kinetic energies purely at the level of experimental results. Such an approach seems feasible in the light of the fact that the nuclear dynamics of hydrogen in condensed matter systems are generally believed to be much less successfully described by harmonic lattice dynamics calculations based on the notion of classical (point-like) nuclei. Hydrogen nuclei are rather to be thought of as delocalized wave-packets on effective potential surfaces and their dynamics should be described within the framework of Born-Oppenheimer Molecular Dynamics (BOMD) with a quantum-thermostat or Path Integral Molecular Dynamics (PIMD) [42,56,[65][66][67][68][69].
Due to the disordered nature of BZCY72 structure [17], the supercell approach was adopted in modelling, as detailed in supplementary material. The supercell approach has enabled the calculation of both NCS-specific and thermodynamic observables. The traditional approach to computing thermodynamic observables, such as C T , V ( ) is based on the vibrational contribution calculated from the total VDOS. However, the nuclear massselective nature of the NCS method allows in principle to extend this approach by computing mass-selective variants of the molar heat capacity and Debye temperature by effectively replacing the total by partial VDOS by atom-projected (partial) VDOS of a given nucleus. This, in turn, allows fitting modified equations for C T V M , ( ) with an effective mass-dependent Debye temperatures, Q .
, ( )and Q D M , can then be benchmarked against the values of the NMD widths measured in an NCS experiment. In other words, despite the fact that those 'partial molar heat capacities' are not directly amenable to experimental scrutiny, one can still validate and benchmark them indirectly as they share the same underlying pVDOSes with the measured NMDs. In this context the NCS technique serves effectively as a 'mass-resolved nuclear quantum calorimeter', thereby providing a novel, previously untapped material characterization dimension.

Neutron diffraction
The powder neutron diffractograms of BZCY72 are shown in figures 1 and 2, for the samples measured below (at 70 K) and above (at 300 K) the orthorhombic (Imma) to rhombohedral (R c 3 ) phase transition respectively After the exclusion of the d-spacing regions between 1.38 and 1.48 and between 2.23 and 2.40 corresponding to the scattering of the aluminium sample container, the quality of the refinement is very good (Rwp=2%). As expected based on the pre-existing XRD data [17], both ND refinement procedures for BZCY72, performed above and below the orthorhombic (Imma) to rhombohedral (R c 3 ) phase transition, resulted in typical AMO3 perovskite structures made of corner-sharing MO6 octahedra (where M is the octahedrally coordinated cation in the AMO3 perovskite). Rietveld refinement, within agreement factor, Rwp, equal to 2.0% yielded two phases: • For the data recorded at 70 K, the orthorhombic phase with the space group Imma and lattice parameters a=5.990 Å, b=8.508 Å, and c=6.0193 Å; • For the data recorded at 300 K, the rhombohedral phase with the space group R c 3 and lattice parameters a=6.0246 Å, b=6.0246 Å, and c=14.724 Å.
In both cases, the refinement yielded mixed occupancies of 0.7, 0.2, and 0.1 for Zr, Ce and Y, respectively. In spite of numerous attempts to include the positions of the hydrogen atoms in the refinement of the ND data, no improvement in agreement factors has been achieved. Moreover, the obtained solutions had a tendency to locate hydrogen atoms pointing towards centres of the MO6 octahedra. This result seemed unphysical in the light of the results of the recent concurrent XRD-ND study on deuterated BZCY72 [17]. This work has shown that the low-temperature orthorhombic phase is characterised by two independent deuteron sites lying close to a plane bisecting the O−O octahedral edge and the rhombohedral phase has got only one deuteron site with the O−D bonds forming angles of 51°and 37°with the octahedral edges Thus, the results of our final Rietveld refinement did not include the hydrogen positions in the unit cells of BZCY72 both at 70 K and 300 K.
The comparison of the final arrangements of the O, Y, Zr, and Ba framework atoms, as obtained from the Rietveld refinement of the ND data of hydrated BZCY72 and XRD data of the pristine (free of hydrogen) BZCY72, is shown in figures 3 and 4, for the structures obtained from data recorded at 70 K and 300 K respectively. Root-mean-square deviations (RMDs) of atomic positions in the unit cells obtained from XRD and ND refinement of BZCY72 have been calculated. The result of this calculation is shown in tables 1 and 2 for data recorded at 70 K and 300 K, respectively. It is evident from tables 1 and 2 that the RMSD values are dominated by the contributions from oxygen atoms in both cases. The RMSDs of oxygen atoms are ca. five times bigger than the values for the other BZCY72 framework atoms at 70 K and ca. four times bigger at 300 K. The RMSDs of all heavier framework atoms are almost identical at both temperatures. Thus, the presence of the hydrogen in the BZCY72 lattice apparently causes only a very small amount of lattice distortion and expansion and a a slight displacement of the oxygen atoms placed at the corners of the octahedra. Taken together, the analysis of the structural data presented above allows to provide initial parameters (as specified in the respective CIF files) for the ab initio modelling of the material. Importantly, in the case of both crystallographic phases, the initial structures used during theoretical modelling are those specifying only the positions of the framework atoms in BZCY72 with hydrogen atoms positions not taken into account directly but rather through a subtle effect of lattice distortion reflected mainly in the displacement of the framework oxygen atoms.

Ab initio modelling
As shown in detail in the Supplementary Information, the total and the atom-projected VDOS distributions are largely insensitive to the choice of a particular randomly generated supercell. Thus, here only the results of calculations for one particular supercell for each structure will be presented here.
Our ab initio DFT calculations are performed within the harmonic approximation, with the electronic structure, the equilibrium geometry and vibrational properties computed in the zero-temperature limit. The    account for the temperature effect has been made by scaling the zero-temperature vibrational densities of states by the Boltzman population factor without accounting for the temperature-dependent lattice expansion. Thus, in what follows the DFT calculation performed for the orthorhombic structure optimised at 0 K will provide the base for the predictions of nuclear and thermodynamic observables in the entire range of temperatures below 85 K. The DFT calculation performed for the rhombohedral structure optimised at 0 K will form the base for predictions above 85 K down to 300 K. Figures 5 and 6 show the electronic band structures and the electronic densities of states (DOS) for the geometry-optimized orthorhombic and rhombohedral structures, respectively. The electronic structure is very similar in both cases with a bandgap of ca 1 eV visible in the DOS, reflecting the insulating nature of the system both below and above the orthorhombic to rhombohedral phase transition at 85 K.
The total VDOS and the pVDOS distributions for the orthorhombic and rhombohedral phase of BZCY72 are shown in figures 7 and 8, respectively. In both cases, the total VDOS distributions are dominated by two components, the pVDOS of Ba and O. The pVDOS of Ba almost entirely consists of a single band centred at 10 meV, which is of the same shape in the case of both structures. On the contrary, the pVDOS of O is distributed more evenly across the entire range of energies. The pVDOS distributions of Zr, Ce and Y are centred in the region of 20-30 meV but stretch towards 60 meV. This energy band becomes even more populated in the case of the rhombohedral structure, thereby shifting the centres-of-gravity of the pVDOS distribution of Ce, and, to a lesser extent, the pVDOS distributions of O, Zr, and Y towards slightly higher energies. Moreover, the appearance of vibrational states in the region of ca 60 meV of the simulated pVDOS of Ce, may be due to the partially reduced character of Ce atoms present in the BZCY72 framework, an effect already observed partially reduced ceria molecular clusters [70].
The molar heat capacity, C v , calculated from equation (7) of SI using as input the total VDOS, g E , ( ) and pVDOS distributions, g E , M ( ) shown in figure 7 (for temperatures below 85 K) and figure 8 (above 85 K and up     to 300 K), are shown in figure 9. As discussed above, the only temperature effects included in the calculations of the C v were the Boltzman population factors, scaling the zero-temperature limit pVDOS and VDOS distributions. Two immediate observations encourage further discussion of the plots shown in figure 9. Firstly, there is a marked jump of the value of the C v right at 85 K, the temperature of the phase transition. Secondly, the value of the C v obtained from the calculation at 300 K is 109. 6 [71].
Proceeding to the discussion about the molar heat capacities obtained from partial (atom-projected) VDOS distributions, shown in figure 9 as colour-coded shaded areas, we first notice one quite striking feature, namely the sudden jump in the value of C v in the case of Ce. The value of C v obtained for Ce jumps by ca. 7.1 [Joule/ mol/K], which is ca. 88% of the jump in the calculated value of the total C v at 85 K In order to shed more light on the observed behaviour of the C v calculated for individual nuclear species in BZCY72 we resort to yet another popular thermodynamic observable, the Debye temperature. A classical molecular dynamics study on the parent compound of BZCY72, the barium zirconate, has shown that the trend of the molar heat capacity as a function of temperature corresponds to the trend predicted by the Debye model, which approaches a classical limit above the Debye temperature [71]. Motivated by this result, we have conducted a similar analysis for the BZCY72 system, enriching it with the possibility of dissection of the C v into contributions from pVDOS distributions. The results of this analysis are shown as solid dashed lines in figure 9. Firstly, we note that the approximation of a single value of the Debye temperature does not work very well. The average of the values of Debye temperatures below and above the phase transition is 468.7 K. This value is smaller than the one calculated for barium zirconate (544 K) [74]. Moreover, it is quite apparent that the single temperature-independent Debye temperature model describes the C v values for Y, Zr and Ba quite well. Yet, the same model fails in the case of O and Ce. In the latter case, the disagreement is quite dramatic, leading to almost doubling of the Debye temperature above the phase transition. The Debye temperature can be interpreted as the temperature of a crystal's highest normal mode of vibration, and it correlates the elastic properties with the thermodynamic properties such as phonons, thermal expansion, thermal conductivity, specific heat, and lattice enthalpy [75]. Generally, compounds with large Debye temperature are hard with large sound velocity and have high thermal conductivity [75]. From less-known correlations, the enthalpy change of destruction of the crystal lattice of oxygen-containing salts during thermal decomposition is proportional to the Debye temperature of the metal element. In other words, reducing the Debye temperature leads to the lower stability of crystalline salts [76]. In the case of barium zirconate, the simulations have demonstrated maximum thermal conductivity at room temperature with a descending trend towards almost half its value at 2000 K [71], a result that matches very well with the experimental values obtained from Vassen et al [72]. Thus, above the orthorhombic to the rhombohedral phase transition, the thermal conductivity of BZCY72 most likely also increases. Apart from thermal stability and conductivity, much-increased values of the partial Debye temperature, calculated from the pVDOS of Ce and O above the phase transition, would suggest that in the rhombohedral crystal phase BZCY72 the local, effective binding strength of Ce and O is increased, thus shifting their respective highest normal modes of vibration towards higher energies. Increased degree of binding would then translate into increased localisation of Ce and O in the BZCY72 lattice in the rhombohedral crystal phase.

Neutron compton scattering 3.3.1. Framework BZCY72 atoms
We start our discussion of the results of the NCS experiments on BZCY72 by showing the data recorded in backscattering after the application of all necessary corrections, as per step 5 of the data correction procedure described in the Supplementary Information. The data recorded at 70 K and 300 K together with the respective fitting curves describing the shapes of individual recoil peaks is shown in figures 10 and 11, respectively. Due to kinematics of the neutron Compton scattering, the single scattering events off protons present in the sample are not present in backscattering. The inverse geometry of the VESUVIO instrument dictates the order in which the recoil peaks of atomic species of masses higher than hydrogen appear in the TOF spectrum [29,30,44,47,[59][60][61]77]. Namely, from left to right of the TOF spectrum axis one can clearly observe the recoil peaks of nuclear masses in increasing order: oxygen, aluminium from the sample container, yttrium, zirconium, barium and cerium. The mass resolution of the NCS technique can be approximately described as depending inversely on the square of the recoiling nuclear mass [29,30,44,47,[59][60][61]77]. Thus, whereas the mass resolution of O and Al and the group of heavier nuclei (Y, Zr, Ba, Ce) is clearly demonstrated, the resolution within this latter group of recoil peaks is much lower. However, owing to the application of the technique of stoichiometric fixing, the widths of the recoil peaks of all nuclei present in BZCY72 can be fitted with satisfactory precision. As evidenced by figure 12, the relative error of the Gaussian fits of the heavy masses in the backscattering regime, measured as the ratio of the single-standard deviation (1-STD) of the fitted peak width to this peak width value, is of the order of 7% or better. As expected, the values of the widths are greater for larger nuclear masses and increase with increasing temperature in an almost monotonic manner, thus reflecting the Boltzmann population of partial vibrational densities of states contributing to respective peak Doppler broadening, a situation commonly occurring in the case of heavyweight nuclear species in condensed matter systems [29,30,46,57,63,78,79].
Apart from stoichiometric fixing, listing the measured NMD widths and atomic kinetic energies of metallic Al (the sample container) at 70 K and 300 K and comparing them with theoretical predictions can serve as a selfconsistency check with values of these observables obtained for other atomic constituents of the sample, especially in the case of BZCY72 framework atoms. The measured NMD widths (atomic kinetic energies) at 70 K and 300 K of Al in the sample container were 9.8±0.6 Å −1 (22.2±2.7 meV) and 13.7±0.9 Å −1 (43.4±5.7 meV), respectively. These values can be contrasted with theoretical predictions based on the model of isotropic harmonic lattice and using the concept of the Debye temperature. Using the values of the Debye temperature for Al of 433 K in the low-temperature limit of 0 K and 390 K at room temperature [80], one obtains the NMD width (atomic kinetic energies) of 9.6 Å −1 (21.6 meV) and 13.4 Å −1 (41.9 meV), at 70 K and 300 K respectively. Thus, the experimental values agree within one standard deviation with the theoretical predictions. This result may thus serve as an encouragement for further analysis of NMD widths (atomic kinetic energies) of framework BZCY72 atoms.
The measured widths of oxygen momentum distributions seem to slightly depart from the trend given by DFT-based predictions calculated from pVDOS distributions for the orthorhombic and rhombohedral phases  . NCS data (black data points with error bars) from BZCY72 recorded at 300 K in backscattering after the application of all necessary corrections, as per step 5 of the NCS data correction procedure described in Supplementary Information. The recoil peaks present in the corrected NCS data were fitted using Gaussian profile functions, shown as shaded areas with labels pertaining to respective atomic species present in the sample. Additionally, the recoil peak fit of the aluminium (Al) sample container is plotted as a pink shaded area. The sum of the fitted recoil peaks is shown as a solid red line. On the experimental front, it would be desired to be able to distinguish oxygen types by their bond type and bond-distribution using the neutron Compton scattering technique. However, the present instrumental resolution and precision of the VESUVIO spectrometer do not allow for such a detailed analysis for atomic species heavier than hydrogen and further analysis must be carried out at the level of observed mean values of KE(O) and NMD widths.
In order to compare and contrast the DFT-predicted values and trends for the mean values of KE(O) and NMD widths of the oxygen with the experimental values, we mention first that, the DFT-predicted jump in these values at 85 K is indeed reflected in the experimental values obtained fro oxygen in the NCS experiments. What happens after the temperature of the phase transition is perhaps even more interesting. Namely, the measured NMD widths of the oxygen seem to start building up a trend, whereby values systematically departing from their respective theoretical predictions are building up as temperature increases from 100 K towards 300 K. The biggest departure is observed at 200 K. A systematic increase in the width of the momentum distributions, by virtue of the uncertainty principle, must be accompanied by a concomitant decrease in the respective nuclear position distributions, a sign of a systematically increasing degree of binding or confinement of the nuclear species in a solid-state system [29,30,44,47,[59][60][61]77]. Literature studies on BZCY-type proton ionic conductors [9] as well as the prototypic barium zirconate [71] strongly emphasize the special role played by oxygen in those systems. A detailed molecular dynamics study of the BZCY72 prototype, the barium zirconate, has shown that the typical interatomic potential well of Ba-O and Zr-O atomic pairs are due to the Coulombic attraction that could support bound states [71]. The authors conclude that further study of local oxygen dynamics in this system using the Schroedinger equation is required to relate these bound states to the observed properties [71]. Thus, our direct measurement of the width of oxygen momentum distribution in a closely related BZCY72 seems to be very much in order. It is the NCS method applied here that is in principle capable of solving the reverse problem of the nuclear Schroedinger equation in nuclear momentum space representation and reconstructing the shape of the local, effective potential of the mean force acting on a nucleus from the measured underlying nuclear momentum distribution [29][30][31]47]. The emerging picture of the local nuclear dynamics of the oxygen in BZCY72, as inferred from our NCS backscattering data, seems to reflect a subtle interplay between the delocalising in space (localising in momentum) nature of oxygen vacancy creation with temperature and space-confining (delocalising in momentum) nature of the local interatomic Ba-O and Zr-O potentials. This interplay, in combination with Boltzmann statistics, leads to a non-monotonic increase of the oxygen momentum distribution widths with temperature beyond the values predicted by the DFT calculation based on pVDOS distributions computed within the harmonic approximation.
Whereas a non-monotonic trend in the case of width of the momentum distribution of the oxygen is visible, similar trends for Zr and Ba are more difficult to assess as they may be partially obscured by the mass-dependent NCS resolution and the recoil peak overlap in the group of heavier nuclei (Y, Zr, Ba, Ce). Notwithstanding the above, an inspection of Zr and Ba panels in figure 12 yields qualitative assessment on how well the harmonic DFT prediction accounts for the observed behaviour of NMD widths of Y, Zr, Ba and Ce with a number of very interesting observations emerging. To start with, we mention that the Ba NMD widths seem to follow the DFT prediction. However, the Zr NMD widths follow a trend whereby their experimental values, starting from below their respective DFT-predictions, seem to systematically converge towards them as the temperature is increased towards 300 K. A conclusion would follow that degree of binding/confinement of Zr is increasing towards the room temperature. A similar, perhaps even more pronounced than in the case of Zr, the trend is visible in the case of Y.
Importantly, a trend whereby the values of the NMDs widths are observed to be systematically lower than the DFT predictions within the harmonic lattice approximation may signify the presence of anharmonicity. One has to bear in mind, however, that the anharmonicity is sensed by the NCS method in a slightly more subtle manner than it is the case for the infrared or Raman spectroscopy for that matter. Namely, as mentioned above, NMD widths can be termed as Boltzman weighted first moments of pVDOS distributions, thus emphasising the highenergy mode contribution to the pVDOS. Thus, the observed systematic shift of the NMD values with temperature may stem from increased anharmonic character of the high-energy modes in the BZCY72 lattice. This, in turn, may signify an increased anharmonic character of the local, effective potentials of these atomic species.
Typically, in a non-magnetic insulating material, the temperature dependence of a phonon energy comes from two primary sources: (i) the lattice contribution originating from the crystal thermal expansion induced changes in the harmonic force constants, and (ii) the anharmonic contribution arising from the coupling of phonon modes through the cubic and quartic anharmonicity [81,82]. The cubic anharmonicity term contributes to the phonon linewidth as well. The quartic anharmonicity to first order only contributes to the phonon frequency [83]. Both the shift and the broadening of the individual phonon modes would then propagate to pVDOS distributions, and their combined action will shift their centres-of-gravity and change the NMD width [83]. Usually, the action of the anharmonicity is identified with lowering (softening) of the mode energy through the cubic term with increasing temperature [81][82][83]. However, as evidenced by the trends for the Zr and Y nuclear species shown in figure 12, the reverse seems to be the case in the case of the BZCY72 system. Even more interestingly, in the case of Ce, a cross-over regime is observed, whereby mode hardening starts from 70 K with the centre-of-gravity of the Ce pVDOS systematically converging from below towards the level given by the harmonic approximation at 200 K and then continuing above 200 K with the centre-of-gravity of the Ce pVDOS distribution possibly shifting towards higher energy to reach values above the level given by the harmonic approximation. Adding to this striking behaviour of Ce, also the fact that, as discussed above in the case of the oxygen, a systematically higher than predicted by the harmonic lattice simulation degree of binding/ confinement is observed throughout the entire temperature range from 70 K to 300 K, a quite puzzling overall picture emerges.
It is instructive to re-examine the different contributions to the anharmonicity within the cubic and quartic term picture as well as the contribution stemming from the sheer lattice expansion. Mathematically, the cubic and quartic anharmonicity will contribute to the overall effect with opposite signs [81][82][83]. These two terms can cause a phonon mode to soften or harden as the temperature increases depending on their relative magnitudes [81][82][83]. Moreover, the strength of the two terms increases at different rates as the phonon frequency increases. Namely, for some modes, the magnitude of the quartic term can be larger than that of the cubic term, which results in mode hardening as the temperature rises. Conversely, when the cubic term is dominant mode softening is observed as the temperature is increased [81][82][83]. There are several physical systems in which this complicated behaviour is observed. One prominent example is served by a thermoelectric and topological class of materials where strong anharmonicity leads to low thermal conductivity due to enhanced phonon scattering [81][82][83]. As far as the lattice expansion with temperature is concerned, its contribution to vibrational mode shifting with temperature is always causing mode softening [81][82][83]. This fact renders the scenario depicted in figure 12 even more striking. Namely, the thermal expansion contribution to mode softening is often quite substantial (at the level between 30 and 60% thermoelectric and topological materials, for instance [83]). Thus, factoring out the expansion term and isolating the purely anharmonic contribution to mode shifting usually reveals even more pronounced mode hardening and/or softening trends [83]. Ideally, such factoring could be performed with the aid of known Grüneisen parameters, or by performing structure elucidation based on diffraction data collected as a function of temperature, followed by DFT work consisting of geometry optimisation and phonon calculation, thus replicating the computational protocol used for modelling in this work for a series of increasing temperatures. However, such a computational protocol, even when performed at a single temperature, presents itself as a formidable task, owing to the large degree of computational complexity related to accounting for the random crystal structure with mixed atomic occupancies when doing the DFT work. The modelling work presented here has thus seemed like a sound compromise between the degree of the qualitative and quantitative description provided and the more complicated, tedious task of disordered lattice structure and dynamics modelling at every given temperature has remained in the realm of the much-needed future work. Notwithstanding the above however and slightly fortunately perhaps, the picture shown in figure 12 does serve very well our characterisation of the degree and nature of local atomic binding in BZCY72. This is owing to the fact that by performing a lattice dynamics simulation for each crystal phase by means of DFT based on the assumption of no accounting for the lattice expansion within this phase and modelling the temperature effect merely by the Boltzmann population factors one naturally overestimates the magnitudes of the NMD widths. Thus, the departures from the computed trends shown in figure 12 can be, in such a computational setting, deemed as lower conservative bounds for the magnitudes of the anharmonic effects. In other words, had the lattice expansion in both crystal phases of BZCY72 been taken into account in calculating the DFT-based predictions of the NMD widths, the observed magnitudes of departures of experimental data from these trends would have been even larger, as the thermal expansion effect would have driven the lattice dynamics-based predictions in the opposite direction. The emerging picture of the degree and nature of average effective binding of nuclear species in BZCY72 lattice as a function of temperature shows marked signatures of anharmonicity in the case of all nuclear species in BZCY72, perhaps with the exception of Ba.
Importantly, the observed departures from the harmonic lattice dynamics predictions for most of the atomic species in BZCY72 appear very interesting in one more context. Namely, in proton-conducting electrolytes, there is an increasing conductivity contribution of oxide ions at higher temperatures and electrons or electron holes under reducing or oxidising conditions, respectively. In BZCY72, there is considerable p-type conductivity under oxidising conditions and oxide ion conductivity at high temperatures, in addition to the proton conductivity [13,84,85]. Moreover, the concomitant increase in the anharmonic character observed for Ce with increasing temperature seems to correlate very well with models and observations present in the literature suggesting that an increased level of cerium doping in BZCY is responsible for increased parameters, water uptake and increased both total and ionic conductivity [9,86,87]. On the whole, it appears that the signatures of anharmonicity, as viewed by oxygen and cerium NMDs measured by NCS in the range of temperatures between 70 K and 300 K, may signal the onset of anharmonicity-driven mixed conductivity at higher temperatures. Usually, the increased degree of anharmonicity leads to decreased material performance in terms of thermal and ionic conductivity at increased temperatures [83]. Thus, our observations based on the assessment of the degree and nature of local atomic binding based on the NCS data recorded in backscattering have potentially important ramifications for the modelling and understanding of proton conduction in BZCY72.

Protons in BZCY72
In order to perform a further critical appraisal of BZCY72 properties, we continue our discussion enriched by the information provided by the NCS forward scattering data. The presence of hydrogen atomic species is explicitly accounted for by the shape of the hydrogen NMD, shown in figures 13 and 14. The figures show isolated in the longitudinal momentum space, focused hydrogen momentum distributions in BZCY72 at 70 K and 300 K, respectively. As clearly seen in figure 15, there is a distinct jump in the values of the widths of the NMD of hydrogen across the phase transition at 85 K, indicating that the change of the symmetry of the crystal lattice from orthorhombic to rhombohedral induces an abrupt change in the degree of hydrogen binding/ confinement. According to a recent concurrent XRD-ND study on deuterated BZCY72, at low temperature (orthorhombic symmetry, 10−85 K) there are two independent deuteron sites, and in the rhombohedral domain, only one deuteron position exists [17]. Thus, the H widths observed below 85 K are expected to result from the averaging between the two types of proton sites and the H widths observed between 85 K and 300 K should correspond to momentum distributions of protons in one type of sites. Figure 16 shows an example of a mass-resolved NCS spectrum recorded for BZCY72 in forward scattering at the temperature of 70 K. The data presented in figure 16 has been corrected for multiple scattering and gamma background, according to the procedure described in steps 1-7 in section 1 of the supplementary material. A clear, well-resolved hydrogen recoil peak is visible in figure 16, enabling the quantitation of the hydrogen in BZCY72 as a function of temperature. Recalling that the integrated intensity of a recoil peak recorded in NCS is proportional to the product of the total bound scattering cross-section of a given nuclear isotope and the number of moles of the nuclear species per formula unit of a sample under investigation [29][30][31]44], one can obtain the hydrogen concentration in BZCY72 as a function of temperature from the ratios of the integral scattering intensities of hydrogen and any other havier-than-hydrogen nuclear species. For the purpose of the hydrogen quantitation in BZCY72 we have chosen to use the scattering integral intensities of barium. The reason for such a choice is twofold. Firstly, unlike Y, Ce and Zr, the occupation of barium atoms in the BZCY72 lattice is not statistical. Secondly, Ba is characterised by the second-largest (after the oxygen) integral intensity of all heavy nuclei in the framework of BZCY72 and the choice of the oxygen as a yardstick for the relative hydrogen concentration would have been related to an error due to the fact that the number of oxygen nuclei per formula unit of BZCY72 is not exactly equal to three. Figure 17 shows the number of moles of hydrogen per formula unit of BZCY72, n , H obtained as a function of temperature using the scattering intensities of H and Ba, fitted from mass-resolved NCS spectra recorded in forward scattering between 70 K and 300 K, and their respective total bound scattering cross-section values tabulated in the literature [88], reading of 82 and 3.38 barn, respectively.
The value of n H is stable across the orthorhombic to a rhombohedral phase transition, with the average, between 70 K and 300 K, of 0.22±0.01. This value is in accordance with concentrations of protons in acceptor doped BaZrO 3 and BaCeO 3 , typically being around 10-20 mol% when fully saturated at lower temperatures, and commensurate with the total concentration of vacancies [89].
The quantitation of hydrogen in BZCY72 can serve one more important purpose. Motivated by the recent work by Bourgeois et al [90], one can perform the calculation of the enthalpy of the phase transition of BZCY72 from the orthorhombic to the rhombohedral phase. In order to obtain this quantity, two assumptions are in order. Firstly, we will assume that, at temperatures equal or lower than 100 K, the vibration contribution to the formation enthalpy of a given phase of BZCY72 is entirely dominated by the vibrational zero-point energy (ZPE), which, in the harmonic approximation, is equal twice the nuclear kinetic energy. Secondly, we will assume that the contributions from the total electronic energy to the formation enthalpies of both the orthorhombic and the rhombohedral phase of BZCY72 are negligibly small. Tables 4 and 5 summarise the ZPE  values obtained from the measured values of the nuclear kinetic energies of all constituent nuclei of the hydrated BZCY72 for the orthorhombic phase at 70 K and rhombohedral phase at 100 K, respectively.
Based on these values, one can anticipate the enthalpy of the phase transition of BZCY72 from the orthorhombic to the rhombohedral phase of −52 meV±9 meV (−5.0±1 kJ mol −1 ). To the best of the authors' knowledge, this is the first time this thermodynamic quantity has been obtained for BZCY72.
Finally, one more aspect of the nuclear kinetic energy of the protons in both phases of the hydrated BZCY72 is worth mentioning. As mentioned in the Introduction, in BZCY72 a large redshift of the OH-stretch infrared absorption band is observed in BZCY72, suggesting that strong hydrogen bond interactions mediate a fast proton transfer process [21,22]. At first glance, this situation seems to be similar to the case of the superprotonic conductor material trirubidium hydrogen bisulfate, Rb 3 H(SO 4 ) 2 , where a huge redshift of the OH stretch occurs combined with a low KE value of the H-atom, KE(H)= 113 meV, as verified by NCS measurements [91] and by DFT calculation [92]. The values of KE(H) = 113 meV, reported in [91], may be contrasted with values of KE(H) obtained, in a harmonic approximation as KE(H)=0.5 ZPE(H), from the values of ZPE(H) listed in tables 3 and 4, being 147±4 meV for the orthorhombic and 173±6 meV for the rhombohedral phase of BZCY72, respectively. A general trend has been reported [92], that the KE(H) values of the M 3 H(SO 4 ) 2 (M=Rb, K) compounds, in which the hydrogen bonds are centro-symmetric, are much lower than those of the KDP-type Figure 16. A mass-resolved NCS spectrum of BZCY72 in forward scattering at 70 K. The black data points-corrected NCS data. The solid red line-the total fit to the data. The coloured shaded areas represent the contributions from individual nuclear mass species: (i) green for H, (ii) blue for the Al container, (iii) magenta for Ba, (iv) dark-yellow for Zr, (v) navy blue for Ce, (vi) violet for Y, and (vii) light-grey for O. See text for details. crystals, in a direct consistency with the oxygen-oxygen distance, being a measure of the HB strength. In such centrosymmetric cases, the KE(H) values are usually low and the proton potential is localized, characterized by a single well where no proton tunnelling occurs [92]. In such a case, the KE(H) value is in fact much smaller than that incorporating one hydrogen bond and a covalent bond, as occurs in water (KE(H) ∼152 meV) where the hydrogen bonding is to a neighbouring H 2 O molecule [92,93]. The KDP crystal was also experimentally studied using NCS Compton scattering [54]. KDP has a transition temperature, Tc=123 K, below which it is in a ferroelectric phase where the structure is tetragonal. Above Tc, it is in a paraelectric phase whose structure is orthorhombic. The KE(H) in KDP was deduced to be 139.3 meV and 138.6 meV at 90 K and 130 K respectively [53]. Thus the KE(H) was found remaining almost constant across the para to the ferroelectric phase transition. It is worth mentioning that, in KDP, the origin of ferroelectricity is attributed to the proton off-centring and the   proton potential is more likely to have a double-well of an asymmetric shape [92].  [92]. As ROO decreases, the local, effective potential, felt by the proton was predicted to evolve from a double-well with a high barrier (compared to kT) to a single-well for the shortest bonds [94][95][96]. A critical ROO distance, Rc, can be defined below which the local potential for the proton is predicted theoretically to be in a single well potential [95,96]. This critical distance, as predicted by a two-state diabatic model by McKenzie et al [95], is in the order of 2.3 to 2.4 Å, depending on the coupling strength between the diabatic states describing two limiting cases of a proton completely transferred either at the site of the proton donor or the site of the proton acceptor. Interestingly, this theoretical prediction is in a very good agreement with the ROO value measured in Rb 3 H(SO 4 ) 2 [97] and the DFT-based prediction for the same system [92]. Moreover, in the NCS work on the KDP [91], it was concluded that the proton is probably under the influence of a double-well potential with the value of the KE(H). It is thus tempting to discuss the case of the proton in BZCY72 In the context of the correlations between the hydrogen bond strength, the hydrogen bond length, and the kinetic energy of the proton in the hydrogen bond widely reported in the literature and mentioned above. To this end, a statistical analysis of bond-length distributions in BZCY72 has been performed using the CrystalMaker software using as input the atomic coordinates corresponding to the structures of the orthorhombic and rhombohedral phase of BZCY72 refined from our neutron diffraction data presented above. The distribution of the O-O distances for both the orthorhombic phase at 70 K and in rhombohedral phase at 100 K was found to be bi-modal with two types of oxygen-oxygen distances distributed around peaks at values of 2.998 Å and 3.033 Å, thus with minimal change with temperature between 70 K and 100 K. Using the two diabatic state model by McKenzie, the proton in BZCY72, both in the orthorhomic and rhomohedral phase at 70 K and 100 K respectively, should be classified as participating in hydrogen bonds of moderate strength (ROO>2.6 Å). Moreover, these ROO values would correspond to the case of a proton experiencing a doublewell local, effective potential in the hydrogen bond. Thus, using the classification by McKenzie our result would suggest that BZCY72 in both the orthorhombic and rhombohedral phase at 70 K and 100 K respectively falls into the cathegory of the KDP-type crystals. This conclusion is even more corroborated by the fact that our KE(H) values, being 147±4 meV for the orthorhombic and 173±6 meV for the rhombohedral phase of BZCY72, at 70 K and 100 K respectively, are relatively close to the value of KE(H) in KDP, found to be 139.3 meV and 138.6 meV at 90 K and 130 K respectively [53].

Conclusions
Our combined neutron data analysis, augmented with density functional theory modelling of lattice dynamics, has enabled, for the first time, a mass-selective view of the combined thermal and nuclear quantum effect on the local, effective binding of atomic species present in hydrated BZCY72. We have observed systematic departures from the harmonic lattice-dynamics predicted trends for the widths of nuclear momentum distributions as a function of temperature, with both systematic decrease and increase of the observed NMDs above and below the theory-predicted levels, with the increased anharmonic character of the lattice dynamics above the orthorhombic to rhombohedral phase transition at 85 K. This scenario involves a subtle interplay between mode hardening and softening with increased temperature. The anharmonic effects seem to be most pronounced in the case of oxygen and cerium. Oxygen NMD widths are consistently observed to indicate a substantially increased degree of binding or confinement, compared to the harmonic lattice prediction, across the entire observed temperature range. If entirely caused by anharmonicity, this trend may be reconciled with systematic hardening of the local effective anharmonic potential experienced by the oxygen. This anharmonicity in the local oxygen environment may constitute an important ingredient of the origins of the onset in oxygen conduction observed at higher temperatures. The departure from the harmonic behaviour observed for the cerium in BZCY72 has got an even more complicated character, compared to the oxygen case. Namely, a cross-over regime has been observed, whereby at around 200 K initially lower, compared to the harmonic lattice prediction, degree of binding/confinement of the cerium species seems to gradually increase above the harmonic lattice prediction. This mixed behaviour may partially responsible for the high sensitivity of the thermal lattice expansion of BZCY72 on the level of Ce doping. This, in turn, may also explain why the mixed ionic conductivity of BZCY72 also depends on the Ce doping level in such a sensitive manner. Our analysis of the proton momentum distribution reveals that the concentration of the hydrogen in the BZCY72 lattice is constant across the orthorhombic to rhombohedral phase transition and further down to the room temperature. Moreover, the average hydrogen concentration obtained from our analysis seems to be commensurate with the total vacancy concentration in BZCY72 framework. Based on the obtained values of vibrational ZPEs for all BZCY72 constituent nuclei one can estimate the enthalpy of the phase transition of −3.1±1 kJ mol −1 . To the best of the authors' knowledge, this is the first time this thermodynamic quantity has been obtained for BZCY72. Finally, our analysis of the KE of the hydrogen obtained from NCS and the oxygenoxygen distance distributions obtained from ND allows to conclude that BZCY72 in both the orthorhombic and rhombohedral phase at 70 K and 100 K respectively falls into the category of the KDP-type crystals where proton is probably under the influence of a double-well potential and forms relatively weak hydrogen bonds.
Undoubtedly, more high-precision diffraction and spectroscopic work, especially above the room temperature, and ideally augmented by detailed ab initio modelling, is highly desired in order to fully understand the structure-function relationship in hydrated BZCY72 with the ultimate goal of being able to engineer the conductivity of this technologically important system.