A model for time-dependent grain boundary diffusion of ions and electrons through a film or scale, with an application to alumina

TL;DR

This paper introduces a comprehensive time-dependent model for ion and electron transport through grain boundaries in thin films, validated by alumina membrane oxygen permeation experiments, revealing microscopic processes and conductivity transitions.

Contribution

The paper develops a self-consistent 3D Poisson-Nernst-Planck model for grain boundary transport, applicable to various materials and phenomena, with detailed simulation of alumina oxygen permeation.

Findings

Simulation captures oxygen flux under different pressures.

Reveals transition between p-type and n-type conductivity.

Quantifies validity of simplified permeation models.

Abstract

A model for ionic and electronic grain boundary transport through thin films, scales or membranes with columnar grain structure is introduced. The grain structure is idealized as a lattice of identical hexagonal cells - a honeycomb pattern. Reactions with the environment constitute the boundary conditions and drive the transport between the surfaces. Time-dependent simulations solving the Poisson equation self-consistently with the Nernst-Planck flux equations for the mobile species are performed. In the resulting Poisson-Nernst-Planck system of equations, the electrostatic potential is obtained from the Poisson equation in its integral form by summation. The model is used to interpret alumina membrane oxygen permeation experiments, in which different oxygen gas pressures are applied at opposite membrane surfaces and the resulting flux of oxygen molecules through the membrane is measured. Simulation results involving four mobile species, charged aluminum and oxygen vacancies, electrons, and holes, provide a complete description of the measurements and insight into the microscopic processes underpinning the oxygen permeation of the membrane. Most notably, the hypothesized transition between p-type and n-type ionic conductivity of the alumina grain boundaries as a function of the applied oxygen gas pressure is observed in the simulations. The range of validity of a simple analytic model for the oxygen permeation rate, similar to the Wagner theory of metal oxidation, is quantified by comparison to the numeric simulations. The three-dimensional model we develop here is readily adaptable to problems such as transport in a solid state electrode, or corrosion scale growth.

Equations

$$P\propto {\left(P_{{\text{O}}_2}^{II}\right)}^{3/16},$$

$$P\propto {\left(P_{{\text{O}}_2}^{II}\right)}^{-1/6}.$$

$${\eta }_s(r,t)={\mu }_s(r,t)+z_s{e}_0\varphi (r,t)$$

$$J_s=-\frac{D_s{c}_s}{k_{\text{B}}T}\nabla {\eta }_s$$

$$J_s=-D_s\nabla c_s-\frac{D_s{c}_s{z}_s{e}_0}{k_{\text{B}}T}\nabla \varphi$$

$$\rho (r,t)=\sum _s{z}_s{e}_0{c}_s(r,t)$$

$$\varphi (r,t)$$

$${\int }_{\mathbb{V}}\rho (r,t)d^{3}r=0.$$

$$\frac{\partial }{\partial t}{c}_s=-\nabla \cdot J_s+R_s$$

$${\int }_{{\mathbb{V}}_{\text{TRI}}}\frac{\partial {\mathcal{C}}_s}{\partial t}{d}^{3}r={\int }_{{\mathbb{A}}_{\text{GB}}}{J}_s\cdot \hat{n}{d}^{2}r+{\int }_{{\mathbb{V}}_{\text{TRI}}}{\mathcal{R}}_s{d}^{3}r.$$

$$\frac{\partial }{\partial t}{\mathcal{C}}_s=\frac{1}{{\mathbb{V}}_{\text{TRI}}}{\int }_{{\mathbb{A}}_{GB}}{J}_s\cdot \hat{n}{d}^{2}r+{\mathcal{R}}_s$$

$$r$$

$$t$$

$$\frac{\partial }{\partial t}{c}_s$$

$${\kappa }^2\varphi (r,t)$$

$$\kappa =\frac{{\ell }_{\text{D}}}{L},{\ell }_{\text{D}}={\left(\frac{{\epsilon }_0{\epsilon }_r{k}_{\text{B}}T}{e_0^2{c}_{\text{ref}}}\right)}^{1/2}$$

$${\lambda }_{\mathrm{D}}={\left(\frac{{\epsilon }_0{\epsilon }_r{k}_{\text{B}}T}{e_0^2{\sum }_s{z}_s^2{c}_s^0}\right)}^{1/2}$$

$$J_s^{i+1/2,j,n}=$$

$$-D_s\frac{c_s^{i+1,j,n}\text{B}\left(-z_s\Delta {\varphi }^{i,j,n}\right)-c_s^{i,j,n}\text{B}\left(z_s\Delta {\varphi }^{i,j,n}\right)}{\Delta x_1^i}$$

$$c_s^{i,j,n}-c_s^{i,j,n-1}+\frac{\Delta t_n}{\Delta {\overline{x}}_1^i}\left(J_s^{i+1/2,j,n}-J_s^{i-1/2,j,n}\right)+$$

$$\frac{\Delta t_n}{\Delta {\overline{x}}_3^j}\left(J_s^{i,j+1/2,n}-J_s^{i,j-1/2,n}\right)-\Delta t_n{R}_s^{i,j,n}=0$$

$$\varphi (r,t_n)={\kappa }^{-2}\text{sideset}{}^{\prime }{\sum }_j^{\mathbb{V}}\frac{{\sum }_s{z}_s{c}_s^{j,n}{\nu }_j}{|r-r_j|}$$

$$\varphi (r_i,t_n)=\sum _j{q}^{j,n}{g}_{ij}.$$

$${\rho }_g(r,r_j,t_n)=\rho (r_j,t_n){\left(\frac{\beta }{\pi }\right)}^{\frac{3}{2}}\exp (-\beta |r-r_j{|}^2)$$

$${\varphi }_g(r,r_j,t_n)={\kappa }^{-2}\frac{q^{j,n}}{|r-r_j|}\text{erf}(\sqrt{\beta }|r-r_j|).$$

$${\varphi }_g(|r-r_j|\to 0,t_n)=2{\kappa }^{-2}{\left(\frac{\beta }{\pi }\right)}^{\frac{1}{2}}{q}^{j,n},$$

$$\frac{1}{2}{\text{O}}_2(g)+2{\mathrm{e}}^{-}$$

$$|z_{\nu }|{\mathcal{C}}_{\nu }\simeq {\mathcal{C}}_{\overline{\nu }}$$

$${\mathcal{C}}_{\nu }=f_{\nu }{P}_{{\text{O}}_2}^{n_{\nu }}$$

$$n_{\nu }=\frac{-z_{\nu }}{2|z_{\text{O}}|\left(|z_{\nu }|+1\right)}$$

$$f_{\nu }={\left({\left|z_{\nu }\right|}^{-|z_{\nu }|}{K}_{\text{S}}^{3-|z_{\nu }|}{K}_{\text{eh}}^{|z_{\nu }|(|z_{\nu }|-2)}{K}_{\text{Al}}^{-z_{\nu }/2}\right)}^{\frac{1}{|z_{\nu }|+1}}$$

$$P$$

Full text

††journal: Acta Materialia

1 Introduction

Thin films of insulating material at metal – gas and metal – liquid interfaces accomplish a range of service functions in materials technology. Common examples are functional ceramics in electronics, energy related applications and sensors. Thin films formed by surface oxidation of a metal can have either beneficial or corrosive effects. Alumina and chromia formed by thermal oxidation are examples of protective oxide films, which find application in thermal barrier coatings, and can be engineered for durability by additions of rare earth elements [1, 2, 3, 4]. They can also grow in an uncontrolled manner, adhering weakly to the metal and allowing corrosion to proceed. Metal oxidation is a heterogeneous process, consisting of multiple steps, involving the dissociation of molecular oxygen, transport through the growing oxide layer, and the reaction between oxygen and metal atoms. Generally speaking a growing oxide layer on the metal surface requires either metal atoms, normally ions, to be transported to the oxide – gas interface to sustain the oxidation reaction, or alternatively oxygen atoms or ions to be transported through the oxide to the metal – oxide interface to sustain an internal oxidation reaction. Furthermore, the transport is usually thought to be effected by diffusion of cation or anion vacancies. Both processes may proceed, depending on the material under consideration and on the environmental conditions such as temperature and oxygen partial pressure, $P_{\text{O}_{2}}$. While the slowest process is rate-determining for consecutive processes, like dissociation of oxygen molecules and their transport through the oxide, the fastest process is rate-determining for parallel processes, like bulk and grain boundary diffusion through the oxide. The measured oxygen and aluminum diffusion coefficients in $\alpha$-alumina are found to be several orders of magnitude greater at the grain boundaries than in the bulk material [5, 6, 7]. The fact that grain boundaries provide the dominant transport mechanism underlines the importance of including their geometric and transport properties in a realistic model of the process. Since vacancies in strongly ionic oxides such as alumina or chromia are charged species relative to the perfect crystal, the fluxes of these species carry an electric current, which in the usual scenario of steady-state growth is not sustainable, unless compensated by an equal and opposite current of electrons or holes, as described by the classic model of Wagner [8, 9]. The prediction of the growth behaviour of thin films, and its influence on the material or device performance, requires us to describe the mixed ionic, electronic transport through the films, while taking their grain boundary structure into account.

Because of the widespread importance of alumina films [10], and since it is a relatively well characterized material, we focus on alumina films for the validation of our modelling approach. Moreover, a recent series of permeation experiments for $\alpha$-alumina polycrystal membranes, e.g. [11], conducted for different combinations of applied oxygen gas pressures at high temperatures, provides an ideal test case for our transport model. These experimental results will be summarized briefly in the following section.

1.1 Brief review of oxygen permeation and diffusion experiments

Permeation rates of oxygen through a polycrystalline membrane of alumina have been reported in the literature [11, 12], and cover a range of oxygen partial pressures. Scanning Electron Microscope (SEM) imaging of the films prepared under different applied pressures strongly suggests that mass transfer occurs along grain boundaries. The thermodynamic driving force in these experiments is the difference in the oxygen chemical potential between the two membrane surfaces $\Delta\mu_{\text{O}_{2}}=k_{\text{B}}T\ln\left(P_{\text{O}_{2}}^{II}/P_{\text{O}_{2}}^{I}\right)$. Figure 1 shows a schematic of membrane permeation experiments.

The experiments included nominally pure $\alpha$-alumina polycrystals [11, 13], doped $\alpha$-alumina polycrystals [14, 15, 16, 17], and nominally pure $\alpha$-alumina bicrystals [18], with temperatures of $\sim 1700-2000\,\text{K}$. A simple analysis of the permeation rate data in reference [11], assumed a model of one-dimensional, steady-state diffusion, in which either Al or O is transported by vacancy migration, depending on the absolute magnitude of the applied oxygen pressure.

In the non-doped polycrystalline alumina experiments [11], when applying high oxygen pressures at surface ($II$), $P_{\text{O}_{2}}^{\text{hi}}=10^{3}-10^{5}\,\text{Pa}$, while keeping surface ($I$) at $P_{\text{O}_{2}}^{I}=1\,\text{Pa}$, grain boundary ridges formed on the $P_{\text{O}_{2}}^{\text{hi}}$ surface and grain boundary trenches were observed on surface ($I$). Applying a low oxygen pressure at surface ($II$), $P_{\text{O}_{2}}^{\text{lo}}=10^{-5}-10^{-8}\,\text{Pa}$, while keeping surface ($I$) at $P_{\text{O}_{2}}^{I}=1\,\text{Pa}$, no grain boundary ridges are formed and only grain boundary trenches are observed. Since the oxygen permeation rates of a single-crystal alumina wafer were below the measurable limit and as the visible surface growth and “dissolution” proceeds at grain boundaries it is reasonable to assume that grain boundaries dominate the transport [11]. The oxygen permeation rates, $P$, for fixed $P_{\text{O}_{2}}^{I}=1\,\text{Pa}$ were found to follow distinct power laws [11]; in the limit of $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$,

[TABLE]

and in the limit of $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{lo}}$,

[TABLE]

The power laws and the pressure dependent formation of the grain boundary ridges have led to the interpretation of the experiments in terms of, aluminum vacancy transport being dominant in oxygen chemical potential gradients with high oxygen pressure magnitude, $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$, and oxygen vacancy transport being dominant in the case of $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{lo}}$.

Indeed, the rational power laws appear in the theory as a direct consequence of the +3 and -2 ionic charges of the ions, assuming that the negative of these charges is carried by each vacancy, with a counter-current of electrons or holes, and no time-dependence of the fluxes (the steady-state assumption) or net local charge densities within the grain boundaries. These assumptions are discussed further in section 4.1. Furthermore, only aluminum vacancy transport can lead to ridge formation, which is observed by SEM imaging in the case of $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$, supporting the above interpretation. The switch-over in dominant point defect species in the grain boundary has been termed “$p-n$ transition” in the literature [5, 19].

The diffusion coefficients determined from alumina bicrystal experiments with $P_{\text{O}_{2}}^{II}/P_{\text{O}_{2}}^{I}=10^{5}\,\text{Pa}/1\,\text{Pa}$ for several distinct grain boundary types have been found compatible with those measured in polycrystalline samples [18]. The bicrystal diffusion coefficients were calculated from the grain boundary ridge volume and a caveat regarding this approach is that the formation of the ridges on the $P_{\text{O}_{2}}^{\text{hi}}$ side does not necessarily imply an exactly equivalent mass transport from the opposite side of the membrane, since oxide can be displaced by formation of internal pores in the subsurface region of the crystals, as has indeed been observed in some polycrystal permeation experiments [14].

1.2 Scope of the paper

The above level of analysis leaves several open questions, e.g. does the grain boundary diffusion mechanism, with the associated inhomogeneity of fluxes and electric fields, map accurately onto a 1D diffusion problem? And what are the magnitude and roles of the surface and interface charges, the electric fields, currents, space charges, and transients that are all believed to be present in a three-dimensional film, traversed by grain boundaries?

To address these questions we have modelled the transport of oxygen through a planar film by making an idealised representation of the grain structure, in which we suppose the grains to be columnar, with identical and perfectly hexagonal cross sections, see figure 2. We describe in this paper the set of coupled reaction-diffusion equations we have used to model the oxygen permeability across the membrane, our method of solving them, and the results we have obtained for the model alumina membrane. Within this geometry, grain boundary transport through the film of charged cation and anion vacancies, electrons, and holes is simulated, while reactions with the environment constitute the boundary conditions. We expect our 3D model to be applicable to different materials and to be readily adaptable to the problem of oxide scale growth, which is very similar to the problem of oxygen diffusion through a membrane, the most significant difference being the boundary conditions.

The equations describing time-dependent diffusion of ions, electrons and holes through a polycrystalline film, driven by electric fields and defect concentration gradients, cannot be solved analytically in general, and numerical methods must be applied. Models including numerical computations to describe transport through films have been developed for homogeneous films [20, 21, 22], for which a one dimensional model may be a suitable approximation. The symmetry of the present model of idealized columnar hexagonal grains is used to reduce the problem to 2D boundary diffusion along the rectangular boundaries of the hexagonal grains, while explicitly taking account of the long-ranged electrostatic interactions between the charged species within the 3D structure. In order to be able to describe transient behaviour and time-dependent environments, time-dependent boundary conditions are taken into account, but the movement of the boundaries is neglected. This means surface charges can build up or be depleted as a function of time, due to the reactions with oxygen in the environment, in our case a prescribed oxygen partial pressure, and the delivery of charged species to the surfaces from the grain boundaries.

The system of equations used in our model to describe the fluxes of point defects, electrons and holes and their Coulomb interaction, is mathematically equivalent to the drift-diffusion (DD) equations applied in semiconductor device simulations for electron and hole transport [23, 24], and to the Poisson-Nernst-Planck (PNP) system, which is used for ion channel simulations [25, 26] and other electrochemical applications [27]. Unlike most computational methods of solution for the DD and PNP equations, which solve the Poisson equation in differential form and often only consider the steady-state solution, the solution method developed here allows for time-dependent calculations and the Poisson equation is solved in its integral form, taking into account the long-range Coulomb interaction within the 3D structure.

2 The hexagonal cell model

2.1 2D-periodic hexagonal prism grain structure

The model is developed for films with columnar grain structure. The columnar grains are idealized as hexagonal prisms, with the rectangular faces denoting the grain boundaries, see figure 2. The hexagonal prisms are periodically repeated to construct a slab of infinite extent in two dimensions.

Transport through the slab is presumed to be dominated by that through the rectangular grain boundaries, located at the interfaces between the cells. The grain boundaries are assumed to be composed of a very thin homogeneous and isotropic medium of finite width $\delta$. The width of the boundary is not a physical width, but rather a theoretical construct, which allows concentrations to be expressed per unit of volume or per atomic site rather than per unit of area.

From the symmetry of the system only the “irreducible zone” of the hexagonal cell, shown in figure 4, needs to be considered as a domain for calculation. It includes a triangular piece of the surface hexagon and half of a grain boundary rectangle. The 2D-periodic tiling enables an explicit calculation of the long-range Coulomb interaction between the charged point defects.

2.2 Equations for grain-boundary transport

The point defect concentrations within the boundary are assumed to be continuous functions of space and time. Based on the local equilibrium hypothesis the electrochemical potential, $\eta_{s}$, of point defects, and electrons or holes of species “$s$”, is given by

[TABLE]

where, $\mu_{s}$ is the chemical potential of species $s$, $\phi$ is the electrostatic potential, $z_{s}$ is the charge (integer number), $e_{0}$ is the positive elementary charge, and $r$ is the position vector. Phenomenologically, the particle flux due to a gradient in the electrochemical potential may be written as,

[TABLE]

which in the ideal solution approximation [28] is equivalent to

[TABLE]

where $D_{s}$ is the diffusion coefficient, and $c_{s}(r,t)$ is the concentration (number per unit volume) of species $s$. Equation 5 is referred to as the Nernst-Planck flux [29, 30], and combines Ohm’s law of conduction and Fick’s law of diffusion. It is used here as the constitutive equation for the description of the point defect transport; magnetic field effects are not considered.

The local charge density is given by

[TABLE]

where the sum is performed over all charged species present. The instantaneous electrostatic potential can be calculated from the charge density by solving Poisson’s equation, which is given here in integral form for a linear dielectric material,

[TABLE]

where $\varepsilon_{0}$ and $\varepsilon_{r}$ are the vacuum and relative permittivity respectively, and the domain of integration $\mathbb{V}$ of the “Coulomb integral” includes the entire system, in which charge densities are non-zero within a slab of infinite extent in two dimensions, composed of identical hexagonal cells and their surfaces. Overall charge neutrality holds for the domain $\mathbb{V}$, at any instant of time

[TABLE]

The continuity equations for the individual species are given by

[TABLE]

where a reaction term, $R_{s}$, has been added to enable processes such as electron and hole recombination within the grain boundary to be described. The continuity equations for the different mobile species are used in the simulations to evolve the defect concentrations in time; they are solved self-consistently with the equation for the electrostatic potential 7, which depends on the charge density $\rho(r,t)$.

To model coupling effects in the species transport we could formulate the dynamics in terms of linear irreversible thermodynamics [31, 32], but this theory does not provide explicit expressions for the constitutive equations including the transport coefficients. The linear constitutive equations, like Fick’s first law, have a phenomenological basis; however, they can also be thought of as laws of inference based on probability theory [33]. For electron and hole transport in semiconductor device simulations a system of equations, mathematically equivalent to equations  5, 7, 9 was first introduced by van Roosbroeck [34]. They are referred to as drift-diffusion equations, and can be derived from the Boltzmann transport equation by either the Hilbert expansion or the moment method [23, 24].

2.3 Time dependent boundary conditions

The boundary conditions are formulated to describe oxide creation and dissolution at the slab surfaces by reaction with the environment. We refer to our system, which includes its upper and lower surfaces, as ‘the slab’. The thin surface layers of the slab are treated as a homogeneous and isotropic medium of thickness $\delta$. The flux of defects between the surfaces and the grain boundary, and the reactions between the slab and the environment, change the concentrations of species $s$ in the surface layers, $\mathcal{C}_{s}$, as expressed by

[TABLE]

Reactions between the slab surface and the environment produce species $s$ in the surface layer at a rate $\mathcal{R}_{s}$. The unit normal, $\hat{n}$, and the integration domains are defined in figure 4. Equation 10 holds for all species “s” and separately for the surfaces ($I$) at $x_{3}=0$ and ($II$) at $x_{3}=L$.

Assuming transport to be much faster across the surface than between surface and grain boundary, uniform defect concentrations are used on the surfaces. This leads to the following simplification of equation 10

[TABLE]

which is used as the boundary condition in the simulations reported below. If the surface transport mechanism and parameters are known it is straightforward to relax the above assumption.

No separate boundary condition is needed for the electrostatic potential since it is calculated by a summation technique, only requiring the instantaneous charge distribution as a function of position. The reaction rates in the grain boundary, $R_{s}$, and on the surface, $\mathcal{R}_{s}$, depend on the application of the model. For the present purpose of describing the alumina permeation experiments the rate equations are derived with the law of mass action, and discussed in more detail in section 4.2.1. The application to oxide scale growth will require a separate boundary condition at the interface between oxide and metal.

The initial conditions also depend on the application. For the oxygen permeation experiments they are discussed in section 4.3.1.

3 Method of solution

An object oriented C++ code has been developed to solve the system of coupled partial differential equations 5, 9, and 11, self-consistently with the Coulomb integral, equations 6, and 7, which is approximated by the summation technique described in section 5. Since the system of equations is nonlinear, involves vastly different rates in the diffusion processes and the reactions, and involves different length scales characterized by the Debye length, defined below, and the system size, this is a challenging computational problem.

3.1 Dimensionless equations

To obtain dimensionless equations for the numerical calculations the variables, parameters, and fields are scaled as follows [24]:

[TABLE]

where L is the thickness of the scale; $D_{\text{ref}}=\max\{D_{s}\}$ is the largest diffusion coefficient for all species; $c_{\text{ref}}$ is a suitable reference concentration; and $e_{0}$ is the elementary positive charge. To avoid unnecessarily heavy notation, the original quantities were denoted with the “tilde” mark over the symbol in the definition of the scaling, and the dimensionless quantities without the tilde are used in the following.

With the scaling defined by equations 12, the system of transport equations and the Coulomb integral can be brought into dimensionless form

[TABLE]

where the dimensionless parameter, $\kappa$,

[TABLE]

and the reference screening length, $\ell_{\text{D}}$, are introduced. $L$ denotes the thickness of the slab, see figure 2. The Debye screening length is defined by

[TABLE]

where $c_{s}^{0}$ are the spatially uniform concentrations obtained in the limit of $T\rightarrow\infty$, while holding the total number of each species, $N_{s}$, constant. In the application considered in this work the oxygen gas chemical potentials are fixed at the surfaces, the concentrations, $c_{s}(r,t)$, are independent variables, and the total number of each species depends on time, $N_{s}(t)$; however, $c_{\text{ref}}$ is chosen such that $\ell_{\text{D}}$ has similar magnitude to $\lambda_{\mathrm{D}}$, and $\ell_{\text{D}}$ is therefore referred to as the reference screening length.

Overall charge conservation and zero total charge within the system are maintained during the evolution of the concentrations, while charge is redistributed within the grain boundary and moved in or out of the surfaces by the fluxes.

3.2 Discretization of the transport equation

The finite difference method is used to discretize the continuum equations in space and time [24, 35]. For the spatial discretization a rectangular mesh is used, and it covers the irreducible zone of the hexagonal cell, indicated in figure 4. The concentrations and electrostatic potential at mesh node $(i,j)$ at position $r(x_{1}^{i},x_{3}^{j})$ in the grain boundary ($x_{2}=0$) and discrete time $t_{n}$ are denoted by, $c_{s}(x_{1}^{i},x_{3}^{j},t_{n})=c_{s}^{i,j,n}$, and $\phi(x_{1}^{i},x_{3}^{j},t_{n})=\phi^{i,j,n}$, respectively, where $i=1,2,...,N_{1}$ and $j=1,2,...,N_{3}$.

The mesh spacings near the surfaces need to be significantly smaller than the reference screening length to resolve the behaviour near the surfaces correctly. However, the whole grain boundary cannot be meshed with such a fine spacing in the $\hat{e}_{3}$ direction, because the summation technique used for calculating the Coulomb interaction would become too computationally expensive. Therefore, a layer-adapted mesh is used in the $\hat{e}_{3}$ direction.

The Nernst-Planck flux in equation 5 is approximated with the Scharfetter-Gummel discretization scheme [36]; in one dimension for fixed $j$ it is given by

[TABLE]

where $J_{s}^{i+1/2,j,n}$ is the flux of species $s$ between node $i$ and $i+1$ at time $t_{n}$, $\Delta x_{1}^{i}=x_{1}^{i+1}-x_{1}^{i}$ is the length of the interval, $\Delta\phi^{i,j,n}=\phi^{i+1,j,n}-\phi^{i,j,n}$ is the potential difference between the mesh nodes at time $t_{n}$, and $\text{B}(x)=x/(\exp(x)-1)$ is the Bernoulli function [24].

For the time stepping the continuity equation 9 is discretized in implicit form,

[TABLE]

where $\Delta\overline{x}^{i}_{1}=(\Delta x^{i}_{1}+\Delta x^{i-1}_{1})/2$, $\Delta\overline{x}^{j}_{3}=(\Delta x^{j}_{3}+\Delta x^{j-1}_{3})/2$, and $\Delta t_{n}=t_{n}-t_{n-1}$. The initial time step, $\Delta t_{1}$, is chosen such that the discretizations in space and time have a similar order of accuracy, hence $\Delta t_{1}\approx(\Delta x_{\mathrm{min}})^{2}$ where $\Delta x_{\mathrm{min}}$ denotes the smallest mesh element. $\Delta x_{\mathrm{min}}$ is chosen to have the same size as the surface layer thickness, $\delta$, therefore typically $\Delta t_{1}=10^{-6}$ and the time step is subsequently increased adaptively, to improve efficiency of the procedure while ensuring convergence at each time step. The maximum time step is typically $\Delta t_{\mathrm{max}}\sim 0.1$ which allows calculations to be performed long enough to reach steady-state conditions without the need for exceptional computational resources.

3.3 Calculation of the long-range Coulomb interaction

Since the functions $c_{s}(r,t)$ are discretized in space and time the Coulomb integral in equation 13b can be approximated with a summation technique. In this section the concentration at the mesh nodes is abbreviated as $c_{s}(r_{j},t_{n})=c_{s}(x_{1}^{j_{1}},x_{2}^{j_{2}},x_{3}^{j_{3}},t_{n})=c_{s}^{j,n}$, with the composite index $j=(j_{1},j_{2},j_{3})$, and in the same way $c_{s}(r_{i},t_{n})=c_{s}^{i,n}$. The volume element corresponding to mesh node $j$ is denoted by, $\nu_{j}$. The charge density in the volume, $\nu_{j}$, around mesh node $j$ is turned into a point charge, $q^{j,n}=\sum_{s}z_{s}\,c_{s}^{j,n}\,\nu_{j}$, placed at position $r_{j}$. With this definition the Coulomb integral is converted into the Coulomb sum

[TABLE]

where $\phi(r,t_{n})$ is still a continuous function of space, and it can only be evaluated at the discrete times, $t_{n}$, since the $c_{s}^{j,n}$ are only known at discrete times $t_{n}$. The prime indicates that the possible term $r=r_{j}$ is excluded from the summation. The summation index $j$ runs over all volume elements of the infinite slab, and since the potential decays with $r^{-1}$, the sum is only slowly and conditionally convergent; it cannot be trivially truncated. Therefore, the Parry summation technique [37, 38], which is an Ewald summation technique for 2D periodic systems, is used. In this technique the sum is split into a real and a Fourier (reciprocal) space part, the short-range interactions are evaluated in real space and the long-range interactions are evaluated in Fourier space. The advantage is rapid convergence compared to the direct summation.

The rhombus shown in figure 6 is defined as the repeat unit of the 2D periodic tiling, and the rhombohedral prism shown in figure 6 is a repeat unit of the infinite slab, used to carry out the Parry summation. The symmetry of the hexagonal prism, with distinct hexagonal surfaces on the top and bottom surface, is used to increase the computational efficiency of the evaluation of the Coulomb sum.

During a calculation the charge density changes with time and the potential needs to be updated with the changing charge density; however, the geometry of the hexagonal structure and the mesh do not change. Therefore, the Parry summation is performed only once, at the beginning of the calculations, to determine the Green’s function for the given mesh and periodic cell structure. The Green’s function is calculated for unit charges on the individual mesh nodes, and is stored as a matrix, $g_{ij}$, for the nodes of the irreducible zone, where $g_{ij}$ denotes the potential at node $i$ due to a unit charge and all its images, generated by symmetry and periodicity, at node $j$. During the time-dependent calculations, the potential is calculated by summing the discretized charge density, multiplied by the Green’s function, over the nodes of the irreducible zone,

[TABLE]

This strategy greatly reduces the computational cost of evaluating the Coulomb integral.

A stretched mesh with variable mesh spacings, which are different in the $\hat{e}_{1}$ and $\hat{e}_{3}$ direction, is found to lead to divergence problems when performing the summation over the point charges, $q^{j,n}$. Therefore, “Gaussian smearing” is used for the charge densities at the individual mesh nodes. The local charge density at mesh node $j$ and discrete time $t_{n}$, $\rho(r_{j},t_{n})$, is replaced by a normalized Gaussian distribution

[TABLE]

where $\beta=(2\sigma^{2})^{-1}$ and $\sigma$ is the width of the Gaussian centered around $r_{j}$. With the Gaussian distribution the electrostatic potential contributed by the charge density corresponding to mesh node $j$ takes the form

[TABLE]

In the limit of $|r-r_{j}|\rightarrow 0$ the electrostatic potential contributed by the charge density centred on mesh node $j$ itself becomes

[TABLE]

which has to be added to the potential at this node generated by all the other charges in the system. The width $\sigma$ of the Gaussians is chosen to be that of the smallest mesh spacing in the slab.

3.4 Iterative solver for the self-consistent calculations

A multidimensional Newton method is used to solve the nonlinear system of discretized equations, including the continuity equations for all species on all mesh points, their boundary conditions, and the Coulomb summation to obtain the self-consistent electrostatic potential value at each mesh point. Depending on the number of mesh points and the number of species the system, which needs to be solved in every time step, can contain several thousand variables. Therefore, we have implemented a Jacobian-Free-Newton-Krylov (JFNK) method from the NOX package of the Trilinos Project [39] in our code. JFNK methods are nested iteration methods and can achieve Newton-like convergence without the cost of forming and storing the true Jacobian required for ordinary Newton methods [40].

4 Diffusion of oxygen through an alumina film

The effective charges of the point defects are defined with respect to the ions of the pristine lattice, and the Kröger-Vink notation for point defects is used, although for generality we use integers to denote the charge states rather than the original superfix ∙ or ′ symbols, but retain the × notation for neutral species.

Oxygen membrane permeation or diffusion through an oxide film is a non-equilibrium process, involving oxygen exchange and electron transfer reactions at the oxide surfaces and the transport of defects between the surfaces.

4.1 1D analytic model for the permeation rate

Our derivation of an analytic model for the oxygen permeation rate emphasizes the role of the electric field for the establishment of approximately stoichiometric proportions for the dominant vacancy and electronic species of the surface concentrations, $\mathcal{C}_{s}$, which determine the permeation rate, and the appearance of an effective (ambipolar) diffusion coefficient, see also B, in the presence of the self-consistent electric field.

The concentrations of the defect species on the surfaces exposed to the oxygen gas environment, characterized by the oxygen partial pressure and the temperature, evolve over time and depend on the reactions taking place and the transport to and from the surfaces. At the oxide – gas interface the absorption (desorption) of oxygen and creation (annihilation) of aluminum vacancies can be described by

[TABLE]

For the purpose of a simple analytical treatment we assume a steady-state is reached at the surfaces, in which Schottky and electron-hole equilibrium, see reactions 37 and 38, characterized by the equilibrium constants $K_{\text{S}}$ and $K_{\text{eh}}$, respectively, are attained. The reaction 23 could be formulated with holes instead of electrons, however, assuming instantaneous equilibration between electrons and holes both formulations yield identical results. At high applied oxygen gas partial pressures, aluminum vacancies are formed predominantly at the surface and oxygen vacancies are annihilated, with the reverse scenario applying for low oxygen partial pressures. Since instantaneous Schottky equilibration is assumed at the surfaces no additional reaction is required involving oxygen incorporation by oxygen vacancy annihilation. In reaction 23 for $P_{\text{O}_{2}}^{\text{hi}}$ applied electrons are consumed by the oxygen atoms creating oxygen sublattice sites and aluminum vacancies, and due to the electron hole equilibration the concentration of holes increases simultaneously; therefore, in the limit of $P_{\text{O}_{2}}^{\text{hi}}$ applied, aluminum vacancies and holes are the dominant species on the $P_{\text{O}_{2}}^{\text{hi}}$ surface. Similarly, oxygen vacancies and electrons are expected to dominate on a $P_{\text{O}_{2}}^{\text{lo}}$ surface. The electric field due to the charged defect species modifies the transport to and from a particular surface in such way that local charge neutrality holds approximately between the dominant vacancy species, $\nu$, and the charge compensating electronic species, $\overline{\nu}$, on the surface

[TABLE]

where the $\{\nu,\overline{\nu}\}$ pair denotes $\{\operatorname{V_{Al}^{3-}},\operatorname{h^{+}}\}$ at high oxygen partial pressure, and $\{\operatorname{V_{O}^{2+}},\operatorname{e^{-}}\}$ at low oxygen partial pressure. This approximation does not hold in general for intermediate pressures, or if the vacancies are much more numerous than the electronic defects in equilibrium. Applying the law of mass action to the reaction given in 23 with equilibrium constant $K_{\text{Al}}$, together with Schottky and electron hole equilibrium at the surfaces, the concentration of the dominant vacancy species, $\nu$, is found to adhere to the simple power law given by

[TABLE]

where $z_{\text{O}}$ is the integer charge number of the oxygen ion, $z_{\nu}$ is the charge of the dominant vacancy species, and the power law exponent, $n_{\nu}$, which in our case takes the values, $n_{\operatorname{V_{O}^{2+}}}=-1/6$ and $n_{\operatorname{V_{Al}^{3-}}}=3/16$, and the prefactor, $f_{\nu}$, are introduced. It is noteworthy that the realisation of these power laws in the experiments supports the model of vacancies carrying nominal ionic charges, not the fractional charges that are usually estimated in electronic structure calculations, which typically vary from 1.0 to 1.6 in DFT calculations for oxides [41, 42]. It is also far from obvious that simple point defect diffusion, well understood in bulk crystals, is the mechanism of diffusion in grain boundaries, in which the prefactors of diffusion coefficients are anomalously large [6, 43].

In general, the flux of aluminum ions, $J_{\text{Al}}$ implies the take up or release of $|3J_{\text{Al}}/4|$ molecules of $\text{O}_{2}$ per unit time and unit area of surface, and the opposite flux of oxygen ions leads to $|J_{\text{O}}/2|$ molecules $\text{O}_{2}$ taken up or released. The ionic charges have the opposite sign and the fluxes have the opposite sign due to their different gradients, therefore the permeation rate in molecules of oxygen per unit area per second is given by

[TABLE]

and hence by using $I_{\text{ion}}=I_{\text{vac}}$

[TABLE]

where the current density carried by the ions, $I_{\text{ion}}$, and vacancies, $I_{\text{vac}}$, represent the same physical current density, which has to be balanced by the electron and hole currents.

To interpret and compare with results of our fully time-dependent approach, we introduce a simple one-dimensional, steady-state treatment, the derivation of which is given in B. By considering only the flux of the dominant vacancy species, denoted by suffix $\nu$, and its charge-compensating electronic species, $\overline{\nu}$:

[TABLE]

the permeation rate is approximated by

[TABLE]

Using equation 25 the permeation rate can be written as

[TABLE]

where again $n_{\operatorname{V_{O}^{2+}}}=3/16$, and $n_{\operatorname{V_{Al}^{3-}}}=-1/6$. Equation 31 shows that the permeation rate becomes a power law only for $\left(P_{\text{O}_{2}}^{II}\right)^{n_{\nu}}\gg\left(P_{\text{O}_{2}}^{I}\right)^{n_{\nu}}$. It is also worth noting that the permeation rate is not proportional to the oxygen pressure difference across the membrane, but rather to the difference in the vacancy concentration between the two surfaces.

The oxidation rate “$J$” of aluminum in atoms per second per unit area in terms of the permeation rate of oxygen in molecules per second per unit area is given by

[TABLE]

We can compare the above results with the Wagner theory [8, 9, 44], in which the oxidation rate is given by

[TABLE]

where $z_{\text{Al}}$ is the aluminum ion charge and $\mu_{\text{O}}$ is the oxygen component chemical potential. By assuming that the conductivity of the dominant vacancy species, $\sigma_{\nu}=D_{\nu}c_{\nu}z_{\nu}^{2}e_{0}^{2}/k_{\text{B}}T$, is much larger than that of the other vacancy species, and similarly for the dominant electronic species, $\overline{\nu}$, the oxidation rate simplifies to

[TABLE]

Applying the charge neutrality approximation to this equation and considering relation 32 the permeation rate given in equations 30 follows, which makes our 1D model consistent with the Wagner theory. The derivation of this form of the Wagner theory, however, requires Schottky and electron-hole equilibrium at any position in the film or scale,

[TABLE]

neither of which is enforced in our treatment. Strictly speaking, the above conditions only need to be met by the spatial variations of the chemical potentials, but if equations 35 hold at the surfaces both formulations are equivalent and there is no arbitrary additive constant.

Equations 30 are equivalent to the formulas given by [11] for the permeation rate in the limit of $D_{\nu}\ll D_{\overline{\nu}}$. However, the details of the derivation are different such as the assumption of purely conductive transport of electrons and holes, $I_{s}=-\sigma_{s}\nabla\phi$, which is made in [11].

4.2 Time-dependent 3D calculations

In the calculations within the hexagonal slab model transport of vacancies is simulated for the 3D grain structure and grain boundaries with finite width $\delta$. The time-dependent calculations prove to have a similar steady-state limit to the 1D model described above, in which the divergences of the fluxes and current densities vanish. The vacancy flux is calculated directly from the concentration and electrostatic potential gradients, and the permeation rate follows from equation 28. Before the steady-state is achieved, an average permeation rate can be calculated by numerical integration of the vacancy flux through the grain boundary, which converges faster to the steady-state permeation rate than the permeation rate calculated from the vacancy fluxes on individual mesh nodes.

4.2.1 Reaction equations

Exchange of oxygen between the gas phase and the oxide surface includes multiple steps, namely adsorption, dissociation, surface diffusion, charge transfer, and incorporation into the oxide surface, each of which might be the rate limiting step. In the simulations the assumptions of instantaneous equilibration between the vacancy species inducing the Schottky equilibrium, and the instantaneous equilibration between electrons and holes are relaxed, and the respective processes are formulated with rate dependent reaction equations. In addition to the oxygen incorporation mechanism given in 23, we consider oxygen absorption (desorption) by annihilation (creation) of an oxygen vacancy

[TABLE]

The Schottky reaction is included at the alumina - oxygen gas surfaces

[TABLE]

Electron and hole recombination and generation,

[TABLE]

is considered at the grain boundary and on the surfaces.

In this application the surface reactions, $\mathcal{R}_{s}$, depend on the time-dependent surface concentrations, $\mathcal{C}^{I/II}_{s}(t)$, and the oxygen partial pressures, $P_{\text{O}_{2}}^{I}$ or $P_{\text{O}_{2}}^{II}$, on the respective surfaces, ($I$) or ($II$), and the reaction rate constants, $k_{i,f}$ and $k_{i,b}$

[TABLE]

The law of mass action is used to derive expressions for the reactions $\mathcal{R}_{s}$, and $R_{s}$. The coupled system of reactions, $\mathcal{R}_{s}$, at the alumina - oxygen gas surfaces, is used for the boundary conditions of the transport equations. The equations are provided in A and they describe reactions that may proceed in either direction, depending on the species concentrations and the oxygen gas pressure applied.

4.3 Simulation results

Results are presented here to address some of the questions posed in section 1.2 of the introduction, regarding the membrane experiments [11] and their interpretation with 1D diffusion models based on the Wagner theory. Time dependent calculations are performed and characteristic aspects of the dynamics are highlighted and explained. A permeation calculation reproduces the power laws found in [11].

4.3.1 Initial conditions and choice of parameters

The system at $t=t_{0}$ is assumed to be of strictly stoichiometric composition, $2c_{\operatorname{V_{O}^{2+}}}(t_{0})=3c_{\operatorname{V_{Al}^{3-}}}(t_{0})$, and the set of initial values of the concentrations, $\{c_{s}(t_{0})\}$, and $\{\mathcal{C}_{s}(t_{0})\}$, is assumed to be in equilibrium with oxygen partial pressure $P_{\text{O}_{2}}^{\text{eq}}$, which is defined as the reference pressure so that only ratios $P_{\text{O}_{2}}/P_{\text{O}_{2}}^{\text{eq}}$ enter the equations. The initial concentrations in the surface layer and in the grain boundary are chosen to be equal, $c_{s}=\mathcal{C}_{s}$, and independent of the position in the grain boundary and on the surface. This choice also constrains the electron and hole concentrations, $c_{\operatorname{e^{-}}}(t_{0})=c_{\operatorname{h^{+}}}(t_{0})$, in order for overall charge neutrality to be maintained. The initial values for the concentrations, $c_{s}(t_{0})$, are used to specify the equilibrium constants, $K_{i}$, of the reactions, which are also related to the reactions rates

[TABLE]

where $\mu_{s}^{0}=-k_{\text{B}}T\,\text{ln}(c_{s}(t_{0}))$ is the reference chemical potential with $c_{s}(t_{0})$ in units of $c_{\text{ref}}$, and $w_{i,s}$ is the stoichiometric coefficient of species $s$ in the $i$-th reaction.

This choice of initial parameters only leaves undetermined the ratio of point defect to electronic defect initial concentrations, which is a function of the difference between the vacancy formation (segregation) energy and the Fermi level. For bulk $\alpha$-alumina the ionic and electronic disorder has been analysed from experimental thermodynamic data [45] and more recently from first principles calculations [46]. At temperature $1900\,$K electronic disorder is found dominant for the bulk material. However, we are unaware of data for the surface and interface equilibrium defect concentrations and aim to justify our choice by comparison to the experimental permeation data, see section 4.3.5. We expect the point defects which are favourably formed at the surfaces and interfaces to introduce states in the band gap, and therefore assume much higher vacancy, electron, and hole concentrations at the surfaces, and interfaces than in the bulk material. The reference concentration for the simulations is chosen as $c_{\text{ref}}=10^{18}\,\text{cm}^{-3}$, which would correspond to about $4\times 10^{-5}$ defects per formula unit of bulk $\operatorname{Al_{2}O_{3}}$, and will only be reached in the high or low pressure limits. The concentrations are initialized at $t=t_{0}$ with $c_{\operatorname{h^{+}}}=0.1\,c_{\text{ref}}$, and $c_{\operatorname{V_{O}^{2+}}}=5\times 10^{-3}\,c_{\text{ref}}$.

For $T=1900\,$K and with $\varepsilon_{r}^{\operatorname{Al_{2}O_{3}}}=9.8$, the reference screening length, see equation 14, becomes, $\ell_{\text{D}}=9.4\,$nm. The thickness of the slab is set to, $L=1\,\mu$m, the side length of the hexagon to $a_{\text{hex}}=1\,\mu$m, the grain boundary and surface layer thickness to $\delta=1\,$nm. The diffusion coefficients of the species are chosen as: $D_{\operatorname{h^{+}}}=D_{\operatorname{e^{-}}}=1\,D_{\text{ref}}$, and $D_{\operatorname{V_{O}^{2+}}}=D_{\operatorname{V_{Al}^{3-}}}=0.01\,D_{\text{ref}}$ where $D_{\text{ref}}$ is unknown and used to define the time scale of the simulations. One reaction rate constant in each reaction has to be estimated and they are set to: $k_{\text{O},b}=10^{3}$, $k_{\text{Al},b}=10^{2}$, and $k_{\text{eh},b}=k_{\text{S},b}=10^{3}$.

In sections 4.3.2, 4.3.3 and 4.3.4 calculations are discussed for which the oxygen partial pressure is raised at surface ($II$), the pressure at surface ($I$) is kept constant at $P_{\text{O}_{2}}^{I}=P_{\text{O}_{2}}^{\text{eq}}$, and the pressure ratio is given by $P_{\text{O}_{2}}^{II}/P_{\text{O}_{2}}^{I}=10^{5}$. Snapshots of the evolution at two times, $t_{n}$, and $t_{N}$, with $t_{n}\ll t_{N}$, are shown. The time $t_{n}$ is chosen to capture characteristic behaviour in the initial transient, and $t_{N}$ is the time at which steady-state conditions are achieved.

4.3.2 Evolution of the concentrations

Figure 7 shows snapshots of the concentrations; the coordinate system is defined in figure 4. The pressure $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$ and the $\operatorname{V_{Al}^{3-}}$ and the $\operatorname{h^{+}}$ concentrations are the dominant species. Local charge neutrality holds approximately for the dominant species, $3c_{\operatorname{V_{Al}^{3-}}}\approx c_{\operatorname{h^{+}}}$, throughout most of the grain boundary and at surface ($II$) in the steady-state $t_{N}$ limit but is clearly violated at surface ($I$).

4.3.3 Evolution of the charge density and the electrostatic potential

Figure 8 shows the charge density, $\rho$, and the electrostatic potential, $\phi$, corresponding to the concentrations shown in figure 7. Charge is accumulated near the high pressure surface ($II$) initially until time $t_{n}$, see figures 8 (a) and 9; this change in the local charge density is due to the different magnitudes of the diffusion coefficients of the mobile species and the requirement of charge conservation within the irreducible zone of the structure. If all species had the same diffusion coefficient the charge density would be identically zero, $\rho(t)\equiv 0$, for all times in calculations with the present model. The negatively charged surface ($II$) generates a linearly decreasing potential with increasing $x_{3}$ ($\phi=-|a|x_{3}+b$); however, at $t_{n}$ it is screened to an almost constant value by the charge density that has accumulated near surface ($II$). At time $t_{N}$, see figure 8 (b), some of charge in the grain boundary has propagated to surface ($I$) and the grain boundary has become weakly negatively charged, see also figure 9. However, the electrostatic potential has become almost linear as a function of $x_{3}$, and is dominated by the surface charges; the remaining departure from linearity is due to the charge density within the grain boundary.

The quantity $\kappa$ denotes the scaled reference screening length $\ell_{\text{D}}/L$, defined in equation 14, and characterizes the decay in the local charge density from the surface into the grain boundary. $\kappa^{2D}$ is the equivalent of the scaled reference screening length for two dimensional systems and is defined here as

[TABLE]

The simulations show that $\ell_{\text{D}}^{2D}$ is a better estimate for the spatial extent of the variations in the charge density within the grain boundary than $\ell_{\text{D}}$. It should be pointed out here that simulations in which the grain boundaries and surfaces are idealized as planes without finite thickness $\delta$ would yield the same results, and $c_{\text{ref}}\delta$ in equation 41 would be replaced by a reference concentration per unit area. The only requirement for this to hold is that the total numbers of each species present in the irreducible zone initially, are chosen equal for the simulations with and without finite thickness $\delta$.

4.3.4 Surface and grain boundary charges

The simulations are initialized with zero total charge and charge conservation requires the total charge in the irreducible zone of the structure, $q_{\text{tot}}$, to remain equal to zero. This is satisfied for the calculations performed, see figure 9, which provides a useful check on the numerical accuracy and stability of the solution.

The charge integrated over the grain boundary volume changes with time, $q_{\text{GB}}(t)$, and the boundary carries an excess of electrons in the steady-state limit. This demonstrates that local charge neutrality, which is an assumption of the simple 1D models, is not consistent with this 3D model. Two characteristic time scales for the charging and discharging of the grain boundary are found. The time constants can be estimated by analogy to “$RC$”-circuits, with time constant $\tau=RC$, where $R$ is the resistance, and $C$ is the capacitance. The initial charge built up within the grain boundary is due to the diffusion of the faster of the dominant species, for $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$ the positively charged holes, $\overline{\nu}=\operatorname{h^{+}}$, and the time constant is estimated as follows from appropriate values of $R$ and $C$:

[TABLE]

The time constant characterizing the discharging of the grain boundary is estimated from the diffusion of the slower of the dominant species, for $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$ the aluminum vacancies, $\nu=\operatorname{V_{Al}^{3-}}$:

[TABLE]

Where the angle brackets denote the spatial average and the physically meaningful conductivity is the product $\sigma\delta$. A third time scale characterizes the time elapsed until reaching the steady-state, and is estimated here from the diffusion length using the effective diffusion coefficient defined in equation 30b. For the parameters uses here, $D_{\nu}=D_{\overline{\nu}}/100$ and in the limit of $D_{\nu}\ll D_{\overline{\nu}}$,

[TABLE]

The three time scales are indicated in figure 9.

4.3.5 Membrane permeation calculations

The permeation rate normalized by the thickness of the slab, $PL$, can be calculated from the current densities of the mobile defect species. Figure 10 shows the permeation rate in the steady-state limit as a function of the oxygen pressure ratio $P_{\text{O}_{2}}^{II}/P_{\text{O}_{2}}^{I}$ for a range of $P_{\text{O}_{2}}^{II}$ values, while $P_{\text{O}_{2}}^{I}$ is held fixed. The individual contributions of the $\operatorname{V_{Al}^{3-}}$ and $\operatorname{V_{O}^{2+}}$ are indicated by $PL_{\operatorname{V_{Al}^{3-}}}$ and $PL_{\operatorname{V_{O}^{2+}}}$, respectively. The parameters are the same as those in the previous section.

Comparing the simulated permeation rate with the experimental values [11, 13] for $T=1900\,$K and $P_{\text{O}_{2}}^{II}/P_{\text{O}_{2}}^{I}=10^{5}\,\text{Pa}/1\,\text{Pa}$ the reference diffusion coefficient is calculated, $\delta D_{\text{ref}}=3.4\times 10^{-13}\,\text{m}^{3}\,\text{s}^{-1}$. From the average aluminum vacancy concentration in the grain boundary at $P_{\text{O}_{2}}^{II}/P_{\text{O}_{2}}^{I}=10^{5}$, $\langle c_{\operatorname{V_{Al}^{3-}}}\rangle=10^{17}\,\text{cm}^{-3}$, and the concentration of alumina formula units, $c_{\operatorname{Al_{2}O_{3}}}=2.26\times 10^{22}\,\text{cm}^{-3}$, the aluminum diffusion coefficient is estimated to be $\delta D_{\text{Al}}=\delta D_{\operatorname{V_{Al}^{3-}}}\,\langle c_{\operatorname{V_{Al}^{3-}}}\rangle/2c_{\operatorname{Al_{2}O_{3}}}=7.6\times 10^{-21}\,\text{m}^{3}\,\text{s}^{-1}$, which is close to the value reported in [11], $\delta D_{\text{Al}}=4.5\times 10^{-21}\,\text{m}^{3}\,\text{s}^{-1}$. The agreement indicates that the reference concentration is a reasonable choice provided the assumptions of $c_{\nu}(P_{\text{O}_{2}}^{\text{eq}})\ll c_{\overline{\nu}}(P_{\text{O}_{2}}^{\text{eq}})$ and $D_{\nu}\ll D_{\overline{\nu}}$ are applicable. Indeed, if electrons and holes would diffuse slower than the vacancies, $D_{\nu}\gg D_{\overline{\nu}}$, the permeation rate would be limited by the electronic defects and the diffusion coefficients determined from the permeation experiments [11] would reflect the electronic diffusion coefficient rather than the ionic one.

The variation of the logarithm of the permeation rate $P$, see equation 31, with the oxygen chemical potential at surface ($II$) in the limits of high and low applied oxygen partial pressure, $P_{\text{O}_{2}}^{II}$, yield the power law exponent corresponding to the dominant defect species, $\nu$,

[TABLE]

in the limit of high and low $P_{\text{O}_{2}}^{II}$.

An asymptotic electronic current density, $I_{\text{el}}$, and vacancy current density, $I_{\text{vac}}$ can be calculated similarly to the permeation rate. The current density of vacancies, $I_{\text{vac}}=I_{\operatorname{V_{O}^{2+}}}+I_{\operatorname{V_{Al}^{3-}}}$, is equal and opposite to the electronic current, $I_{\text{el}}=I_{\operatorname{e^{-}}}+I_{\operatorname{h^{+}}}$, at all pressure ratios in the long time limit, see figure 11, this means the net current, $I_{\text{net}}=I_{\text{vac}}+I_{\text{el}}$ is zero.

The plot for the average conductivities, see figure 12, is equivalent to a Kröger-Vink (Brouwer) diagram [47] for the mobile species present, except that the average concentrations are scaled by the corresponding diffusion coefficients. The averaged conductivities of the dominant point defect species, $\nu$, adhere to

[TABLE]

The crossover in the dominant point defect species in the permeation rate and the average conductivity is observed at the same value of $P_{\text{O}_{2}}^{II}$.

4.3.6 The Schottky equilibrium

As discussed in section 4.1, Schottky equilibrium and electron-hole equilibrium are required conditions for the validity of the Wagner model and equation 33. In this section we examine to what extent these equilibria are attained in our grain boundary calculations as the steady-state limit is approached.

Figure 13 shows an example of the voltages, $\eta_{s}/z_{s}e_{0}$, in the long time limit, the component chemical potentials are shown in figure 14. The calculations are performed for $P_{\text{O}_{2}}^{I}=P_{\text{O}_{2}}^{\text{lo}}$, and $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$. In these calculations Schottky equilibrium does not hold, see figures 13, and 14, except at the surfaces where the Schottky reaction is included in the equations. The reaction rates at the surfaces are higher than the transport between the surfaces and the grain boundaries, therefore Schottky equilibrium is attained at the surfaces. The defect chemical potentials are calculated from the ideal dilute solution approximation, $\mu_{s}=k_{\text{B}}T\,\text{ln}(c_{s}/c_{s}(t_{0}))$, where $c_{s}(t_{0})$ are the equilibrium concentrations. The spatial variation of the electrochemical, and chemical potentials would be significantly different if internal Schottky equilibrium was imposed.

The component chemical potential distributions $\mu_{\text{Al}}(x_{3})$, and $\mu_{\text{O}}(x_{3})$ shown in figure 14 are different from those shown in figure 8 of reference [13] for polycrystalline alumina membrane permeation experiments. Schottky equilibrium together with electron-hole equilibrium implies $3\mu_{\text{O}}+2\mu_{\text{Al}}=0$, which is not satisfied within the grain boundary, see figure 14, because $3\mu_{\operatorname{V_{O}^{2+}}}+2\mu_{\operatorname{V_{Al}^{3-}}}\neq 0$. The qualitative discrepancy in the chemical potential distributions between reference [13] and this treatment is likely due to the assumption of Schottky equilibrium in [13].

A Schottky reaction term could be added to the transport equations 9, which would lead to Schottky equilibrium in the grain boundaries depending on the reaction and transport rates. However, the internal Schottky reaction with formation and dissolution of oxide at the grain boundaries would induce stress and cause a non-trivial modification in the permeation rate. The simulations discussed in this work are thought to correspond to the limit in which stress at the grain boundaries prohibits the formation and dissolution of oxide internally, similarly to its effect in the bulk material.

5 Discussion

The power laws for the oxygen permeation rate found experimentally in the limits of high and low applied oxygen pressures, are confirmed in the calculations with the slab model in the steady-state limit. The experiments are performed on $0.25\,$mm thick polycrystalline alumina membranes while the geometry of the slab in the calculations with $L=a_{\text{hex}}=1\,\mu$m is chosen to resemble more closely the situation of planar films growing with a columnar grain structure. The variation of the permeation rate, and the average conductivity with the oxygen chemical potential are related to the power law exponent $n_{\nu}$, see figures 10 and 12, and equations 48 and 49. The power law exponent in turn depends on the stoichiometry of the quasi-chemical reactions at the surfaces, see equation 26. The applied pressure at which the transition between $p$-type and $n$-type ionic conductivity of the grain boundaries takes place depends on the ratio of the vacancy diffusion coefficients in the grain boundary, $D_{\operatorname{V_{O}^{2+}}}/D_{\operatorname{V_{Al}^{3-}}}$, given that $\operatorname{V_{O}^{2+}}$, and $\operatorname{V_{Al}^{3-}}$ are the mobile vacancy species.

The time-dependent calculations not only elucidate the initial transient behaviour but also help to clarify the steady-state. Varying the applied oxygen partial pressure from $P_{\text{O}_{2}}^{\text{eq}}$ to $P_{\text{O}_{2}}^{\text{hi}}$ or $P_{\text{O}_{2}}^{\text{lo}}$ on one of the surfaces changes the rate of creation and annihilation of vacancies, electrons and holes in stoichiometric proportions on the surface, the resulting chemical potential gradients between the surfaces along the grain boundaries drive the transport of species through the slab. The resulting fluxes of the mobile species depend on their diffusion coefficients and will therefore not necessarily preserve the stoichiometric proportions of the surface populations of the species, except when all species have the same diffusion coefficient. The local charge density thereby becomes non-zero, and generates an electric field that retards the diffusion of the fastest dominant species, but enhances the diffusion of the slower dominant species which has an effective charge of opposite sign.

Three time scales involved in the dynamics of species concentration evolution are identified, $\tau_{c}$ (see equation 44) characterizes the rate of charge build up on the surface and in the grain boundary, $\tau_{d}$ (see equation 46) characterizes the discharging of the grain boundary, and $\tau_{s}$ (see equation 47) characterizes the time to reach the steady-state at which $\nabla\cdot J_{s}\simeq 0$.

During the initial transient, $t\lesssim\tau_{c}$, the different diffusion coefficients lead to the accumulation of charge in the grain boundary and on surface ($II$). If, for example in the $P_{\text{O}_{2}}^{\text{hi}}$ case, holes and aluminum vacancies are generated on surface ($II$), the holes, which are assumed to be the faster species, $D_{\operatorname{V_{Al}^{3-}}}\ll D_{\operatorname{h^{+}}}$, diffuse away from the surface, generating positive charge density in the grain boundary, and surface ($II$) becomes negatively charged due to the aluminum vacancies. Once the holes reach surface ($I$) it becomes positively charged, and as the slower aluminum vacancies diffuse into the grain boundary it is discharged over a time $t\sim\tau_{d}$. In this example the electric field effectively enhances the aluminum vacancy diffusion, and retards the diffusion of the oxygen vacancies and holes. In the steady-state limit the electrostatic potential difference between the surfaces, $\Delta\phi=\phi^{II}-\phi^{I}$, is on the order of the thermal voltage, $k_{\text{B}}T/e_{0}$, in the high and low pressure limits considered.

For the calculation shown in figure 13 with $P_{\text{O}_{2}}^{II}=P_{\text{O}_{2}}^{\text{hi}}$ and $P_{\text{O}_{2}}^{I}=P_{\text{O}_{2}}^{\text{lo}}$ the voltage between the surfaces observed in the simulations is $V_{\operatorname{e^{-}}}=-(\eta_{\operatorname{e^{-}}}^{II}-\eta_{\operatorname{e^{-}}}^{I})/e_{0}=50\,\text{mV}$. $V_{\operatorname{e^{-}}}$ depends strongly on the species diffusion coefficients.

Physically, the mobile species distributions reconfigure to screen the electric field arising due to non-zero local charge density, and the Debye screening length (see equation 15) characterizes the spatial extent of the variations in the local charge density. Due to screening effects the local charge density and electrostatic potential are challenging to resolve numerically near the surfaces. For spatial variations in the charge density within the grain boundary the two-dimensional equivalent of the reference screening length, see equation 41, is found to be appropriate. Apart from the initial transient the electrostatic potential in the simulations is dominated by the contributions from the surface charges. In the parameter regimes investigated the non-zero charge density in the grain boundary does not affect the defect fluxes significantly, the 3D calculations can therefore be mapped onto a 1D model. In the intermediate pressure regime the dominance of one vacancy species is less pronounced and in particular for $c_{\nu}(r,t)\sim c_{\nu}(t_{0})$ in the steady-state the local charge neutrality approximation becomes invalid. For small species concentrations the Debye length gets larger, and the electric field gets weaker. Both of the later facts are considered responsible for the failure of the analytic model in the intermediate pressure regime, see figure 10.

The simulations with the fixed oxygen chemical potential difference between the surfaces reach a stationary non-equilibrium state with a constant rate of entropy production. This means the populations of vacancies, electrons, and holes reach a dynamic equilibrium state, in which they are created and annihilated at the same rate and $\frac{\partial c_{s}}{\partial t}\simeq 0$, while the oxygen chemical potential difference sustains non vanishing fluxes between the surfaces.

6 Conclusions

A model for time-dependent grain boundary diffusion of ions and electrons through a film of polycrystalline oxide has been constructed, in the form of a slab comprising hexagonal columnar grains. The long-range Coulomb interactions between the grain boundary planes affect the mass transfer dynamics through the slab significantly only during the initial transient. The electric field generated by the evolution of the charge density influences the transport significantly; in the long time limit it is dominated by the surface charges.

Four mobile defects, charged aluminum and oxygen vacancies, electrons and holes, were considered, and simulations were made to compare with the behaviour observed in alumina oxygen permeation experiments. The power laws for the permeation rate in the alumina membrane calculations depend on the stoichiometry of the quasi-chemical reactions at the membrane surfaces.

Work is in progress to extend the model; for example by including a Schottky reaction in the grain boundary we can couple the fluxes to the development of internal stress.

We have introduced a simplified, one-dimensional analytic model that employs the approximation of a single dominant defect, which appears from the 3D calculation to be justified. It does not assume internal Schottky equilibrium which is often assumed in Wagnerian treatments of oxidation. The analytic model is shown to agree well with the simulation results under certain limiting conditions, but fails in the intermediate oxygen partial pressure regime.

7 Acknowledgements

M.P.T was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the Engineering and Physical Sciences Research Council under grant EP/G036888/1. The authors would like to acknowledge the funding and technical support from BP through the BP International Centre for Advanced Materials (BP-ICAM) which made this research possible. We are grateful for discussions with Profs. Arthur Heuer, Matthew Foulkes and Brian Gleeson, and Dr. Paul Tangney. The authors acknowledge support of the Thomas Young Centre under grant TYC-101.