The The Effect of Ionic Defect Interactions on the Hydration of Yttrium-doped Barium Zirconate

The use pair for the that defects in the v O∙∙ , proton H i∙ , and ion Y Zr′ . The Y Zr′ − Y Zr′ interactions were not modelled ABSTRACT This supplement extends the main article “ The Effect of Ionic Defect Interactions on the Hydration of Yttrium-doped Barium Zirconate ”. The pair interaction models calculated with density functional theory are quantified for the short-range and long-range cases, including symmetry reduced geometries and energies. Information on the applied in-house Monte Carlo software MOCASSIN concerning the Metropolis annealing routine is provided. Additionally, the results using the long-range interaction model for BZY10, BZY20, and BZY30 not contained in the main article, the derivation for the mass action law as applied in the main article, the derivation of the calculation of 𝐹 from MMC results, and Δ𝐹 int′ raw data tables are provided.


BZO structure
The exact BZO unit cell description used for the MMC simulations and in the following pair interaction tables is summarized in Table 1. The affiliated space group of the perovskite BZO is Pm3 ̅ m (No. 221) with 48 symmetry operations and the unit cell parameter is 4.23 Å. The undoped barium sites were omitted during simulation and the chosen H i • site is at 1.00 Å distance to the oxygen site.

Pair interaction models
The applied DFT interaction models use pair interactions for the MMC active species that are defects in the ideal BZO host matrix; namely, oxygen vacancy v O •• , interstitial proton H i • , and yttrium ion Y Zr ′ . The Y Zr ′ − Y Zr ′ interactions were not modelled as the Zr sublattice occupants were uniformly distributed at random prior to the simulations. No MMC equilibration of the cation sublattice was conducted and thus the Y Zr ′ − Y Zr ′ interactions had no influence on the results. The relevant interactions are summarized in Table 2 to Table 6. Here, interactions with ID values containing an "a,b,c,…" label mark interactions at the same distance but with different reference geometries, and IDs extended with a second "1,2,3,…" label mark interactions at the same distance where the affiliated extended vector sets in the unit cell are chiral to each other. This means, they cannot be fully mapped onto each other by the space group symmetry operations without an additional chirality-breaking operation and it is not possible to calculate distinct values for the two sets using DFT. Each table has a fixed occupation pair, a fixed site for the first position and a variable position for the second site. Geometric data is provided as (a,b,c) coordinates in the affine coordinate system of the unit cell. The last three columns of the tables define the number of interactions per site after P1 translation of the symmetry reduced structure for the origin site (#Origin) and the partner site (#Partner), and the model affiliations with (S) for short-range and (L) for long-range model, respectively. The last decimal for each energy value is uncertain.     Supplement defining the probability of changing from one system state to another was modified with a factor ( ) as shown in expression (1). Here, min is the minimal target temperature and is the current temperature. In this formalism, = 0 and = 1 represent infinite temperature ∞ and the final temperature T min , respectively, and the solver gradually simulates [0 → 1] for one system. A description of the implementation can be found in the affiliated publication.

Interactions between interstitial protons
A minimal C pseudocode implementation of the distribution recording routine is illustrated in Code 1. Here, the 'DoMmcCycle(…)' function refers to the usual MMC cycle processing the supplied alpha according to expression (1).
Code 1. Minimal C pseudocode of the MMC annealing routine that logs MMC distributions. The SITE_BLOCKED flag signals that no transition rule was found, these cases are omitted from the total counts so that only ACCEPTED and REJECTED attempts are recorded.

Results for and
The results of Δ int and Δ int for BZY10, BZY20, and BZY30 using the long-range interaction model are illustrated in Figure  1, Figure

4.2
Results for ′ and The results of Δ int ′ and int for BZY10, BZY20, and BZY30 using the long-range interaction model are illustrated in Figure   4 and Figure 5.

Results for [ • ] against
The results of [OH O • ]( H 2 ) for BZY10, BZY20, and BZY30 using the long-range interaction model are illustrated in Figure 6.

Derivation of the mass action law with
Please note that the given derivation of the int formalism uses the incorporation of one [OH O • ] defect according to the incorporation reaction: In contrast, the hydration reaction in the main manuscript is normalized to two [OH O • ] defects as this is the usually used reaction formalism in the literature. The resulting description is equivalent and is normalized to two defects in the last step.
The change in Gibbs free energy of the yttrium-doped barium zirconate system due to the hydration process as written in reaction (2) .
The ideal mixing entropy can be described using statistical mechanics. The only relevant contribution in this context is the change in the mixing entropy of the oxygen sublattice as the mixing entropy of the cation lattices does not change due to hydration. Two subcontributions can be distinguished; the first is the mixing entropy of the hydrated system as described by eq. (2).
The second subcontribution according to eq. (3) is negative and corrects the mixing entropy already present in the dry doped system.
Here, the quantity Ω( ) describes the total number of states with energy . Using the Metropolis Monte Carlo algorithm, it is possible to obtain samples of state distributions ( , ) that fulfill relation (23).
The direct evaluation of and consequently is not possible as the proportionality factor is unknown. When taking a uniform random sample Ω ′ ( ) of the configurational space using MMC, that is, a simulation where each exchange attempt is performed without checking the probability function, the equality relation (24) must however hold true as the ration between Ω( ) and Ω ′ ( ) does not depend on the energy .
Combining equations (24) and (23) yields the relation between the random sample Ω ′ ( ) and the finite temperature MMC distribution ( , ) as shown in expression (25).
Rearranging and integrating the distributions leads to a new factor ′ as shown in (26).
Assuming Ω ′ ( ) and ( , ) contain the same number of state samples, the partition function can be rewritten into equation (27) where Θ is the total number of configurations of the system as defined in expression (28).

= ⋅ ∫ ( )d = ′ ⋅ Θ (27) Θ = ∫ Ω( )d (28)
This allows to redefine the free energy as shown in (29) where the first term on the right-hand side is the free energy of the ideal system (Θ) and the second term is the interaction contribution int .
Again, the direct evaluation of ′ is usually not possible as it rarely possible to obtain Ω ′ ( ) and ( , ) that show a statistically significant overlap. At this point, the multistage sampling solution suggested by Valleau and Card 2 that calculates ′ from a series of bridging distributions with small temperature deltas as shown in the main manuscript can be used. This additionally solves an issue that ′ might be large, causing issues with the standard IEEE754 double precision numbers.  12.2 10.8 9.9 9.6 9.5 9.6 9.7 9.7 9.6 9.4 9.1