Atomistic scale investigation of cation ordering and phase stability in Cs-substituted Ba1.33Zn1.33Ti6.67O16, Ba1.33Ga2.66Ti5.67O16 and Ba1.33Al2.66Ti5.33O16 hollandite

The titanate-based hollandite structure is proposed as an effective ceramic waste form for Cs-immobilization. In this study, quantum-mechanical calculations were used to quantify the impact of A-site and B-site ordering on the structural stability of hollandite with compositions BaxCsy(MzTi8-z)O16, where M = Zn2+, Ga3+, and Al3+. The calculated enthalpy of formation agrees with experimental measurements of related hollandite phases from melt solution calorimetry. Ground state geometry optimizations show that, for intermediate compositions (e.g., CsBaGa6Ti18O48), the presence of both Cs and Ba in the A-site tunnels is not energetically favored. However, the decay heat generated during storage of the Cs-containing waste form may overcome the energetics of Ba and Cs mixing in the tunnel structure of hollandite. The ability of the hollandite structure to accommodate the radioparagenesis of Cs to Ba is critical for long term performance of the waste. For the first time, B-site ordering was observed along the tunnel direction ([001] zone axis) for the Ga-hollandite compositions, as well as the intermediate Al-hollandite composition. These compositionally dependent structural features, and associated formation enthalpies, are of importance to the stability and radiation damage tolerance of ceramic waste forms.

waste-forms containing more than one type of B-site cation. Studying the potential for cation and vacancy ordering in hollandite structures will contribute to the understanding of the stability of hollandite over geologic timescales associated with waste storage.
In our previous work, we investigated the compositional effect of the A site of Ga-doped hollandite, where increasing Cs concentration corresponds to increased structural and thermodynamic stability of the hollandite phase 5,6 . Costa et al. 7 discussed B-site dopant effects on the thermodynamic stability of hollandite for a single composition of Ba 1.04 Cs 0.24 (M, Ti) 8 O 16 . In this paper, both the A-site composition and B-site dopant effects are systematically explored using density functional theory calculations.
Density functional theory (DFT) calculations have been used to investigate cation incorporation [8][9][10][11][12] and defect ordering [13][14][15] in ceramic waste form materials. Fluorite-related structures, such as pyrochlore, have proven resistant to radiation damage 16,17 , which has been attributed to the flexibility of the system to withstand disorder 18,19 . Such disorder has been observed experimentally and both confirmed and predicted computationally. Here, DFT calculations are used to quantitatively compare the energetic favorability of different A-site and B-site configurations across a range of hollandite compositions with a fixed A-site occupancy. Ordering observed through theoretical simulations can guide future experimental exploration into preferential cation ordering in hollandite with implications on radiation damage tolerance and long-term stability related to radioparagenesis.
Hollandite compositional range. Hollandite is a titanate structure composed of octahedral B-site cations that form tunnels, which incorporate A-site cations. The A-site is composed of mono and/or divalent cations, specifically Cs + and Ba 2+ . The occupancy of the A-site ranges from one-half to one mole fraction depending on both the A-site and B-site dopant 22 . Hollandite with divalent A-site cations generally have tunnel site occupancy near one-half and less than two-thirds, while hollandite with monovalent cations, such as Na, K or Cs, have tunnel site occupancy of two-thirds to three-quarters, depending on the cation size. For this study the occupancy was 1.33/2 (two-thirds filled), which is achievable for both divalent Ba and monovalent Cs. For the 1 × 1 × 3 supercell used in this study, 4 of the 6 A-sites are occupied. Additional calculations were performed for 2 × 1 × 3 supercells of a subset of hollandite compositions to clarify energetic and structural trends according to increasing Cs concentration. The B-site is predominately composed of Ti 4+ , as well as di and/or trivalent cations for charge-balance. For this study, Zn 2+ and Ga 3+ were used. Al-hollandite was also included for comparison to experimental calorimetry measurements 7 . Table 1 lists the 13 different compositions evaluated in the present work. A range of configurations were studied for the end-member (pure Ba and pure Cs-hollandite) and intermediate member compositions, leading to an understanding of cation ordering at both the A-and B-sites. The configurations evaluated during the current study are available from the corresponding author on reasonable request.
Calculation of thermodynamic properties. The total energy from the DFT calculations does not include any contribution from entropy (either configurational or vibrational); therefore, the total energy can be considered an enthalpy. In fact, the total energy from the CASTEP calculations is a sum of the energies of the zero valent components. The total CASTEP energies are used only for comparison of different hollandite configurations with the same composition. To gleam translatable meaning from the calculated total energy, the enthalpy of formation based on the simple oxides and the enthalpy of mixing across a binary solid solution are computed 23 . In total, three different energetic terms are highlighted for the evaluation of ordering and stability in the hollandite system − 1. Change in total energy for systems with the same composition, but differing configuration (ΔE conf ), 2. Enthalpy of formation for a given configuration based on formation from the simple oxides (ΔH f oxides ), and The enthalpy of mixing for the solid solutions are determined based on a linear combination of the lowest enthalpy end member configurations, where a positive enthalpy of mixing value indicates a preference for a mechanical mixture of the end members 24 . The lowest enthalpy configurations for the end-members are set to zero, and the enthalpy for the other configurations are with respect to the lowest enthalpy configuration. The enthalpy for the intermediate composition is actually the enthalpy of mixing with respect to the lowest enthalpy end-member configurations.
Configurational entropy is considered in the final evaluation of cation ordering for hollandite compositions containing trivalent B-site cations. The mathematical expression for the configurational entropy for the hollandite system is complicated due to the coupled nature of the solid solution, which includes two-site mixing (i.e., mixing on both A and B sites) and an interdependency of the A-site and B-site composition. Following the derivation for the configurational entropy for the plagioclase [(Ca 1−x Na x )(Al 2−x Si 2+x )O 8 ] coupled-substitution solid solution 25 , the configurational entropy for the 1 × 1 × 1 unit cell of hollandite, with the number of mixing components N = 4, is defined in Equation 2. . The configurational entropy is divided into two sets of terms, the first describing the A-site mixing and the second describing the B-site mixing, where the total number of occupied A-site cations is 1.33 and the total number of B-site cations is 8. An additional term may be added to consider the ordering of the vacancy in the A-site, such that N = 5 and = + + x mole fraction of Cs Ba Cs vancany . The vacancy concentration is fixed in this series of calculations; therefore, the vacancy ordering term simply increases the magnitude of the configurational entropy by a factor of 9.481 J/molK according to Equation 3, where the A-site occupancy is 2/3.

Results and Discussion
B-site ordering. B-site ordering was evaluated across the Ba-Cs solid solution with varying B-site dopants (Zn 2+ , Ga 3+ , and Al 3+ ). Figure 1 is a schematic representing a single tunnel layer formed by the B-site cations, as well as the tunnel based on the stacking of the tunnel layers. The A-site cations are removed for clarity. Within a tunnel layer, the B-sites are numbered 1-8, while the tunnel layers within the 1 × 1 × 3 supercell are denoted A, B, and C. The configuration chart in Fig. 1 shows the general format used in describing the B-site configuration, where only the B-site dopants (i.e., non-Ti atoms) are noted. Figure 2a shows the enthalpy of mixing for Zn-hollandite across the Ba-Cs solid solution, where the spread in the data for each composition shows the enthalpy range for varying B-site configurations. Vibrational and configurational entropy are not directly computed in these calculations; therefore, the enthalpy is equivalent to the total calculated energy. The range in energy associated with changing B-site configurations is detailed in Fig. 2b 16 , respectively. Following the configuration chart described in Fig. 1, the lowest energy configuration, highest energy configuration, and an intermediate configuration are highlighted for each composition.
Generally, the lowest energy configuration is one in which the B-site dopants are arranged with the highest symmetry. For example, the Ba-end member for the Zn-hollandite, which contains four Zn atoms in the 1 × 1 × 3 supercell, has a lowest energy configuration in which the Zn are evenly spaced around the eight B-sites of a single tunnel layer (Fig. 2b). The Zn dopants prefer to order within a single tunnel layer (e.g., lowest energy configurations for all Zn-doped compositions tested in Fig. 2b-d) rather than along the tunnel direction (e.g., highest energy configurations for all Zn-doped compositions tested in Fig. 2b-d). Interestingly, Cs 1.33 Al 1.33 Ti 6.67 O 16 prefers the same B-site configuration, where the B-site dopant is symmetrically spaced around one tunnel layer (e.g., Fig. 2b). Meanwhile, the Cs 1.33 Ga 1.33 Ti 6.67 O 16 composition prefers have Ga ordered along the tunnel direction (i.e., the [001] zone axis). The difference in the B-site ordering for these systems stems from the difference in B-site dopant ionic radii, where the Ga 3+ (R Ga = 0.76 Å) is most similar to Ti 4+ (R Ti = 0.745 Å), while the Al 3+ (R Al = 0.675 Å) and Zn 2+ (R Zn = 0.88 Å) differ from the Ti 4+ by 9% and 18%, respectively 26 .  To maintain charge balance, the trivalent B-site dopants are incorporated with a higher concentration than the divalent B-site dopants. The increase in dopant concentration leads to changes in the favorability of the B-site configurations. Optimal configurations still maintain symmetry; however, some deviations are observed related to the size of the dopant cation. For example, for the Ba-end member (eight M 3+ dopants), the lowest energy Ga-hollandite configuration prefers Ga dopants stacked in adjacent layers, while the Al-hollandite prefers the Al dopants to be shifted in adjacent layers (Fig. 3). The larger change in ionic radius between the Al 3+ and Ti 4+ (|∆R| = 0.07 Å), as compared with the Ga 3+ and Ti 4+ (|∆R| = 0.02 Å), likely causes the Al 3+ in adjacent layers to be shifted.
The intermediate trivalent cation compositions contain six B-site dopants, and therefore, the symmetry is inherently broken. Six cations cannot be symmetrically spaced across eight sites. In these cases, the B-site dopants align along the tunnel direction, maintaining two-fold rotational symmetry, rather than the standard 4-fold rotation in the 4-dopant systems.
A-site ordering. As expected, the A-site cations are distributed equally between tunnels within a unit cell; that is, for an occupancy of 1.33, two cations and one vacancy are located in each tunnel within the 1 × 1 × 3 and 2 × 1 × 3 supercells. For intermediate compositions, Cs and Ba prefer to arrange homogeneously within the tunnels of the 1 × 1 × 3 supercell, such that one tunnel contains two Cs atoms and the other contains two Ba atoms. For example, the intermediate composition Ba 2 Cs 2 Ga 6 Ti 18 O 48 indicates that Ba and Cs prefer to be arranged homogeneously, such that each tunnel contains only one cation type; the homogenously ordered configuration is 0.0356 eV lower in energy than the configuration with A-site cations heterogeneously arranged. This energy difference is relatively low, as discussed below, and may be overcome with slightly elevated temperatures.
Cation mobility along the tunnel, and eventually out of the solid phase, is largely dependent on the distribution of the A-site cations and A-site vacancies, as well as the relationship between the A-site cations/vacancies and the B-site dopants. In general, cation-vacancy ordering on the A site is a function of site occupancy with higher occupied tunnels exhibiting no superlattice peaks in diffraction studies 27 . The minimum energy configuration for the CsAl-hollandite (Fig. 4a) has A-site ordering of Cs V Cs/V Cs Cs and a B-site ordering in which the four Al atoms are evenly spaced around a single tunnel layer. Shifting two Al such that the Al remain evenly spaced but span two tunnel layers results in an increase (less favorable) in energy of 0.0695 eV (Fig. 4b). More significantly, when the evenly spaced Al atoms are all located in a different tunnel layer (i.e., the middle layer rather than the outside layer) the total energy for the system increases by 0.1304 eV (Fig. 4c). This increase in energy with change in Al tunnel layer filling suggests that the relationship of the B-site substitutions with the A-site vacancy may play a noticeable role in changing the total energy of the system. The observation of the impact of vacancy ordering in the A-site on the overall system energy is supported when comparing different A-site ordering with the same B-site configuration. For example, for the case of Al evenly spaced in alternating tunnel layers (Fig. 4b), the A-site configuration Cs V Cs/V Cs Cs is favored over Cs Cs V/V Cs Cs by 0.0920 eV. Compositional controls on tunnel diameter. The tunnel diameter, as described here by the diagonal O1-O1 distance, is a key factor that determines the mobility of tunnel cations 6 . The diagonal O1-O1 distance forms a bottleneck that inhibits cation mobility and is affected by both A-site and B-site composition. In general, larger A-site cations will expand the tunnel framework and increase the O1-O1 distance. The average O1-O1 distance increases with increasing Cs concentration (Fig. 5) because the number of tunnels containing Cs increases relative to the number of tunnels containing Ba. Both the ionic radius and oxidation state of the A-site cation contributes to the O1-O1 distance. The ionic radii of cubically-coordinated Cs + and Ba 2+ are 1.88 Å and 1.56 Å 26 , respectively. The lower formal charge of Cs (1+), as compared with Ba (2+), leads to longer Cs-O1 distances than Ba-O1 distances, and thus, longer overall O1-O1 distances across tunnels containing Cs. As noted from neutron diffraction studies, and confirmed with DFT calculations, the O1-O1 distance remains short enough to trap Cs in the hollandite tunnel, limiting cation diffusion 5 .
The oxidation state and ionic radius of B-site dopants also affect the O1-O1 distance. With increasing Cs + concentration, fewer B-site dopants are required for charge balance, decreasing the impact of the B-site dopant on the O1-O1 distance. For example, the ionic radius of octahedrally-coordinated Zn 2+ (0.880 Å) is 18% larger than Ti 4+ (0.745 Å) 26 , but due to the large difference in oxidation state, only one Zn 2+ dopant is needed for charge balance in the 1 × 1 × 3 supercell of the Cs-end member. Thus, the average O1-O1 distance for the Cs-end member of the Zn-hollandite, is similar to the Al-and Ga-hollandite (Fig. 5).
For cations with the same oxidation state (i.e., Ga 3+ and Al 3+ ), the impact of the B-site dopant on the O1-O1 distance is limited, yet consistent. While the difference in the average O1-O1 distance between the Ba-and Cs-end members of the Ga-hollandite solid solution is 0.386 Å, the difference between the Ga-and Al-hollandite for any A-site composition is, at most, 0.078 Å. Table 2 shows the average O1-O1 distances for O1 atoms bonded to two Ti atoms (O1-M-O1 Ti/Ti ), as well as to one Ti and one B-site dopant (O1-M-O1 Ti/B ), for the lowest energy configuration for the 50% Cs composition. The differing trends in the M-O1 vs. Ti-O1 distances for the Al-and Ga-hollandite are due to the size of the ionic radius of the B-site dopant, where the smaller Al dopant leads to shorter M-O1 distances, and in turn, a longer average O1-O1 Ti/M distance. The structure of the Ga-hollandite is more distorted than the Al-hollandite, as seen by the difference in the O1-O1 distances across the Ba-containing tunnel (Fig. 6). Noteably, the O1-O1 distance between O1 bonded to both Ga and Ti are shorter than the O1 bonded to only Ti. The same effect is not observed for the Al-hollandite intermediate system. Therefore, while the A-site cation (Cs vs. Ba) is the dominant control on the symmetry and stability of the hollandite tunnel structure, the B-site cation contributes to localized deviations in symmetry.

Enthalpy of formation.
The enthalpy of formation based on the oxide components was calculated for the lowest energy configuration for all compositions studied and compared with available calorimetry measurements for similar compositions (Fig. 7). Early calculations on the enthalpy of formation for Ga-hollandite 5 are refined in this study with improved computational parameterization, particularly for the simple oxide phases, and increased number of cation configurations evaluated. The calculated enthalpy of formation for Ba 1 16 is −207.80 kJ/mol. The difference between the calculated and measured enthalpy of formation is about 15%, which is in good agreement considering the difference in composition and complexity of the crystalline system. Similar differences have been observed between calculated and measured oxides in literature. For example, the calculated and measured enthalpy of formation for metal uranates differ by about 10-40%, where limitations due to the computational theory are noted only for the system with 40% deviation from measurements 11,28 . Further, the calculated enthalpy of formation for the Ba-Al-hollandite end member ranged from −200.75 kJ/mol to −239.55 kJ/mol depending on the A-site and B-site configuration.
The enthalpy of formation decreases with increasing Cs concentration, where the total system charge is balanced by increasing the concentration of B-site dopants. Therefore, with higher Cs concentration, the deviation of the average B-site radius from 0.745 Å (ionic radius for octahedral Ti 4+ ) increases. For systems with different types of B-site dopants, one might expect the enthalpy of formation to increase (become less favorable) with a greater difference in ionic radii between M and Ti 4+ . Further, the concentration of M-substitutions on the B-site contributes to the overall lattice strain and ultimately impacts the enthalpy of formation. The percent change between the M and Ti 4+ ionic radii is 2%, 9%, and 18% for Ga 3+ , Al 3+ , and Zn 2+ , respectively; however, the enthalpy of formation is highest for Ga-hollandite and lowest for Zn-hollandite. This is likely due to the B-site valance, where the Zn-hollandite system requires fewer B-site substitutions to maintain charge neutrality.    Implications for Long-term Disposal of Nuclear Waste. The B-site dopant has been noted for impacting the tunnel diameter 7 , where the lattice parameter increased with increasing average B-site ionic radius. The B-site ordering has the largest impact on structural stability of the system, where the energy difference between the most favorable and least favorable B-site configuration for the Cs-hollandite end members range between 0.26 eV and 0.36 eV. The propensity for B-site ordering also impacts the expected radiation damage tolerance of the material. Structures that are initially ordered, but have the ability to accommodate lattice disorder are known to possess larger radiation stability 17,19 . Limited studies are available on the radiation damage behavior of hollandite, however ordered B-site structures are expected to be more resistant to damage 29,30 .
Recently published neutron diffraction data for hollandite with high Cs-loading shows that the Cs concentration has a larger effect on the tunnel diameter and overall system symmetry than the B-site dopant type 5 . The calculated enthalpy of formation supports the conclusion that Cs stabilizes the formation of the hollandite system. Overall, the calculated variance of the enthalpy across the Ba-Cs solid solution decreases with increasing Cs concentration, further indicating the stabilization of the hollandite structure due to Cs in the tunnel cavity. The stabilization of the hollandite structure with increased Cs concentrations has positive implications for Cs-waste loading in hollandite. Thermodynamic comparisons with competitive Cs-host phases is necessary to quantify the expected maximum concentration of Cs in hollandite formed in a multi-phase ceramic.
Over time, 137 Cs undergoes beta decay (τ ½ = 30.08 yrs) to become 137 Ba; therefore, understanding the arrangement of Cs and Ba on the A-site provides insight into the long term stability of hollandite due to the radioparagenesis of the waste components 31 . Ba and Cs prefer to segregate into different tunnels by about 0.036 eV. The implication of this energetic favorability can be estimated using a Boltzmann conversion of energy to temperature (E = kT, where E is energy, k is the Boltzmann constant of 8.617 × 10 −5 eV/K, and T is temperature in K), where 0.036 eV equates to about 418 K (145 °C). While elevated, this temperature is far below the temperature used during waste form processing (1000-1500 °C) 2,32,33 . Commercially-derived high-level waste could generate enough self-heat from decay of fission products to result in an initial storage temperature as high as 600 °C. Even 100 years after emplacement, the waste temperature may be as high as 300 °C 34 . Therefore, enough energy may be supplied by the elevated temperature such that, upon decay, Ba may remain stable in the Cs-dominant tunnel of hollandite. Table 3 highlights the quantum-mechanically determined energy for different structural conditions examined in this work and the associated temperature based on the aforementioned Boltzmann conversion. Waste processing and storage temperature ranges are included for reference.
The smallest energetic difference for A-site vacancy ordering is 0.091 eV (~780 °C), which is associated with A-site ordering for the Cs-Al hollandite system. This energy difference is large enough to have implications for A-site cation mobility since it falls above the temperature range expected for storage of these waste forms. The occupancy explored in this study (1.33/2) is slightly higher than previous experimental work (1.24/2) in order to  compare detailed computational geometry optimizations with experimental structural and mobility data. Thus, further analysis of A-site occupancy is critical to fully understand the cation mobility advantages of Cs-doped hollandite recently indicated 6 . The ordering of A-site (tunnel) cations has been both predicted 35 and confirmed computationally and experimentally 27,[36][37][38] .
The entropic contributions to thermodynamic stability were evaluated based on Equations 2 and 3. Figure 8 shows that the highest configurational entropy (57.3 J/mol K for N = 5 system) is associated with 0.296 mole fraction of Cs due to the coupling of the A-site and B-site mixing in this system. Further consideration of the vacancy ordering does not change the coupling of the A-site and B-site mixing for fixed occupancy systems. Thus, the N = 5 curve results in an increase in the magnitude of the configurational entropy associated with the additional mixing options on the A-site. The enthalpic favorability of B-site ordering and vacancy ordering are on the order of magnitude of the configurational entropy, particularly in the temperature range of waste form storage (below 600 °C); thus, the ordering in these systems can contribute to the overall stability of the ceramic waste form.
This study is the first to predict B-site (framework) ordering in Cs-loaded hollandite phases. The models presented here can guide future experimental efforts to confirm Ga ordering along the [001] zone axis (i.e., the tunnel direction). Following the trends in ionic radii, B-site dopants with ionic radii close to octahedral Ti 4+ (0.745 Å) may also order along the tunnel direction. Further, the present work shows the coupling of B-site and A-site ordering, which, in turn, impacts the mobility of the A-site cations. Ultimately, quantification of the A-site cation mobility is critical for predicting long term potential for radionuclide release.  Table 3. Key structural conditions explored in quantum-mechanical calculations, calculated energy difference between most favorable and next most favorable configuration, and temperature (from Boltzman conversion) necessary to overcome the enthalpic favorability for given cation arrangement. Waste processing and storage temperature ranges are included for reference. + Temperature for the different cation configuration scenarios (rows 1-3) are calculated based on a Boltzman relationship of T and energy.