# Solubility and diffusivity of hydrogen and its isotopes in the BeO system

Aurélien Riot, François Virot, Damien Connétable

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

## To cite this version:

Aurélien Riot, François Virot, Damien Connétable. Solubility and diffusivity of hydrogen and its isotopes in the BeO system. Journal of the American Ceramic Society, 2023, 106 (8), pp.5005-5021. 10.1111/jace.19130. hal-04067611

HAL Id: hal-04067611

https://hal.science/hal-04067611

Submitted on 18 Jul 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Solubility and diffusivity of hydrogen and its isotopes in the BeO system

Aurélien Riot<fnref n="1" />,<fnref n="2" />

François Virot<fnref n="2" />

Damien Connétable<fnref n="1" /> ID

## Correspondence

Damien Connétable, CIRIMAT, UMR 5085, CNRS-INP-UPS, ENSIACET, 4 allée Emile Monso, BP 44362, F-31030 Toulouse Cedex 4, France.

Email: damien.connetable@ensiacet.fr

## Funding information

CALMIP, Grant/Award Number: 2021-p0749; GENCI-TGCC, Grant/Award Numbers: A0080911483, A0100911483;

EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Program, Grant/Award Number: 101052200 - EUROfusion

## Abstract

This work revisits hydrogen insertion in the wurtzite beryllium oxide in order to fill the gap in scientific knowledge regarding the diffusion coefficients of the different chemical states of hydrogen (neutral, charged, or molecular). Both first-principle and macroscopic models were used to this end. In the former, two exchange-correlation functionals (PBE and SCAN) were used to compute the properties of interest and accuracy was then discussed. Regarding the behavior of interstitial insertion of hydrogen, this work is slightly different from previous works; hydrogen was mainly found in the form of charged $H^-$ or $H^+$ ions, depending on experimental conditions. In regard to diffusivity properties, a complex migration pathway was found for hydrogen cation. Finally, the present study succeeded in producing a reliable set of diffusion coefficients for neutral, charged, and molecular hydrogen.

## KEYWORDS

BeO, diffusivity, hydrogen, solubility, tritium

# 1 | INTRODUCTION

Beryllium oxide (BeO), known as beryllia, crystallizes under a hexagonal wurtzite structure at atmospheric pressure and room temperature. It is therefore unique as compared to all other alkaline-earth oxides, which have a rock salt structure. This ceramic has many interesting properties. For example, because of its high melting point (2843 K<fnref n="1" />) and chemical resistance, beryllia is often used as a crucible. In addition, it combines high thermal conductivity (290–330 W/m/K at 298 K,<fnref n="2" /> slightly lower than copper) with high electrical resistivity, which makes this material a promising insulator or heat sink for the electrical and microelectronic industries. Furthermore,

its dielectric constant is interesting (higher than that of SiO₂) and its interfacial thermodynamic stability with silicon is good. This may lead to opportunities in achieving the challenges related to the downscaling of field effect transistors.

In the nuclear field, its large neutron cross-section makes beryllium oxide suitable for applications as a neutron moderator and reflector.³ Moreover, this ceramic has also been suggested as an additive substance to uranium dioxide in the accident tolerant fuel (ATF) concept.⁴ Because of its interactions with neutrons, a possible application as a tritium breeder blanket in the field of nuclear fusion was explored. Beryllium oxide has also been suggested as a permeation barrier applied

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2023 The Authors. Journal of the American Ceramic Society published by Wiley Periodicals LLC on behalf of American Ceramic Society.

<footnote n="1">¹CIRIMAT, UMR 5085, CNRS-INP-UPS, ENSIACET, 4 allée Emile Monso, Toulouse Cedex 4, France</footnote>

<footnote n="2">²Institut de Radioprotection et de Sûreté Nucléaire (IRSN, PSN-RES/SAM), Saint Paul Lez Durance, France</footnote>
for hydrogen isotopes.⁵ Furthermore, in ITER facilities, plasma-facing components will be made of metallic beryllium blocks because of its interesting properties like its oxygen affinity which causes the formation of a beryllium oxide layer thereby purifying the plasma⁶

In both the microelectronic and the fusion fields, the behavior and properties of hydrogen insertion in beryllium oxide are of importance. In the former, hydrogen could have a detrimental effect on the gate insulator,⁷ while in the latter, tritium retention is one of the main safety issues in the ITER project. Indeed, due to the radiological consequences in the event of a loss of confinement in this type of facility, a safety limit of 700 g of tritium in the vacuum vessel has been established. As ITER will be fueled with deuterium and tritium associated with a low combustion efficiency, a part of the fuel can be trapped in the plasma-facing components made out of beryllium. Consequently, the tritium trapped in the wall should be periodically removed by performing a baking process.⁸ Knowledge on tritium behavior in beryllium has been improving recently from a theoretical point of view.⁹⁻¹¹ Since beryllium walls can be oxidized, the behavior of tritium in BeO should be first assessed by studying elementary properties like solubility and diffusivity.

Even though there are numerous data on beryllium oxide properties, those related to hydrogen isotopes in BeO are unfortunately scarce, especially in the experimental literature. The work of Fowler et al.¹² can be considered as a reference regarding tritium diffusivity in the different states of beryllium oxide (single-crystal, sintered and powdered); since the end of the 1970s, no additional experimental investigation has been performed on this topic. Moreover, as regards the solubility of hydrogen isotopes, the only experimental estimate was carried out on sintered beryllium oxide samples loaded with deuterium.¹³ Nevertheless, a few theoretical DFT-based studies investigated the insertion or diffusion properties of hydrogen in the wurtzite structure.¹⁴⁻¹⁷ The recent work of Hodille¹⁵ extends the work of Marinaupoulos et al.¹⁶ by publishing the solution energies for each chemical state of hydrogen (neutral, charged, molecular) but only the diffusion coefficient of molecular hydrogen was evaluated,¹⁵ while the diffusivity of neutral hydrogen was previously computed by Allouche et al.¹⁴ This work aims at filling the gap in scientific knowledge on the diffusion coefficients of hydrogen isotopes in beryllium oxide by means of DFT computation. It also seeks to compute their solution energies in order to perform an estimate of their solubility.

The paper is organized as follows. Section 2 introduces the computational approach. Section 3 presents BeO properties and discusses the improvements obtained from using the strongly constrained and appropriately normed (SCAN) functional. Sections 4 and 5 are dedicated to

results and discussions on the solubility and diffusivity of hydrogen isotopes.

## 2 | METHODOLOGY

Density functional theory (DFT)¹⁸,¹⁹ calculations were performed using the Vienna Ab initio Simulation Package (VASP).²⁰ Simulations were carried out under the Perdew–Burke–Ernzerhof (PBE)²¹ generalized gradient approximation for the exchange-correlation functional, as well as under the SCAN semi-local density functional.²² Magnetism (spin polarized) was taken into account in calculations. Projector augmented wave (PAW) pseudopotentials²³ were used, and a kinetic energy cut-off of 720 eV for the plane-wave basis set was chosen to compute energies (convergence criterion of 1 meV/atom). Energy calculations were performed on 3×3×3 supercells, i.e., 108 atoms. The corresponding k-mesh grids were then sampled with Γ-centered 2×2×2.²⁴ To optimize ideal structures, atomic coordinates and lattice parameters were fully relaxed.

To study uncharged hydrogen insertion (H and H₂), lattice parameters were fully relaxed. For H⁺ and H⁻, pure BeO lattice parameters were selected and the residual stress was corrected as described below. The formation energy of hydrogen, E<sub>f</sub>(H, q), in the charge state q, is computed as:

$$E_f(H, q) = E_o[\text{BeO}+\text{H}] - E_o[\text{BeO}] - \mu_H + q(\mu_e + E_{\text{VBM}}) + E_{\text{elas}} + E_{\text{elec}}, \quad (1)$$

where E<sub>o</sub>[BeO+H] and E<sub>o</sub>[BeO] are the DFT energies of the BeO supercell, respectively, with and without one H atom. μ<sub>H</sub> is the chemical potential of the H atom in its reference state, equal to 1/2E<sub>o</sub>[H<sub>2</sub>] or E<sub>o</sub>[H<sub>2</sub>] for atomic or molecular hydrogen, respectively; E<sub>o</sub>[H<sub>2</sub>] being the DFT energy of H<sub>2</sub>. μ<sub>e</sub> + E<sub>VBM</sub> is the energy of the Fermi level in the gap. Finally, E<sub>elas</sub> and E<sub>elec</sub> correspond to the elastic and electrostatic energies due to the periodic images. E<sub>elas</sub> was computed using the Aneto code, developed by Varvenne et al.,²⁵ from elastic constants of the BeO system and the residual stress when hydrogen is inserted (calculations were thus performed at fixed lattice parameters). The image-charge correction scheme proposed by Freysoltd et al.²⁶ was used to compute E<sub>elec</sub>. It uses planar-averaged electrostatic potential for aligning the defect-induced potential, which is obtained by subtracting bulk supercell potential from defective supercell potential, to the model charge potential.

To study migration paths, five intermediate positions were used for each nudged elastic band (NEB) calculation.²⁷ The energies of the transition states were
<figure><img src="image_3.png" /><figcaption>FIGURE 1 (Color online) Schematic representation of the BeO structure: Be and O atoms are in green and red, respectively. Gold balls represent the interstitial hydrogen sites considered here. The figure on the left (right) is a lateral (top sliced) view of the primitive cell. Figures were built with VESTA.³⁰.</figcaption></figure>

then calculated using the same criteria of convergence as for insertion sites. Additional calculations were carried out to compute interatomic force constants (IFC). The finite displacement method using the *phonopy* package²⁸ was performed on 2×2×2 supercells, that is, 32 atoms. Helmholtz free energies, $F_v$,

$$F_v(T) = k_B T \sum_{\nu=1}^{3N} \int n_{\nu}(\mathbf{q}) \ln \left[ 2 \sinh \left( \frac{\hbar \omega_{\nu, \mathbf{q}}}{2k_B T} \right) \right] d\mathbf{q}, \quad (2)$$

were then calculated. $\omega_{\nu,q}$ the frequencies in wave-vector **q** and mode $\nu$, and $n_{\nu}(\mathbf{q})$, the phonon density of state, were computed on 20×20×20 **q**-mesh grids. Frequencies were used to analyze the stability of configurations and Helmholtz free energies were used to compute jump rates, $\Gamma_{ab}$, zero-point energies (ZPE) and Gibbs insertion energies.

# 3 GROUND-STATE PROPERTIES OF THE BEO SYSTEM

## 3.1 Structural properties

The BeO system is in the hexagonal space group 186 (hP4, P63/mc), a wurtzite structure. Wurtzite is defined as a nesting of two compact hexagonal structures of Be and O atoms, which corresponds to an ABAB stacking of hexagonal planes, see Figure 1. BeO has two Be atoms and two O atoms per primitive unit cell, where both atoms occupy 2b Wyckoff sites (1/3, 2/3, z).²⁹ The structure of BeO can thus be defined by two lattice parameters *a* and *c*, and one parameter, $z_O$, which corresponds to the shift of an O-plane with respect to the Be plane, that is, $z_{Be} = 0$ and $z_O \approx 0.378$.

<table><caption>TABLE 1 Lattice parameters, *a* and *c*, in Å, and $z_O$ parameter of the oxygen in the BeO system.</caption><thead><tr><th></th><th></th><th>a (Å)</th><th>c (Å)</th><th>c/a</th><th>z<sub>O</sub></th></tr></thead><tbody><tr><td>Exp.</td><td>31–36</td><td>2.698</td><td>4.379</td><td>1.623</td><td>0.378</td></tr><tr><td rowspan="3">Theo.</td><td>PBE<sup>15,16,37–41</sup></td><td>2.726</td><td>4.417</td><td>1.620</td><td>0.378</td></tr><tr><td>HSE06<sup>16</sup></td><td>2.680</td><td>4.350</td><td>1.623</td><td>0.378</td></tr><tr><td>LDA<sup>31,42–44</sup></td><td>2.688</td><td>4.366</td><td>1.624</td><td>0.378</td></tr><tr><td></td><td>PBE†</td><td>2.713</td><td>4.404</td><td>1.623</td><td>0.378</td></tr><tr><td></td><td>SCAN†</td><td>2.686</td><td>4.363</td><td>1.624</td><td>0.378</td></tr></tbody></table>

*Note:* The data obtained in this work are marked with a † symbol. LDA, Local Density approximation; PBE, Perdew–Burke–Ernzerhof; SCAN, strongly constrained and appropriately normed.

Compared to the ideal wurtzite structure ($c/a = 1.633$ and $z_O = 0.375$), the BeO system has a lower $c/a$ ratio and a higher $z_O$ parameter, see Table 1.

Theoretical and experimental values are also summarized in this table. Experiments³¹⁻³⁶ give the same values, $a = 2.698$ Å, and $c = 4.379 \pm 0.003$ Å. Theoretical studies using the PBE functional show, at 0 K, slightly overestimate *a* and *c* values (about 0.6%), whereas the HSE06 or the LDA functionals tend to underestimate them, which is consistent with other literature results¹⁵,¹⁶,³⁷⁻⁴¹ and are thus in excellent agreement with experimental data. The SCAN (Meta-GGA) functional tends to slightly underestimate lattice parameters, as is the case for hybrid or LDA functionals.

Cohesive, $E_{coh}$, and formation, $E_f$, energies of the BeO system were then computed. The energy of formation corresponds to the difference between the energy of the crystal and the energy of its constituents under standard pressure and temperature conditions. $E_{coh}$ is thus given by:

$$E_{coh} = E_o[\text{BeO}] - E_o^{\text{at}}[\text{Be}] - E_o^{\text{at}}[\text{O}], \quad (3)$$

where $E_o^{\text{at}}[\text{Be}]$ and $E_o^{\text{at}}[\text{O}]$ are the DFT energies of atoms. The formation energy is given by

$$E_f = E_o[\text{BeO}] - E_o[\text{Be-hcp}] - \frac{1}{2}E_o[\text{O}_2], \quad (4)$$

where the reference states are the Be-hcp system and $O_2$ molecule. ZPE was also included in $E_{coh}$ (see the discussion on vibrational properties below). Theoretical and experimental data are gathered in Table 2.

The formation energy, taking into account the ZPE, was found equal to −2.73 eV/atom with the PBE functional and equal to −3.10 eV/atom with the SCAN functional. The
<table><caption>TABLE 2 Cohesion energy, $E_{coh}$, in eV per atom, and formation energy, $E_f$ in eV per atom, without and with zero-point energy correction (ZPE): comparison of theoretical and experimental data.</caption><thead><tr><th></th><th></th><th>E<sub>coh</sub></th><th>E<sub>f</sub> (+ZPE)</th></tr></thead><tbody><tr><td>Exp.</td><td>0 K</td><td>6.1 ± 0.1<sup>45</sup></td><td>-3.13<sup>46</sup></td></tr><tr><td rowspan="4">Theo.</td><td>PBE<sup>37,38</sup></td><td>5.2-5.5</td><td></td></tr><tr><td>PW91<sup>43</sup></td><td>6.5</td><td></td></tr><tr><td>LDA<sup>44,47-49</sup></td><td>6.8-7.9</td><td></td></tr><tr><td>PBE†</td><td>6.34</td><td>-2.77 (-2.73)</td></tr><tr><td></td><td>SCAN†</td><td>6.39</td><td>-3.15 (-3.10)</td></tr></tbody></table>

PBE, Perdew-Burke-Ernzerhof; SCAN, strongly constrained and appropriately normed. †: Present work.

PBE functional overestimates the experimental value of -3.13 eV/atom by 0.4 eV,⁴⁶ which corresponds to a difference of 13%. This value is an improvement compared to the energy previously calculated by Byggmästar et al.³⁷ with interatomic potentials. The value obtained with the SCAN functional is, on the other hand, in better agreement with experimental data with an error below 1%.

The cohesive energies calculated with the two functionals are closer to experimental values than previous studies. Here, differences of about 5% were found, whereas the values reported in the literature show larger differences of about 10%, even 30% in some cases, see Table 2. Contrary to the formation energy, the cohesion energy found using the PBE functional is closer to the experimental value than the one found using the SCAN functional. However, the difference remains small. Indeed, the SCAN functional allows a better description of the molecules than the PBE functional, but improves only slightly the description of isolated atoms.²²

## 3.2 | Electronic properties

Figure 2 represents the calculated electronic band structure of the BeO wurtzite which exhibits a direct band gap at Γ point, in agreement with experimental result⁵⁰⁻⁵² and theoretical calculations.¹⁶,⁵³⁻⁵⁶ Results are summarized in Table 3.

As done by Duman et al.,⁴¹ valence bands are separated into two groups. The first one, between -16 and -19 eV, corresponds to the 2s orbitals of the O atom with a small contribution from the 2p and 2s orbitals of beryllium, and not only from 2s as indicated in Ref. [41]. The second group, between -7 and 0 eV (top of the valence band), corresponds to the 2p orbitals of oxygen, with a weak contribution of the 2p orbitals of Be atoms. Apart from the difference between the 2s and 2p orbitals of beryllium, results obtained in this work are similar to those presented by Duman et al.⁴¹

<figure><img src="image_4.png" /><figcaption>FIGURE 2 Electronic band structure of the BeO wurtzite (PBE functional) along high symmetry points of the first Brillouin zone. The Fermi level was set at the top of the valence band, and $\Delta E_F = E_F - E_{vbm}$. The 4.2 eV value (blue line) corresponds to a shift of the Fermi level position induced by the assumption of the electrically neutral system for charged hydrogen solubility, whereas the other values (red and orange lines) are selected for the illustrative solubility calculations. PBE, Perdew-Burke-Ernzerhof.</figcaption></figure>

<table><caption>TABLE 3 Electronic band gap, $E_{gap}$ in eV, of the BeO system.</caption><thead><tr><th></th><th></th><th>E<sub>gap</sub></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="3">Exp.</td><td>50</td><td>10.6</td><td>Direct</td><td>Γ</td></tr><tr><td>52</td><td>10</td><td>-</td><td>-</td></tr><tr><td>51</td><td>9.9</td><td>-</td><td>-</td></tr><tr><td rowspan="6">Theo.</td><td>PBE<sup>16,53</sup></td><td>8.4/7.4</td><td>Direct</td><td>Γ</td></tr><tr><td>HSE06<sup>16</sup></td><td>10.6</td><td>-</td><td>-</td></tr><tr><td>SIC<sup>54</sup></td><td>10.5</td><td>Direct</td><td>Γ</td></tr><tr><td>B3LYP<sup>55</sup></td><td>10.4</td><td>Direct</td><td>Γ</td></tr><tr><td>GW<sup>56</sup></td><td>10.7</td><td>Direct</td><td>Γ</td></tr><tr><td>PBE†</td><td>7.45</td><td>Direct</td><td>Γ</td></tr><tr><td></td><td>SCAN†</td><td>7.5</td><td>Direct</td><td>Γ</td></tr></tbody></table>

Note: PBE, Perdew-Burke-Ernzerhof; SCAN, strongly constrained and appropriately normed. †: Present work.

The calculated electronic band gap, equal to about 7.45 eV with the PBE functional, is in agreement with previous theoretical results.¹⁵,¹⁶ A similar value is obtained with the SCAN functional (approximately 7.5 eV). For the band gap energy, this functional is still worse than the HSE06 hybrid functional, which gives more accurate experimental values.¹⁶,¹⁷ These values are 30% lower than the experimental value; a well-known phenomenon due to the fact that DFT is a ground-state theory which tends to underestimate the value of the experimental band gap, which is equal to 10.6⁵⁰,⁵² or 9.9 eV.⁵¹ This “low” value of the gap should not affect our results since in all the cases studied here, the insertion of H leads to bands that are always located in the gap (large enough).
<table><caption>TABLE 4 Bulk modulus, *B*, and elastic constants, *C*<sub>ij</sub> of beryllium oxide.</caption><thead><tr><th></th><th></th><th><i>B</i></th><th><i>C</i><sub>11</sub></th><th><i>C</i><sub>12</sub></th><th><i>C</i><sub>13</sub></th><th><i>C</i><sub>33</sub></th><th><i>C</i><sub>44</sub></th><th><i>C</i><sub>66</sub></th></tr></thead><tbody><tr><td rowspan="5">Exp.</td><td><sup>31</sup></td><td>208</td><td>454</td><td>85</td><td>77</td><td>488</td><td>155</td><td>184</td></tr><tr><td><sup>57</sup></td><td>242</td><td>469</td><td>131</td><td>119</td><td>499</td><td>159</td><td>169</td></tr><tr><td><sup>57</sup></td><td>224</td><td>461</td><td>127</td><td>89</td><td>492</td><td>147</td><td>167</td></tr><tr><td><sup>58</sup></td><td>-</td><td>461</td><td>127</td><td>89</td><td>492</td><td>148</td><td>167</td></tr><tr><td><sup>59</sup></td><td>252</td><td>470</td><td>168</td><td>119</td><td>494</td><td>153</td><td>151</td></tr><tr><td rowspan="2">Theo.</td><td>PBE<sup>37–39,41</sup></td><td>206–218</td><td>424–452</td><td>97–134</td><td>72–97</td><td>466–483</td><td>127–155</td><td>145–178</td></tr><tr><td>LDA<sup>42,60,61</sup></td><td>204–223</td><td>411–472</td><td>116–129</td><td>85–121</td><td>446–530</td><td>126–138</td><td>143–154</td></tr><tr><td></td><td>PBE<sup>†</sup></td><td>211</td><td>427</td><td>115</td><td>82</td><td>463</td><td>137</td><td>156</td></tr><tr><td></td><td>SCAN<sup>†</sup></td><td>230</td><td>470</td><td>129</td><td>90</td><td>510</td><td>150</td><td>170</td></tr></tbody></table>

Data are expressed in GPa. $C_{66} = (C_{11} - C_{12})/2$. PBE, Perdew-Burke-ErnzerhofA; SCAN, strongly constrained and appropriately normed. †: Present work.

## 3.3 | Elastic properties

For hexagonal structures, there are five independent elastic constants, that is, $C_{11}$, $C_{12}$, $C_{13}$, $C_{33}$, and $C_{44}$. To compute them, a stress-strain method was used. The primitive cell (with high convergence criteria: 800 eV and, 18×18×16 k-mesh grids) was thus employed. The bulk modulus, *B*, was deduced from these values using the following expression:

$$B = \frac{2C_{11} + C_{33} + 2C_{12} + 4C_{13}}{9} \quad (5)$$

and $C_{66}$ deduced from $C_{66} = (C_{11} - C_{12})/2$. Table 4 summarizes the different data of *B* and $C_{ij}$ available for the BeO system.

Our results show a good agreement between calculated and experimental data.<sup>31,57–59</sup> The calculated constants are in the same order of magnitude as the other theoretical results, with similar differences as found in the experimental data. Similarly, for the incompressibility modulus, the calculated value agrees with the other theoretical values, while being in good agreement with the value of Bosak et al.<sup>31</sup>

For elastic properties, the SCAN approximation also produces results that are consistent with the experimental data.

## 3.4 | Vibrational properties

Vibrational properties were computed using the *phonopy* code.<sup>28</sup> Figure 3 shows the vibrational band structure calculated with PBE and SCAN functionals. To include splitting of LO-TO modes, the Born charge and dielectric tensor were included in the calculations (with PBE only). Values are reported in Table 5. They were computed using the density functional perturbation theory (DFPT) as implemented in the VASP.<sup>62,63</sup>

<figure><img src="image_5.png" /><figcaption>FIGURE 3 Phonon dispersion curves (along high symmetry points) calculated with the SCAN (full black line) and PBE (dashed red line) functionals in comparison with the experimental data of Bosak.<sup>31</sup> PBE, Perdew-Burke-Ernzerhof; SCAN, strongly constrained and appropriately normed.</figcaption></figure>

<table><caption>TABLE 5 Dielectric constants and effective Born charges of the BeO system.</caption><thead><tr><th></th><th></th><th>ε<sup>∞</sup><sub>⊥</sub></th><th>ε<sup>∞</sup><sub>∥</sub></th><th>Z<sup>Be</sup><sub>⊥</sub></th><th>Z<sup>Be</sup><sub>∥</sub></th></tr></thead><tbody><tr><td rowspan="2">Exp.</td><td><sup>64</sup></td><td>2.95</td><td>3.00</td><td>-</td><td>-</td></tr><tr><td><sup>65</sup></td><td>2.95</td><td>2.99</td><td>-</td><td>-</td></tr><tr><td rowspan="2">Theo.</td><td><sup>41</sup></td><td>3.06</td><td>3.13</td><td>1.80</td><td>1.86</td></tr><tr><td><sup>31</sup></td><td>3.06</td><td>3.12</td><td>1.79</td><td>1.85</td></tr><tr><td></td><td>PBE<sup>†</sup></td><td>3.13</td><td>3.19</td><td>1.79</td><td>1.84</td></tr></tbody></table>

Note: PBE, Perdew-Burke-Ernzerhof. †: Present work.

The plot of the phonon dispersion curves show an excellent agreement between the two functionals (PBE and SCAN) as well as with the experimental curves,<sup>31</sup> see Figure 3. The frequencies obtained with the SCAN functional are systematically higher than those obtained with the PBE functional, about 1 THz, and thus fit better with the experimental data. This can be explained by the fact that the lattice parameters calculated with the SCAN
functional are slightly smaller than those calculated using PBE. Results thus show that the SCAN functional provides a slight improvement over the PBE functional.

As explained, the dielectric tensor $\epsilon_{\parallel}^{\infty}$ of BeO and Born charge was computed and summarized in Table 5. The dielectric constant tensor has two values, one in the basal plane, $\epsilon_{\perp}^{\infty}$, and one in the c-direction, $\epsilon_{\parallel}^{\infty}$. Similarly, the Born effective charge tensor of the Be atom shows two values: $Z_{\perp}^{\text{Be}}$ and $Z_{\parallel}^{\text{Be}}$. The values for the O atoms are the opposite of those for beryllium. The two experimental studies$^{64,65}$ found similar values for the dielectric constants, and the theoretical studies reproduced these values with a difference of less than 5%. The effective Born charges that can be found in the two theoretical studies, which were also obtained using the DFPT (with ABINIT$^{31}$ and Quantum-Espresso$^{41}$), are not only close to each other but also close to our values, with a difference of less than 1%. Calculated values for the dielectric constants overestimate experimental values by 7%, which is low.

From these interatomic forces, thermodynamic properties of the BeO system (Helmholtz vibration energy and its derivative, that is, $F$, $C_u$, etc.) were computed. They were compared to the JANAF thermodynamic dataset.$^{46}$ Results (not shown here) show that a $2 \times 2 \times 2$ supercell is sufficient to accurately generate thermodynamic functions; heat capacity deviation is less than 3 J/mol/K at 1000 K. Results also show that the difference between SCAN and PBE functionals is negligible when calculating thermodynamic functions. There is no significant improvement between SCAN and PBE results, contrary to what was found for the formation energies. Therefore, in what follows, the thermodynamic functions will not be calculated with the SCAN functional. In addition, Born charges and dielectric constants were included in the calculations when considering charged or uncharged H atoms to determine the stability of each insertion site with the phonon dispersion curve, but they were ignored for the transition states.

# 4 | HYDROGEN SOLUBILITY

## 4.1 | Presentation of sites

In the wurtzite system, as shown by Wahl et al.,$^{66}$ there are 10 potential insertion sites, see Figure 1. There are two non-equivalent octahedral sites ($O_{be}$ and $O_o$) located at the center of the beryllium and oxygen octahedrons, respectively. There are also two tetrahedral sites ($T_{be}$ and $T_o$) located in the center of the beryllium and oxygen tetrahedrons, respectively. The R sites are equidistant from the two first-nearest octahedral (1NN) sites and the E sites from the two 1NN tetrahedral sites. In the BeO lattice, Be and O atoms form two types of bonds: an in-plane bond and a

<table><caption>TABLE 6 Wyckoff positions of the insertion sites tested: 6c is in (x, 2x, z), 2b (1/3, 2/3, z), and 2a (0, 0, z).</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">Position</th><th rowspan="2">Wyckoff's label</th></tr><tr><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td>$AB_{be}$</td><td>1/6</td><td>1/3</td><td>0.063</td><td>6c</td></tr><tr><td>$AB_o$</td><td>1/6</td><td>1/3</td><td>0.32</td><td>6c</td></tr><tr><td>$BC_b$</td><td>1/2</td><td>-1/2</td><td>0.44</td><td>6c</td></tr><tr><td>$BC_c$</td><td>2/3</td><td>1/3</td><td>0.689</td><td>2b</td></tr><tr><td>$E$</td><td>2/3</td><td>1/3</td><td>0.189</td><td>2b</td></tr><tr><td>$T_{be}$</td><td>2/3</td><td>1/3</td><td>0.125</td><td>2b</td></tr><tr><td>$T_o$</td><td>2/3</td><td>1/3</td><td>0.253</td><td>2b</td></tr><tr><td>$R$</td><td>0</td><td>0</td><td>0.685</td><td>2a</td></tr><tr><td>$O_{be}$</td><td>0</td><td>0</td><td>3/4</td><td>2a</td></tr><tr><td>$O_o$</td><td>0</td><td>0</td><td>0.635</td><td>2a</td></tr><tr><td>$X$</td><td>2/5</td><td>3/5</td><td>0.6</td><td>6c</td></tr></tbody></table>

perpendicular bond. At the center of each bond, there is an insertion site labeled $BC_c$ for the site at the center of the c-axis perpendicular bond, and $BC_b$ for the site at the center of the in-plane bond. For each in-plane bond, there are two anti-bond insertion sites, $AB_o$ and $AB_{be}$. Finally, an additional position must be considered: the X site (identified in the migration study of $H^+$, see below). For each site, Table 6 indicates the atomic positions in the primitive cell (space group 186) as well as the Wyckoff position.

## 4.2 | Solubility energies

Allouche et al.,$^{14}$ Marinopoulos et al.,$^{16}$ Hodille et al.$^{15}$, and Tsunoda et al.$^{17}$ already published results on hydrogen insertion in the BeO system. Allouche studied only neutral hydrogen in its ground state, the $O_o$ site. Marinopoulos' work does not include the values of hydrogen formation energies. Therefore, only the stability order can be discussed as well as the energy difference between sites. He studied hydrogen in both forms, that is, neutral and charged forms, $\pm 1$. For neutral hydrogen, the ground state is $O_o$, followed by local minima, $BC_c$, $BC_b$, and E. For $H^+$, the most stable site is $BC_c$ and it is locally stable in $AB_o$, $BC_b$ and $T_{be}$-like sites. $H^-$ is found stable in the $O_o$ site only. In addition, Hodille extended the study to molecular hydrogen. His results regarding the ground states for atomic insertion are in full agreement with previous works.$^{14,16}$ For neutral molecular hydrogen, the ground state is also the $O_o$ site with the molecule oriented along the c-axis. Finally, based on another approach seeking the interstitial configurations of hydrogen atoms in oxide, especially efficient for protons, Tsunoda$^{17}$ confirmed the ground states of each hydrogen state. The formation energy of neutral hydrogen is slightly higher (+0.52 eV) in Allouche's study
<table><caption>TABLE 7 Formation energies ($E_f$, in eV) of stable hydrogen insertion sites in BeO for each functional.</caption><thead><tr><th></th><th>Site</th><th>PBE</th><th>SCAN</th><th>ZPE</th><th>μ<sub>m</sub></th><th>BC</th><th>Ω<sub>f</sub></th></tr></thead><tbody><tr><td rowspan="3">H<sup>0</sup></td><td>O<sub>o</sub></td><td>3.99</td><td>4.10</td><td>-0.01</td><td>1</td><td>1.1</td><td>3.5</td></tr><tr><td>BC<sub>c</sub></td><td>4.60</td><td>4.79</td><td>-0.03</td><td>1</td><td>1.2</td><td>7.3</td></tr><tr><td>BC<sub>b</sub></td><td>4.77</td><td>4.99</td><td>-0.03</td><td>1</td><td>1.2</td><td>7.3</td></tr><tr><td rowspan="3">H<sup>+</sup></td><td>BC<sub>c</sub></td><td>3.18</td><td>3.34</td><td>+0.10</td><td>0</td><td>0.5</td><td>-2.3</td></tr><tr><td>AB<sub>o</sub></td><td>3.66</td><td>3.78</td><td>+0.11</td><td>0</td><td>0.4</td><td>3.2</td></tr><tr><td>X</td><td>3.75</td><td>3.92</td><td>0.12</td><td>0</td><td>0.5</td><td>-2.4</td></tr><tr><td>H<sup>-</sup></td><td>O<sub>o</sub></td><td>3.18</td><td>3.34</td><td>+0.09</td><td>0</td><td>1.9</td><td>13.5</td></tr><tr><td>H<sub>2</sub></td><td>O<sub>o</sub></td><td>4.21</td><td>4.32</td><td>+0.23</td><td>0</td><td>0.93/1.2</td><td>6.9</td></tr></tbody></table>

The formation energy of the charged species is given at the Fermi level, set by $E_f(H, +1) = E_f(H, -1)$ in ground states (see Equation (1)). Calculations were done with 3×3×3 supercells. The ZPE is also given (in eV). In addition, the total magnetization ($\mu_m$ in $\mu_B$/cell), the volume of relaxation ($\Omega_f$ in Å³), and the Bader charge (BC per H atom) are also given. For charged species, the relaxation volume is given by Aneto code²⁵ otherwise it is calculated as the volume variation between structures with inserted hydrogen and the ideal lattice. PBE, Perdew–Burke–Ernzerhof.

than in Hodille's. This discrepancy could be explained by the fact that Allouche's work includes a semi-empirical long-range correction (contracting the lattice) with no elastic correction which should have been taken into account because the lattice parameters are set to the ideal lattice. In addition, for neutral hydrogen, Hodille investigated two other local minima (relaxed positions are not described) for which the corresponding stability energies (from ground state) are similar to the PBE results in Marinaupoulos' work. In this latter work, the energy differences rise by approximately 0.25 eV when calculated with the HSE06 hybrid functional.

From these data, Hodille calculated a solubility of 5 × 10⁻¹² appm for a temperature of 773 K and a H₂ pressure of 13.3 kPa. This is much lower than the experimental values, which are above 100 appm.¹³

We have therefore reviewed the content of these studies. Each position presented in the previous section was systematically tested to make sure no cases were missing. Vibrational calculations were always performed to check the stability of sites. Table 7 summarizes the formation energies of stable hydrogen sites obtained in this work. They are given as a function of the chemical state for both PBE and SCAN functionals, taking into account the ZPE.

Results show similarities with previous literature, but some discrepancies are also highlighted. Our results are in full agreement with former studies regarding the most stable sites: the O_o site for H⁰, H⁻ and H₂ species, and the BC_c site for H⁺. Nevertheless, there are differences in the local minima for H⁰ and H⁺. Each form has three stable sites: O_o, BC_c, BC_b, and BC_c, AB_o and the X sites, respectively. For all cases, the stable configurations identified are indeed local minima (insofar as these calculations were made using a non-relaxed volume). The phonons of the different cases were indeed calculated and confirm this result. On the other hand, the other configurations are all unstable. Our phonon spectra show that the two highest sites for H⁺ found by Marinaupoulos are in fact transition states (see Section 5.1.3). The small discrepancies between the formation energies previously published¹⁵ and our own data obtained with the PBE functional could come from the methodology, specifically the elastic and electrostatic correction terms. For molecular hydrogen, without ZPE, values are in perfect agreement, whereas for H⁰, our computations lead to a slightly higher formation energy (approximately 160 meV).

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Charge density difference, $\Delta\rho$, for the investigated hydrogen states in their most stable site along a lattice plane (only atoms in the plane are drawn). From the upper left to lower right picture: H⁰, H₂, H⁺, and H⁻. Color map starts from -10⁻⁴ (blue) to 10⁻⁴ (red) $e⁻/Å³$, whereas green corresponds to 0. Figures were built with VESTA.³⁰</figcaption></figure>

Tests were made to evaluate the influence of the size of the supercell on the formation energy. Calculations were performed on supercells from 2 × 2 × 2 to 4 × 4 × 4. Results show that whichever the species, the formation energy of its ground state increases slightly when the size of the supercell increases from 3 × 3 × 3 and 4 × 4 × 4. For charged, neutral and molecular hydrogen, the energy difference is about 4, 5, and 8 meV, respectively.

In order to analyze the interactions between H and the BeO network, we analyzed the transfer of charges using the charge difference $\Delta\rho$ ($\Delta\rho = \rho[\text{BeO} + \text{H}] - \rho[\text{BeO}] - \rho[\text{H}]$, calculated with the same grid and simulation box) for the different forms of hydrogen. $\Delta\rho$ are drawn in Figure 4. The Bader charge⁶⁷ and the volume of formation, $\Omega_f$, were also
calculated, see Table 7. The volume of formation, Ω<sub>f</sub> = V[BeO + H] - V[BeO] which represents the overall effect of H insertion on the lattice, and the Bader charge both give information on whether or not hydrogen is truly charged.

In the case of H⁰ and H₂, the interactions with O or Be atoms are weak, see Figure 4 (top). The charge transfer between the BeO system and hydrogen upon insertion is weak indeed. The values of the Bader charges confirm that the charge of the inserted hydrogen corresponds to the state considered: for instance, there is one electron around each H atom in H⁰ and H₂. In the case of H⁺, Δρ show a distinguished charge transfer between H and O atoms, leading to the formation of an O–H bond, confirmed by a bond length of 0.97 Å. In the case of H⁻, it can be noted that hydrogen tends to form an electrostatic bond with the first-nearest beryllium atom. Bader charges agree with these results: for H⁺, hydrogen almost recovers its electron and for H⁻, it gains one electron. Positive charges are attracted by O atoms, whereas negative charges by Be atoms to make bonds/interactions.

Regarding the effect of hydrogen insertion on the bulk, results show that, except for H⁺, insertion causes an increase in the lattice size. In the case of neutral hydrogen, for the same site, O₀, Ωₜ is almost twice as large for H₂ as for H⁰. In the case of H⁻, one additional electron is inserted. This expands the volume of the cell, whereas in the case of H⁺ a contraction of the cell can be noticed. However, an exception can be noted for AB₀: H⁺ insertion contracts the system due to the particular position of AB₀ in the supercell.

It should be noted that our calculations, which take into account the magnetic moment, show a localized electron for neutral hydrogen in the vicinity of the H atom. This results in a magnetic moment, μ<sub>B</sub>.

To calculate the solubility, only the most stable site for each hydrogen form was considered. However, to calculate diffusion coefficients, all stable sites for each form are necessary.

## 4.3 | Solubility

The formation energies of H⁰, H⁺, H⁻, and H₂ were calculated. The relationship between these energies and the concentration of the species will now be established. Two cases are considered: hydrogen, either in its atomic form (H⁰, H⁺, H⁻) or in its molecular form (H₂). The solubility is calculated from the Gibbs free energy, Gᵢ(H⁰). To evaluate the solubility of atomic hydrogen, the following reaction is assumed:

$$ \mathrm{BeO} + \frac{1}{2}\mathrm{H}_2 \rightarrow \mathrm{H}^q \quad \text{in BeO.} \qquad (6) $$

According to this reaction, the Gibbs free energy for the insertion of Hq is given by

$$ G_i(\mathrm{H}^q) = E_f(\mathrm{H}^q) + \Delta F_v - \frac{1}{2}[G_{v,t,r}(\mathrm{H}_2)], \qquad (7) $$

where E<sub>f</sub>(H<sup>q</sup>) is the formation energy without the ZPE (see Table 7 for data at 0 K). The term ΔF<sub>v</sub>, computed with Equation (2), corresponds to the difference in vibrational Helmholtz energy between the supercell with inserted hydrogen and the ideal one, including the ZPE. The last term includes the vibrational, rotational, and translational parts of enthalpy and entropy for the gaseous molecular hydrogen (see details in B). The concentration of hydrogen with a q charge, C[H<sup>q</sup>], is then expressed as:

$$ C[\mathrm{H}^q] = \sqrt{\frac{P}{P^0}} \exp\left(-\frac{G_i(\mathrm{H}^q)}{k_B T}\right) \qquad (8) $$

Similarly, for molecular hydrogen, the reaction and its Gibbs free energy of insertion become:

$$ \mathrm{BeO} + \mathrm{H}_2 \rightarrow \mathrm{H}_2 \quad \text{in BeO} \qquad (9) $$

and

$$ G_i(\mathrm{H}_2) = E_f(\mathrm{H}^q) + \Delta F_v - [G_{v,t,r}(\mathrm{H}_2)]. \qquad (10) $$

The concentration of H₂, is then given by the equation:

$$ C[\mathrm{H}_2] = \frac{P}{P^0} \exp\left(-\frac{G_i(\mathrm{H}_2)}{k_B T}\right) \qquad (11) $$

From Equations (8) and (11), the concentrations of the different species were calculated as a function of temperature for a hydrogen pressure P = P⁰ (1 bar). In first approximation, to compute the H⁺ and H⁻ concentrations, an assumption was made: the system is electrically neutral, meaning that their concentrations are equal. This implies that the insertion Gibbs energies are also equal, which in turn determines the Fermi level position of the BeO system. This constraint was computed at each temperature: at 0 K, ΔE<sub>F</sub> is 4.2 eV, see Figure 2. The different values for the solubility of hydrogen and its isotope (tritium) in BeO are represented in Figure 5, and compared to experimental data.⁵,¹³,⁶⁸ Whichever the chemical state of hydrogen isotope, the calculated solubility is always much lower, by several orders of magnitude, than the experimental data. Present results show that, in beryllium oxide, hydrogen should be charged, H⁺/⁻. The concentrations of neutral and molecular hydrogen are found lower than that of charged hydrogen, regardless of the approximations. Furthermore, a slight concentra-
<figure><img src="image_7.png" /><figcaption>FIGURE 5 (Color online) Solubility of hydrogen (atom or molecule) as a function of its charge state. The constraint $[H^+] = [H^-]$ enforces the Fermi level for the charged species. Concentration of charged hydrogen atom for two Fermi levels are also shown. The experimental data related to poly-crystal sample come from Ref. [13], whereas the upper experimental solubility of the single crystal is invoked in Ref. [5]. To show the isotopic effect, only species may have the highest impact is plotted, that is, $T_2$.</figcaption></figure>

tion increase can be observed in the case of isotopes, without quantitatively changing the results. This disagrees with the previous conclusion made by Hodille et al.¹⁵ stating that hydrogen is mainly dissolved as molecular hydrogen. Unfortunately, corresponding experimental results are scarce; there are no available measures performed on single crystals, which leads to an ineluctable bias. From a deuterium-loaded sintered sample, Macaulay-Newcombe et al.⁶⁸ proposed the following Sieverts law: $S = 10^{18} \exp(+0.8/(K_B T))\text{at.m}^{-3}\text{Pa}^{-0.5}$. Nevertheless, this law should be considered as an estimate, at least because it derives from simulations of thermal desorption spectra which include a postulated diffusion coefficient. There is another uncertainty due to a possible misprint regarding the sign of the solution energy in their article. Indeed, the sign is negative in subsequent publications.⁶⁹,⁷⁰ Earlier, these authors performed measurements on oxidized beryllium samples¹³ and obtained a solubility of about 200 appm at 773 K and 13.3 kPa (roughly in agreement with the Sieverts' law above). In this experiment, the oxide layer (approximately 276 nm) was grown with an oxidation stage at 1073 K. Nonetheless, more recent oxidation experiments on beryllium samples concluded to a rapid oxide layer growth with many cracks if the temperature is higher than 950 K.⁷¹ These defects must affect the apparent solubility measured by Macaulay-Newcombe¹³ and explain the huge gap with our data. Finally, Alexander et al.⁵ tried to measure the solubility of deuterium in a single crystal, but the concentration was below the detection limit of their mass spectrometer. This could establish an upper limit for

hydrogen solubility: $[H] \le 4.1 \times 10^{-6}$ atom/f.u. This limit is still significantly distant from our data.

To conclude, it is rather difficult to confirm whether the gaps result from inadequacies coming from the DFT simulations and the solubility model, or if they can be explained by experimental errors (or both). As explained, concentrations were computed based on the following assumption: the charge neutrality was controlled by the content of $H^+$ and $H^-$ species. However, this neutrality can also be achieved through other species or defects. For instance, charged intrinsic defects (vacancy, anti-sites, Schottky defects, etc.) or substitutional species could also modify the position of the Fermi level, see Appendix C. But a change in the position of the Fermi level could greatly affect the concentration of the various charged species. Figure 2 shows the electronic band structure of BeO with different values of the Fermi level: in the middle of the band gap at 4.2 eV when the concentrations of $H^+$ and $H^-$ are equal at 0 K, but also a level at 1.5 eV close to the valence band and at 7eV near the (DFT) conduction band. It should be noted that the concentration of the neutral species remains unchanged in all cases.

The new concentration of $H^+$ and $H^-$ is then represented for these different values of the Fermi level, Figure 5. Concentration values consistent with the experimental data were thus obtained for a Fermi level of 1.5 eV for $H^+$ or 7 eV for $H^-$. The impact of the Fermi level shift on the concentrations can be observed clearly. It is therefore important to be able to determine how the position of the Fermi level in beryllium oxide can be modified, in order to understand how this affects the solubility of hydrogen and its isotopes in beryllium oxide. At this stage, the conclusion is that the position of the Fermi level, influenced by intrinsic defects but also by impurities or doping of the *n* or *p* type, can change the charge state of hydrogen in the BeO system.

## 5 DIFFUSIVITY

Despite the fact that $H^+/-$ should be the preponderant species, it is important to compare the diffusion coefficients of the other forms of hydrogen in order to check if their mobility (product of concentration by diffusion) is indeed much lower than that of $H^+$ or $H^-$. The results on the diffusion of hydrogen and its isotopes in the BeO system are presented below.

### 5.1 Atomistic process

To calculate diffusion coefficients, it is necessary to first describe the elementary processes of hydrogen atom migration at the atomic scale. From the stable sites
<figure><img src="image_8.png" /><figcaption>FIGURE 6 (Color online) Energy profile computed with the PBE functional and diagram of the migration path of H⁻ between O₀ sites || and ⊥ to the c-axis in BeO. Hydrogen paths are drawn in gold whereas, red (green) balls correspond to oxygen (beryllium) atoms. PBE, Perdew-Burke-Ernzerhof.</figcaption></figure>

obtained previously for each state of hydrogen, migration paths were identified. NEB simulations²⁷ were then used to calculate activation energies and find transition states. The statuses of the latter were always confirmed by a phonon analysis. Only configurations with one imaginary frequency were kept. At this point, a new assumption was made: the electrostatic and elastic correction energies (calculations were always performed at fixed lattice parameters) are identical in both transition and stable sites, implying that no correction was taken into account for the diffusivity-related energy data. This is supported by the fact that, for a given hydrogen chemical state, these terms are similar between stable sites. For example, the highest change in energy was found for H⁺ between X and AB₀: about 2 and 25 meV for elastic and electrostatic correction terms, respectively.

### 5.1.1 | Case of H⁻ and H₂

H⁻ and H₂ species are only stable in the O₀ sites. Therefore, diffusion paths are narrowed to two jumps toward the two first-nearest neighboring O₀ sites, see Figures 6 and 7, for H⁻ and H₂, respectively. For the pathway parallel to the c-axis, both species moves along a straight line. Moreover, as the H₂ molecule is initially aligned in that same direction, there is no modification of its orientation. The transition state is located at the intersection of the oxygen layer, whichever the species. In contrast, for the pathway in the x, y-plane, the molecule tilts slightly but, at the transition state, the molecule recovers an orientation mainly along the c-axis, inducing a significant distortion of the first nearest atoms. The transition states for both species are located in the middle of the O₀-O₀ path.

As expected, computed activation energies gathered in Table 8 show that diffusion is anisotropic. The effect of the functional was investigated on both species. The ener-

<figure><img src="image_9.png" /><figcaption>FIGURE 7 (Color online) Energy profile computed with the PBE functional and diagram of the migration path of H₂ between O₀ sites || and ⊥ to the c-axis in the BeO system. Molecular hydrogen paths are drawn in gold and red (green) balls correspond to oxygen (beryllium) atoms. PBE, Perdew-Burke-Ernzerhof.</figcaption></figure>

<table><caption>TABLE 8 Activation energies, Eₐ, related to H⁻ and H₂ diffusion paths in BeO (in eV), PBE/SCAN results.</caption><thead><tr><th rowspan="2"></th><th colspan="2">PBE</th><th colspan="2">SCAN</th></tr><tr><th>H<sup>-</sup></th><th>H<sub>2</sub></th><th>H<sup>-</sup></th><th>H<sub>2</sub></th></tr></thead><tbody><tr><td>O<sub>o</sub> → O<sub>o</sub> ||c</td><td>1.58</td><td>1.86</td><td>1.67</td><td>1.91</td></tr><tr><td>O<sub>o</sub> → O<sub>o</sub> ⊥c</td><td>1.11</td><td>1.16</td><td>1.28</td><td>1.22</td></tr></tbody></table>

PBE, Perdew-Burke-Ernzerhof; SCAN, strongly constrained and appropriately normed.

gies obtained using the SCAN functional are a bit higher than those computed using the PBE, but the difference is not significant enough to affect the diffusion coefficient, see Figure A.2. Whichever the functional, the difference in activation energy between paths parallel and perpendicular to the c-axis is higher than the previously published theoretical data.¹⁵ The main reason is that our saddle point of the path in the x, y-plane is lower than the one found in previous works, but the computed phonon frequencies show an imaginary vibrational mode along the migration path which tends to corroborate that the present result corresponds to the most favorable path in the potential energy surface.

### 5.1.2 | Case of H⁰

Neutral hydrogen is stable in three insertion sites, that is, O₀, BC_c, and BC_b. All available combinations of migration paths are narrowed down to three paths, see Figure 8.

The diffusion mechanism along the c-axis is driven by a jump between two first nearest neighboring O₀ sites, whereas in the ab-plane this path goes through intermediate sites, BC_c or BC_b. To quit an O₀ site, there are two non-equivalent paths toward BC_b. However, the direct diffusion toward the one above BC_b is impossible because hydrogen would go through an upper O₀ site. For all other
<figure><img src="image_10.png" /><figcaption>**FIGURE 8** (Color online) Energy profile computed with the PBE functional and diagram of the migration path of H⁰ in BeO. Hydrogen paths are drawn in gold and red (green) balls correspond to oxygen (beryllium) atoms. PBE, Perdew-Burke-Ernzerhof.</figcaption></figure>

<table><caption>**TABLE 9** Activation energies, *E*<sub>a</sub> (in eV), for the diffusion of neutral hydrogen computed with PBE.</caption><thead><tr><th></th><th>O<sub>o</sub></th><th>BC<sub>c</sub></th><th>BC<sub>b</sub></th></tr></thead><tbody><tr><th>O<sub>o</sub></th><td>1.24</td><td>0.65</td><td>0.77</td></tr><tr><th>BC<sub>c</sub></th><td>0.03</td><td>-</td><td>-</td></tr><tr><th>BC<sub>b</sub></th><td>0.01</td><td>-</td><td>-</td></tr></tbody></table>

*Note*: PBE, Perdew-Burke-Ernzerhof.

possibilities, paths are narrowed to elementary processes including a site, even in the case of very close BC<sub>b</sub> sites surrounding an oxygen atom.

The energy landscape obtained shows that local minima, BC<sub>b</sub> and BC<sub>c</sub>, are near their transition states from a spatial and energy standpoint. This induces a very low energy barrier (10 and 30 meV, respectively) to exit toward the next O<sub>o</sub> site, see Figure 8. By taking into account the ZPE of both sites (stable and saddle), these activation energies still decrease down to a null value (≤ 0 eV) considering hydrogen diffusion along the BC<sub>b</sub>→O<sub>o</sub> path. Nonetheless, temperature acts oppositely, therefore the activation Helmholtz energies increase the energy barrier of both migration paths. The activation energies of the three migration paths are given in Table 9.

### 5.1.3 | Case of H⁺

H⁺ is stable in three insertion sites: BC<sub>c</sub>, X, and AB<sub>o</sub>. As shown in Section 4.2, a charge transfer between hydrogen and a first nearest-neighboring oxygen atom induces a complex energy landscape. All the possible diffusion paths are narrowed down to five elementary processes, see Figure 9. The migration between two local minima AB<sub>o</sub> is only possible with in-plane sites, where the jump between two sites surrounding an oxygen atom passes through the intermediate local minimum X. The saddle point of this

<figure><img src="image_11.png" /><figcaption>**FIGURE 9** (Color online) Energy profile computed with the PBE functional and diagram of the migration path of H⁺ between stable sites. The identified hydrogen paths are drawn in gold and red (green) balls correspond to oxygen (beryllium) atoms. PBE, Perdew-Burke-Ernzerhof.</figcaption></figure>

<table><caption>**TABLE 10** Activation energies of H⁺ diffusion paths computed with PBE functional, *E*<sub>a</sub> (in eV).</caption><thead><tr><th></th><th>AB<sub>o</sub></th><th>BC<sub>c</sub><sup>1</sup></th><th>BC<sub>c</sub><sup>2</sup></th><th>X</th></tr></thead><tbody><tr><th>AB<sub>o</sub></th><td>0.23</td><td>-</td><td>0.02</td><td>0.13</td></tr><tr><th>BC<sub>c</sub><sup>1</sup></th><td>-</td><td>-</td><td>-</td><td>0.77</td></tr><tr><th>BC<sub>c</sub><sup>2</sup></th><td>0.47</td><td>-</td><td>-</td><td>-</td></tr><tr><th>X</th><td>0.006</td><td>0.18</td><td>-</td><td>0.078</td></tr></tbody></table>

*Note*: PBE, Perdew-Burke-Ernzerhof.

latter path is a second-order transition state according to the multiplicity of the AB<sub>o</sub> site and its location in the middle of the Be-O in-plane bond (BC<sub>b</sub>-like site, see Figure 1). To escape from a BC<sub>c</sub> site, there are only two out-of-plane paths: toward AB<sub>o</sub> and X. Additionally, the latter site has a three-fold axial symmetry along the c-axis meaning that three X sites are surrounding their saddle point (E site, a second-order transition state). The exit paths out of X sites have low activation energies, especially toward AB<sub>o</sub> (6 meV) which is in line with a flat potential energy. Indeed, the two computed Eigen-frequencies related to hydrogen are low: one is complex, whereas the other corresponds to a soft mode close to a few THz. This latter should be found complex if computational limitations were to be overcome (sampling of Brillouin zone beyond Γ). Nevertheless, when taking into account ZPE contributions, this activation energy increases by several meV: from +10 up to +21 meV for hydrogen and tritium, respectively. Thus, as temperature increases, the energy barrier between X and AB<sub>o</sub> sites also increases. In contrast, the second-lowest activation energy, which corresponds to the AB<sub>o</sub> → BC<sub>c</sub><sup>2</sup> path, decreases if ZPE is taken into account (−4 meV and −2.8 meV for H and T, respectively). The other activation energies of migration paths are summarized in Table 10.
## 5.2 | Diffusion coefficients of H and its isotopes

Two approaches were used to compute the diffusion coefficients: explicit expressions and kinetic Monte-Carlo (KMC) simulations. These methods are described thoroughly in Appendix A. Both approaches use the harmonic transition state theory according to the Eyring's theory<sup>72</sup> to compute the jump rate probabilities. These jump rates reflect the probability of an atom leaving site a toward site b, it is expressed as:

$$ \Gamma_{ab}(T) = \frac{k_B T}{h} \frac{Z_{TS}}{Z_a} \exp \left[ -\frac{\Delta E_m^{ab}}{k_B T} \right]. \quad (12) $$

where $\Delta E_m^{ab}$ is the activation energy (difference between the saddle point and the initial local minima), $Z_{TS}$ and $Z_a$ correspond to the (vibration) partition function of the transition state and the initial position, respectively. $Z$ and the Helmholtz free energy are linked through:

$$ F_v(T) = -k_B T \ln Z. \quad (13) $$

The migration mechanism of $H^-$ and $H_2$ is limited to jump between $O_o$ sites, for which only one site is stable. The diffusion coefficients in the ab-plane, $D_{x,y}$, and along the c-axis, $D_z$, can then be expressed as:

$$ D_{x,y}(T) = \frac{3a^2}{2} \Gamma_{OOxy} \quad (14) $$

$$ D_z(T) = \frac{c^2}{4} \Gamma_{OOz}, \quad (15) $$

where $\Gamma_{OOxy}$ corresponds to the jump probability in the ab-plane, $\Gamma_{OOz}$ is the jump rate along $z$, and $a$ and $c$ are the lattice parameters. For $H^0$, diffusion coefficients are given by:

$$ D_{x,y}(T) = \frac{a^2}{4} \frac{\Gamma_{co} \Gamma_{bo} (2\Gamma_{oc} + 3\Gamma_{ob})}{\Gamma_{bo} \Gamma_{co} + \Gamma_{bo} \Gamma_{oc} + 3\Gamma_{co} \Gamma_{ob}} \quad (16) $$

$$ D_z(T) = \frac{c^2}{4} \frac{\Gamma_{bo} \Gamma_{co} \Gamma_{oo}}{\Gamma_{bo} \Gamma_{co} + 3\Gamma_{co} \Gamma_{ob} + \Gamma_{bo} \Gamma_{oc}}, \quad (17) $$

where expressions of jump rates are simplified as: $BC_c \to c$, $BC_b \to b$ and $O_o \to o$. For $H^+$, an explicit expression cannot be obtained due to the complex migration pathways which involve high multiplicity sites ($AB_o$ and $X$). Therefore, it overall diffusion coefficient is evaluated by KMC.

The diffusion coefficients of species $H^0$, $H^-$, and $H_2$ are found highly anisotropic, see Figure 10; species move faster in the ab-plane, while diffusion along the c-axis is several orders of magnitude lower and this is also true for higher

<figure><img src="image_12.png" /><figcaption>FIGURE 10 (Color online) Diffusion coefficients of $T^+$, $T^-$, $T_2$, and $T$ in the BeO bulk compared to the experimental data of tritium diffusion in single BeO crystal.<sup>12</sup></figcaption></figure>

temperatures. In contrast, the diffusion of $H^+$ exhibits an opposite behavior. The trajectories obtained with KMC show a slightly anisotropic behavior, with a displacement along the c-axis that is longer than the displacements in the ab-plane. This behavior increases significantly as the temperature rises, whereas at low temperature, diffusion paths tend to be less anisotropic, see Figure A.1. Results show that diffusion coefficients extend over a large range of values depending on the state of hydrogen. Unexpectedly, the diffusivity of hydrogen in beryllium oxide (tritium species, Figure 10) can be the same order of magnitude as in the metal, especially for neutral hydrogen. Here, the theoretical diffusion coefficient of neutral hydrogen is even higher than in beryllium metal for temperatures above 1 000 K.<sup>11</sup> Thus, the ability of beryllium oxide to act as a permeation barrier becomes questionable due to the interstitial diffusion mechanism found here.

Experimental data on the diffusion of hydrogen isotopes in beryllium oxide are not numerous, but data on tritium diffusion in BeO single crystal are available.<sup>12</sup> Fowler et al.<sup>12</sup> used the recoil injection method to measure the diffusion coefficient of tritium. Single crystal samples were first covered with a lithium salt layer, which was then exposed to a neutron flux to produce tritium. To evaluate the diffusion coefficient of tritium, the time rate of tritium release was measured in post-irradiated conditions, where samples were heated (tritium release was measured with a mass spectrometer). The authors did not specify the crystallographic orientation,<sup>5,12</sup> the data should therefore be considered as an effective diffusion coefficient. In addition, the authors claim that the sample release is apparently not controlled by defects, but rather by a diffusion mechanism. Unfortunately, the chemical state of the diffusing species should be considered as unknown. Even
<table><caption>TABLE 11 Activation energies, *E*<sub>a</sub> (in eV) and *D*<sub>0</sub> (in m<sup>2</sup>/s) obtained from an Arrhenius fit for hydrogen and tritium diffusing in BeO.</caption><thead><tr><th rowspan="2">Species</th><th colspan="2"><i>E</i><sub>a</sub></th><th colspan="2"><i>D</i><sub>0</sub></th></tr><tr><th><i>x</i>, <i>y</i></th><th><i>z</i></th><th><i>x</i>, <i>y</i></th><th><i>z</i></th></tr></thead><tbody><tr><td>H</td><td>0.69</td><td>1.25</td><td>1.38×10<sup>-5</sup></td><td>2.75×10<sup>-6</sup></td></tr><tr><td>H<sup>-</sup></td><td>1.07</td><td>1.56</td><td>2.36×10<sup>-7</sup></td><td>6.93×10<sup>-6</sup></td></tr><tr><td>H<sub>2</sub></td><td>1.11</td><td>1.83</td><td>1.64×10<sup>-5</sup></td><td>1.86×10<sup>-5</sup></td></tr><tr><td>H<sup>+</sup></td><td>0.63</td><td></td><td>2.36 × 10<sup>-7</sup></td><td></td></tr><tr><td>T</td><td>0.68</td><td>1.25</td><td>9.91×10<sup>-6</sup></td><td>1.67×10<sup>-6</sup></td></tr><tr><td>T<sup>-</sup></td><td>1.08</td><td>1.57</td><td>1.80×10<sup>-7</sup></td><td>4.94×10<sup>-6</sup></td></tr><tr><td>T<sub>2</sub></td><td>1.13</td><td>1.84</td><td>1.02×10<sup>-5</sup></td><td>1.14×10<sup>-5</sup></td></tr><tr><td>T<sup>+</sup></td><td>0.68</td><td></td><td>2.69 × 10<sup>-7</sup></td><td></td></tr></tbody></table>

if the neutron transmutation produces T<sup>+</sup>, this species can react to form molecular tritium or perhaps a hydroxyl group. The authors pointed out that the formation of molecular tritium can be expected before atoms reach the free surface of the sample, even with a very low tritium concentration. The comparison between computed and experimental diffusion coefficients is rather weak, whichever the species. In addition, the effective activation energies deduced from the Arrhenius law (Table 11) are not in agreement with the experimental interpolation (2.27 eV). These results tend to confirm Fowler's conclusion¹² stating that this high activation energy suggests that species diffusion involved complex defects.

## 6 | CONCLUSION

The goal of this work was to determine the diffusion coefficients of hydrogen in beryllium oxide. To do so, the properties of the ideal lattice as well as the behavior of hydrogen isotopes in interstitial sites were investigated. Our calculations show that the SCAN functional improves on the determination of electronic properties, especially in the case of formation energies, which can be predicted with an error below 1%. This functional also improves significantly on the determination of elastic constants. Nonetheless, the improvement on vibrational properties is rather small compared to PBE results, justifying the use of the latter to compute vibrational Helmholtz energies.

The most stable sites found in this work are the same as those previously identified in the literature, whichever the chemical state of hydrogen, H<sup>0</sup>, H<sup>+</sup>, H<sup>-</sup>, or H<sub>2</sub>. Nonetheless, our results show that some previously identified local minima are in fact unstable or saddle points (sites *E* and *BC*<sub>b</sub> of H<sup>+</sup>). Therefore, all the insertion sites and their formation energies were computed. The electronic structure analysis highlighted a significant charge transfer in the

case of cation insertion, whereas for the other species, the charge remained localized around the hydrogen atom. A simple solubility model was used to perform an estimate of hydrogen concentration in the ideal lattice. Moreover, our predicted hydrogen concentration falls to several orders of magnitude lower than the experimental data. Nonetheless, our calculations predict that hydrogen in beryllium oxide would be mainly present as cation or anion. By changing the constraint related to the Fermi level, the model displays its effect on solubility and shows the possibility to take into account point defects such as beryllium/oxygen vacancies in order to improve solubility determination.

Regarding migration pathways, all paths were carefully studied to identify every transition state and the corresponding activation energies. Due to the charge transfer, anion migration shows a complex diffusion pathway through five elementary processes, two of which pass through a second-order transition state. Based on present reliable DFT data (paths, activation energies, and Helmholtz energies), two approaches were used to compute diffusion coefficients. Indeed, analytical expressions were established for H<sup>0</sup>, H<sup>-</sup>, and H<sub>2</sub> species, whereas KMC was implemented to calculate the diffusion properties of H<sup>+</sup>. Our study shows an insignificant isotopic effect. Furthermore, the diffusion mechanisms exhibit an anisotropic behavior whichever the chemical state. Nevertheless, hydrogen cation diffuses faster, along the *c*-axis, in the opposite direction to the other species. Finally, comparing our own set of diffusion coefficients against experimental data shows a knowledge gap that could be filled by extending both theoretical and experimental investigations.

## ACKNOWLEDGMENTS

This work was performed using HPC resources from CALMIP (Grant 2021-p0749) and GENCI-TGCC (Grants A0080911483 and A0100911483). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Program (Grant Agreement No.: 101052200, EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. The authors would like to thank Marc Barrachin for the valuable discussions and his comments on this work, as well as Laura Ferry for her support regarding KMC tool.

## ORCID

Damien Connétable ID https://orcid.org/0000-0003-3642-780X
REFERENCES

1. Petrovic J. Los Alamos National Laboratory, editor. Beryllium oxide (BeO) handbook; 2020. LA-UR-20-24561.

2. Walsh KA. Beryllium chemistry and processing. ASM International; 2009.

3. Li K, Qian L, Li XMY, Zhou W. BeO utilization in reactors for the improvement of extreme reactor environments—a review. Front Energy Res. 2021;9:669832.

4. NEA. State-of-the-art report on light water reactor accident tolerant fuels. Paris: OECD; 2018.

5. Alexander C, Causey R, Chari D, Elleman T, Feng P, Fowler J. Tritium diffusion in nonmetallic solids of interest for fusion reactors. North Carolina State University; 1979.

6. Federici G, Doerner RP, Lorenzetto P, Barabash V. Beryllium as a plasma facing material for near term fusion devices. Elsevier; 2012.

7. Houssa M, Pourtois G, Heyns MM, Stesmans A. Defect generation in high *k* gate dielectric stacks under electrical stress: the impact of hydrogen. J Phys: Condens Matter. 2005;17:2075.

8. Shimada M, Pitts RA, Ciattaglia S, Carpentier S, Choi CH, Orco GD, et al. In-vessel dust and tritium control strategy in ITER. J Nucl Mater. 2013;438:S996.

9. Ganchenkova MG, Borodin VA, Nieminen RM. Hydrogen in beryllium: solubility, transport, and trapping. Phys Rev B. 2009;79:134101.

10. Matveev D, Wensing M, Ferry L, Virot F, Barrachin M, Ferro Y, et al. Reaction-diffusion modeling of hydrogen transport and surface effects in application to single-crystalline Be. Nucl Instrum Methods Phys Res, Sect B. 2018;430:23.

11. Ferry L, Virot F, Ferro Y, Matveev D, Linsmeier C, Barrachin M. Diffusivity of hydrogen and properties of point defects in beryllium investigated by DFT. J Nucl Mater. 2019;524: 323.

12. Fowler J, Chandra, Elleman T, Payne A, Verghese K. Tritium diffusion in Al₂O₃ and BeO. J Am Ceram Soc. 1977;60:155.

13. Macaulay-Newcombe RG, Thompson DA, Smeltzer WW. Thermal absorption and desorption of deuterium in beryllium and beryllium oxide. J Nucl Mater. 1992;191-194:263.

14. Allouche A, Ferro Y. First-principles study of hydrogen retention and diffusion in beryllium oxide. Solid State Ionics. 2015; 272:91.

15. Hodille EA, Ferro Y, Piazza ZA, Pardanaud C. Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen. J Phys: Condens Matter. 2018;30:305201.

16. Marinopoulos AG, Vilão RC, Vieira RBL, Alberto HV, Gil JM, Yakushev MV, et al. Defect levels and hyperfine constants of hydrogen in beryllium oxide from hybrid-functional calculations and muonium spectroscopy. Philos Mag. 2017;97:2108.

17. Tsunoda N, Kumagai Y, Oba F. Recommendation of interstitial hydrogen positions in metal oxides. Comput Mater Sci. 2022;203:111068.

18. Hohenberg P, Kohn W. Inhomogeneous electron gas. Phys Rev. 1964;136:B864.

19. Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effect. Phys Rev. 1965;140:A1133.

20. Kresse G, Hafner J. Ab initio molecular dynamics for liquid metals. Phys Rev B. 1993;47:558.

21. Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Phys Rev Lett. 1996;77:3865.

22. Sun J, Ruzsinszky A, Perdew JP. Strongly constrained and appropriately normed semilocal density functional. Phys Rev Lett. 2015;115:036402.

23. Kresse G, Joubert D. From ultrasoft pseudo-potentials to the projector augmented-wave method. Phys Rev B. 1999;59:1758.

24. Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations. Phys Rev B. 1976;13:5188.

25. Varvenne C, Bruneval F, Marinica MC, Clouet E. Point defect modeling in materials: coupling ab initio and elasticity approaches. Phys Rev B. 2013;88:134102.

26. Freysoldt C, Neugebauer J, Van de Walle CG. Electrostatic interactions between charged defects in supercells. Phys Status Solidi B. 2011;248:1067.

27. Henkelman G, Uberuaga BP, Jónsson H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J Chem Phys. 2000;113(22):9901.

28. Togo A, Oba F, Tanaka I. First-principles calculations of the ferroelastic transition between rutile-type and CaCl₂-type SiO₂ at high pressures. Phys Rev B. 2008;78:134106.

29. Aroyo MI, Perez-Mato JM, Capillas C, Kroumova E, Ivantchev S, Madariaga G, et al. Bilbao Crystallographic server: I. Databases and crystallographic computing programs. Z Kristallogr Cryst Mater. 2009;221:15.

30. Momma K, Izumi F. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. J Appl Crystallogr. 2011;44:1272.

31. Bosak A, Schmalzl K, Krisch M, Van Beek W, Kolobanov V. Lattice dynamics of beryllium oxide: Inelastic x-ray scattering and ab initio calculations. Phys Rev B. 2008;77:224303.

32. Sirota NN, Kuzmina AM, Orlova NS. Debye-Waller factors and elastic constants for beryllium oxide at temperatures between 10 and 720 K (I). Anisotropy of ionic mean-square displacements. Cryst Res Technol. 1992;27:703.

33. Vidal Valat G, Vidal JP, Kurki-Suonio K, Kurki-Suonio R. Multipole analysis of X-ray diffraction data on BeO. Acta Crystallogr A: Found Crystallogr. 1987;43:540.

34. Hazen RM, Finger LW. High-pressure and high-temperature crystal chemistry of beryllium oxide. J Appl Phys. 1986;59:3728.

35. Sabine TM, Hogg S. The wurtzite Z parameter for beryllium oxide and zinc oxide. Acta Crystallogr B: Struct Sci. 1969;25:2254.

36. Jeffrey GA, Parry GS, Mozzi RL. Study of the wurtzite-type binary compounds. I. Structures of aluminum nitride and beryllium oxide. J Chem Phys. 1956;25:1024.

37. Byggmästar J, Hodille EA, Ferro Y, Nordlund K. Analytical bond order potential for simulations of BeO 1D and 2D nanostructures and plasma-surface interactions. J Phys: Condens Matter. 2018;30:135001.

38. Cinthia AJ, Priyanga GS, Rajeswarapalanichamy R, Iyakutti K. Structural, electronic and mechanical properties of alkaline earth metal oxides MO (M=Be, Mg, Ca, Sr, Ba). J Phys Chem Solids. 2015;79:23.

39. Li D, Zhang P, Yan J. Ab initio molecular dynamics study of high-pressure melting of beryllium oxide. Sci Rep. 2014;4:4707.

40. Shein IR, Gorbuniva MA, Makurin YN, Kiiko VS, Ivanovskii AL. Magnetization of beryllium monoxide (BeO) without magnetic impurities: a first-principles study. Int J Mod Phys B. 2008;22:4987.

41. Duman S, Sütlü A, Bağcı S, Tütüncü HM, Srivastava GP. Structural, elastic, electronic, and phonon properties of zinc-blende and wurtzite BeO. J Appl Phys. 2009;105:033719.
42. Yu BR, Yang JW, Guo HZ, Ji GF, Chen XR. Phase transition and elastic properties of BeO under pressure from first-principles calculations. Physica B. 2009;404:1940.

43. Park CJ, Lee SG, Ko YJ, Chang KJ. Theoretical study of the structural phase transformation of BeO under pressure. Phys Rev B. 1999;59:13501.

44. Chang KJ, Cohen ML. Theoretical study of BeO: structural and electronic properties. Solid State Commun. 1984;50:487.

45. Haynes WM, Bruno TJ. Handbook of chemistry and physics. vol. 1. 96th ed. CRC press; 2016.

46. Chase MW. NIST-JANAF Themochemical Tables, 4th ed. J Phys Chem Ref Data. 1998;9:1-1951.

47. Lichanot A, Chaillet M, Larrieu C, Dovesi R, Pisani C. Ab initio Hartree-Fock study of solid beryllium oxide: structure and electronic properties. Chem Phys. 1992;164:383.

48. Sorokin PB, Fedorov AS, Chernozatonskiǐ LA. Structure and properties of BeO nanotubes. Phys Solid State. 2006;48:398.

49. Wrasse EO, Baierle RJ. First principles study of native defects in BeO. Phys Procedia. 2012;28:79.

50. Roessler DM, Walker WC, Loh E. Electronic spectrum of crystalline beryllium oxide. J Phys Chem Solids. 1969;30: 157.

51. McLeod JA, Wilks RG, Skorikov NA, Finkelstein LD, Abu-Samak M, Kurmaev EZ, et al. Band gaps and electronic structure of alkaline-earth and post-transition-metal oxides. Phys Rev B. 2010;81:245123.

52. Emeline AV, Kataeva GV, Ryabchuk VK, Serpone N. Photo-stimulated generation of defects and surface reactions on a series of wide band gap metal-oxide solids. J Phys Chem B. 1999;103(43):9190-9.

53. Milman V, Warren MC. Elasticity of hexagonal BeO. J Phys Condens Mat. 2001;13:241.

54. Baumeier B, Krüger P, Pollmann J. Atomic and electronic structure of BeO and the BeO(1010) surface: an ab initio investigation. Phys Rev B. 2007;75:045323.

55. Groh D, Pandey R, Sahariah MB, Amzallag E, Baraille I, Rérat M. First-principles study of the optical properties of BeO in its ambient and high-pressure phases. J Phys Chem Solids. 2009;70:789.

56. Shahrokhi M, Leonard C. Quasi-particle energies and optical excitations of wurtzite BeO and its nanosheet. J Alloy Compd. 2016;682:254.

57. Sirota NN, Kuzmina AM, Orlova NS. Debye-Waller factors and elastic constants for beryllium oxide at temperatures between 10 and 720 K (II). Elastic constants from X-ray diffraction data. Cryst Res Technol. 1992;27:711.

58. Cline CF, Dunegan HL, Henderson GW. Elastic constants of hexagonal BeO, ZnS, and CdSe. J Appl Phys. 1967;38: 1944.

59. Bentle GG. Elastic constants of single-crystal BeO at room temperature. J Am Ceram Soc. 1966;49:125.

60. Duan Y, Qin L, Tang G, Shi L. First-principles study of ground- and metastable-state properties of XO (X = Be, Mg, Ca, Sr, Ba, Zn and Cd). Eur Phys J B. 2008;66:201.

61. Shein IR, Kiįko VS, Makurin YN, Gorbunova MA, Ivanovskii AL. Elastic parameters of single-crystal and polycrystalline wurtzite-like oxides BeO and ZnO: ab initio calculations. Phys Solid State. 2007;49:1067.

62. Baroni S, Resta R. Ab initio calculation of the macroscopic dielectric constant in silicon. Phys Rev B. 1986;33:7017.

63. GajdošM, Hummer K, Kresse G, Furthmüller J, Bechstedt F. Linear optical properties in the projector-augmented wave methodology. Phys Rev B. 2006;73:045112.

64. Hofmeister AM, Keppel E, Speck AK. Absorption and reflection infrared spectra of MgO and other diatomic compounds. Mon Not R Astron Soc. 2003;345:16.

65. Loh E. Optical phonons in BeO crystals. Phys Rev. 1968;166:673.

66. Wahl U, Vantomme A, Langouche G, et al. Emission channeling studies of Pr in GaN. J Appl Phys. 2000;88:1319.

67. Henkelman G, Arnaldsson A, Jónsson H. A fast and robust algorithm for Bader decomposition of charge density. Comput Mater Sci. 2006;36:354.

68. Macaulay-Newcombe RG, Thompson DA. The effects of surface contamination on absorption and desorption of deuterium in beryllium and beryllium oxide. J Nucl Mater. 1994;212-215:942.

69. Zakharov AP, Gorodetsky AE. Behavior of deuterium implanted into W, Be and BeO. Phys Scr. 2001;T94:68.

70. Billone MC, Donne MD, Macaulay-Newcombe RG. Status of beryllium development for fusion applications. Fusion Eng Des. 1995;27:179-90.

71. Roth J, Doerner R, Baldwin M, others M. Oxidation of beryllium and exposure of beryllium oxide to deuterium plasmas in PISCES B. J Nucl Mater. 2013;438:S1044.

72. Eyring H. The activated complex in chemical reactions. J Chem Phys. 1935;3:107.

73. Landman U, Shlesinger MF. Stochastic theory of multistate diffusion in perfect and defective systems. I. Mathematical formalism. Phys Rev B. 1979;19:6207.

74. Landman U, Shlesinger MF. Stochastic theory of multistate diffusion in perfect and defective systems. II. Case studies. Phys Rev B. 1979;19:6220.

75. Scotti L, Mottura A. Interstitial diffusion of O, N, and C in α-Ti from first-principles: analytical model and kinetic Monte Carlo simulations. J Chem Phys. 2016;144:084701.

**How to cite this article:** Riot A, Virot F, Connétable D. Solubility and diffusivity of hydrogen and its isotopes in the BeO system. J Am Ceram Soc. 2023;106:5005-5021.
https://doi.org/10.1111/jace.19130

## APPENDIX A: MULTI-SCALE APPROACH OF DIFFUSIVITY

Here, we present the method for the calculation of the scattering coefficients with their expressions. Two cases are discussed below: the diffusion of H⁰ and H⁻.

### A.1 | Diffusion of H⁻ and H₂

In the case of H⁻ and H₂, there is one position: O₀ (0, 0, 0.635), labeled 2a. Then, the different possible jumps from each site were identified. As shown previously, there are two different jumps, see Figure 6, and therefore two different jump rates: with Γ<sub>ooXY</sub> (6 equivalent jumps, length 2.71 Å), the jump between two 1NN O₀ sites in the plane, and
$\Gamma_{ooZ}$ along the c-axis, (2 jumps, up and down, length 2.46 Å). By following Landman' method,⁷³,⁷⁴ we find:

$$D_{x,y}[T] = \frac{3a^2}{2} \Gamma_{ooXY}[T] \qquad (A.1)$$

$$D_z[T] = c^2 \Gamma_{ooZ}[T] \qquad (A.2)$$

These formulas were obtained by first identifying the number of non-equivalent positions in the primitive cell. Two quantities are thus required: the Laplace transform of the waiting time density matrix, $\underline{\psi}(u)$, and the Fourier transform matrix of the displacements of hydrogen in BeO, $\underline{p}(k)$. $\underline{\psi}(u)$ is given by:

$$\underline{\psi}(u) = \frac{2a_1}{2a_2} \begin{bmatrix} \frac{6\Gamma_b}{K} & \frac{2\Gamma_c}{K} \\ \frac{2\Gamma_c}{K} & \frac{6\Gamma_b}{K} \end{bmatrix} \qquad (A.3)$$

where $K = 6\Gamma_{ooXY} + 2\Gamma_{ooZ} + u$, and $\underline{p}(k)$ by:

$$\underline{p}(k) = \frac{2a_1}{2a_2} \begin{bmatrix} B & A \\ A & B \end{bmatrix} \qquad (A.4)$$

where

$$B = [e^{ik_x a} + e^{-ik_x a} + e^{ik_y a} + e^{-ik_y a} + e^{ia(k_y+k_x)} + e^{-ia(k_y+k_x)}]/6 \qquad (A.5)$$

$$A = [e^{ik_z c} + e^{-ik_z c}]/2.$$

$D_{x,y}$ and $D_z$ were deduced with $\underline{p}(k)$ and $\underline{\psi}(u)$, and by following the procedure suggested by Landmann.⁷⁴

## A.2 | Diffusion of H⁰

In this case, there are three non-equivalent positions: $O_o$ (labeled 2a), $BC_c$ (2b), and $BC_b$ (6c) sites, see Figure 8. The jump rates are: $\Gamma_{ooZ}$ between two 1NN $O_o$ sites along the c-axis, $\Gamma_{bo}$ and $\Gamma_{ob}$, between $O_o$ and $BC_b$ sites, and $\Gamma_{co}$ and $\Gamma_{oc}$, between $O_o$ and $BC_c$ sites. Using the same approach, $D_{x,y}$ and $D_z$ can be expressed as:

$$D_{x,y}(T) = \frac{a^2}{4} \frac{\Gamma_{co}\Gamma_{bo}(2\Gamma_{oc} + 3\Gamma_{ob})}{\Gamma_{bo}\Gamma_{co} + \Gamma_{bo}\Gamma_{oc} + 3\Gamma_{co}\Gamma_{ob}} \qquad (A.6)$$

$$D_z(T) = \frac{c^2}{4} \frac{\Gamma_{bo}\Gamma_{co}\Gamma_{ooZ}}{\Gamma_{bo}\Gamma_{co} + 3\Gamma_{co}\Gamma_{ob} + \Gamma_{bo}\Gamma_{oc}} \qquad (A.7)$$

(A.8)

## A.3 | Diffusion of H⁺ in the BeO system

For H⁺, an analytical formula could not be established due to complex migration pathways involving high-

<figure><img src="image_13.png" /><figcaption>FIGURE A.1 400 000 steps trajectories of the H⁺ species. The top (bottom) trajectory has been computed at 600 K (1 200 K). For clarity, coordinates have been printed each 1 000 steps.</figcaption></figure>

multiplicity stable sites. Therefore, its diffusion coefficients were investigated using the KMC method. There are many KMC techniques, but this work used the method based on the residence time algorithm, which involves a jump at each time step. The steps of a typical simulation are detailed in Ref. [75]. From the three stable sites, $BC_c$, $AB_o$ and $X$, all possible jumps were defined depending on their multiplicity and paths previously introduced (Section 5.1.3). Each jump was initialized with a displacement vector and a jump probability (Equation (12)), while for each site, a residence time was computed as the reciprocal of the sum of jump rates. To determine the diffusion coefficient, KMC simulations were performed at several temperatures. Then, the recorded trajectories as well as the elapsed time were processed. The mean square displacement ($\langle R_k^2 \rangle$) and the elapsed time ($\langle t_k \rangle$) were evaluated by means of sliding average at each time step with a number of subsets ($k$) going up to half of the step numbers to make sure averages were weighted equally. According to Einstein's equation, the diffusion coefficient is expressed as follows:

$$D_k = \frac{\langle R_k^2 \rangle}{6 \langle t_k \rangle} \qquad (A.9)$$
<figure><img src="image_14.png" /><figcaption>FIGURE A.2 (Color online) Diffusion coefficients of T⁻ and T₂ computed with PBE and SCAN functionals. Experimental data were extracted from Ref. [12]. PBE, Perdew-Burke-Ernzerhof; SCAN, strongly constrained and appropriately normed.</figcaption></figure>

Then, an average diffusion coefficient for a given temperature was obtained. Its evolution as a function of the step was used to check when convergence was reached. In the case of H⁺ migration, simulations were extended to 800 000 steps. Finally, from the Arrhenius law, the effective activation energy and the pre-exponential factor were evaluated by fitting the computed diffusion coefficients.

### A.4 | Example of trajectories computed with KMC for the H⁺ species

These plots (Figure A.1) show H⁺ trajectories computed with the KMC method in order to highlight the slightly anisotropic nature of the diffusion.

### A.5 | Comparison of diffusion coefficients computed from PBE and SCAN data

Figure A.2 shows the impact of the chosen functional (PBE vs. SCAN) on the diffusion coefficients of T⁻ and T₂. The effect is more important for charged hydrogen, by one order of magnitude, specially at low temperature. The PBE functional provides accurate data at this stage of approximation.

## APPENDIX B: VIBRATIONAL, TRANSLATIONAL, ROTATIONAL CONTRIBUTIONS TO GIBBS ENERGY FOR H₂(g)

The vibrational enthalpy in the case of molecules, $H_{vib}$, is expressed as

$$H_{\text{vib}}[T] = \frac{h\nu}{2} + \frac{h\nu}{\exp\left(\frac{h\nu}{k_B T}\right) - 1} \quad (\text{B.1})$$

<figure><img src="image_15.png" /><figcaption>FIGURE C.1 (Color online) Formation energies of interstitial hydrogen in BeO as function of Fermi level, calculated with PBE. Only the most stable site for each hydrogen form was considered. The origin of Fermi energy is set at the top of valence band maximum. PBE, Perdew-Burke-Ernzerhof.</figcaption></figure>

and the vibrational entropy, $S_{vib}$:

$$S_{\text{vib}}[T] = \frac{h\nu}{T} \frac{1}{\exp\left(\frac{h\nu}{k_B T}\right) - 1} - k_B \ln\left(1 - e^{-\frac{h\nu}{k_B T}}\right) \quad (\text{B.2})$$

The translation and rotation contributions to enthalpy and entropy are given by:

$$H_{\text{tran,rot}}[T] = H_{\text{tran}} + H_{\text{rot}} = \frac{5}{2}k_B T + k_B T \quad (\text{B.3})$$

$$ \begin{aligned} S_{\text{tran,rot}}[T] &= S_{\text{tran}} + S_{\text{rot}} \\ &= k_B \left( 1 + \ln \left[ \frac{8\pi^2 I k_B T}{\sigma h^2} \right] \right) \\ &\quad + k_B \left( \frac{5}{2} + \ln \left[ \frac{k_B T}{P^0} \left( \frac{2\pi m k_B T}{h^2} \right)^{\frac{3}{2}} \right] \right) \end{aligned} \quad (\text{B.4}) $$

where m is the mass of the H₂ molecule, I its moment inertia, ν its vibrational frequency, and σ = 2 for diatomic molecules.

## APPENDIX C: VARIATION OF HYDROGEN FORMATION ENERGIES WITH THE FERMI LEVEL

Figure C.1 shows the effect of the Fermi level variation on the formation energies of charged interstitial hydrogen in comparison of the neutral cases. Whatever the Fermi level position, the neutral hydrogen atom and molecular hydrogen are not favorable; hydrogen in beryllium oxide exhibits an amphoteric behavior.