Multilayer interface tracking model of pure tungsten oxidation

Under high oxygen partial pressure conditions, the WO₃ layer grows by the transfer of mass from the metal represented by the reaction:

$$\begin{align} {\mathrm{W}}^{6+} + 3{\mathrm{O}}^{2-}\rightarrow {\mathrm{WO}}_3. \end{align} \tag{ 1 }$$

This reaction is enabled by the reduction of molecular oxygen ⁵ on the environment side:

$$\begin{align} {\mathrm{O}}_2+2\textrm{e}^{-}\rightarrow 2{\mathrm{O}}^{2-} \end{align} \tag{ 2 }$$

which receives electrons from the oxidation of metal tungsten within the metal substrate:

$$\begin{align} {\mathrm{W}}\rightarrow {\mathrm{W}}^{6+}+6\textrm{e}^{-}. \end{align} \tag{ 3 }$$

The formation of WO$_{3-x}$ phases follows the same sequence as above but with other oxidation states of tungsten as the starting point, such as +2, +3, +4, or +5.

On the basis of these considerations, in the most general case the oxidation of metallic tungsten will evolve into a structure characterized by five distinct layers. These layers represent different tungsten–oxygen structures ranging from WO₃ on the environment side (where O radicals are produced) to the W–O metallic solid solutions (where the metal is reduced). In between, three more oxide phases are expected to form, namely, WO_2.9, WO_2.72, and WO₂. More details about the structure and properties of the different oxide compounds can be found in appendix A.1. As shown in figure 1(a), WO_2.9 only exists above 484 ^∘C, while WO_2.72 is expected to appear only above 585 ^∘C. WO₃ and WO₂ exist as equilibrium phases in the entire temperature range. All these oxides appear as strong line compounds in the phase diagram, suggesting that they should form with near perfect stoichiometry.

Zoom In Zoom Out Reset image size

It is worth remarking that the redox reactions described above are not rate-limiting, as they are governed by reaction rates that occur over times that are orders of magnitude shorter than diffusive timescales. We demonstrate this quantitatively in appendix A.2.

The model starts with the adsorption of oxygen radicals at the surface of the material. Adsorption is a thermally activated process that sets the rate of oxygen atom ingress in the solid:

$$\begin{align} \dot{c}_0 = c_0\nu_{\mathrm{ads}} = c_0\nu^0_{\mathrm{ads}}\exp\left(-\frac{E_{\mathrm{ads}}}{kT}\right) \end{align} \tag{ 4 }$$

where c₀ is the concentration of oxygen radicals in solution on the environment side, $\nu^0_{\mathrm{ads}}$ is an oxygen pick-up rate, $E_{\mathrm{ads}}$ is the adsorption activation energy, k is Boltzmann's constant, and T the absolute temperature. The value of c₀ is environment-dependent, but as a reference the oxygen concentration in air under standard ambient conditions is $5.64\times10^{24}$ m⁻³ [29]. The activation barrier for adsorption depends on the structure of the surface, and ranges between 0.1 and 1.0 eV [30–32]. Here we take $\nu^0_{\mathrm{ads}} = 10^{12}$ Hz and $E_{\mathrm{ads}} = 1.0$ eV. Note that the model can be solved using both $\dot{c}_0$ or c₀ as boundary conditions on the environmental side.

The model of oxide growth can then be formulated by considering the motion of four interfaces, each characterized by its own evolution equation. Together, the motions of these interfaces determine the rate of oxidation of the underlying metal substrate. This means that the oxide grows into the metal, which is the correct physical assumption due to mass conservation. As well, here we ignore the relative volumetric expansion of the oxide layers relative to the initial metal phase (the implications of this assumption are discussed further at the end of appendix C.1). The evolution equations can be written assuming a sharp-interface model (known as a Stefan model) in one dimension:

$$\begin{align} \dot{s}_i(x) = \frac{\Delta J_i|_{s_i}}{n_i\rho_{\mathrm{W}}} \end{align} \tag{ 5 }$$

where the subindex i refers to a specific interface, x is the depth coordinate, $\Delta J_i|_{s_i}$ is the net flux of (free) oxygen exiting the interface towards the interior, n_i is the O-to-W atomic ratio, and $\rho_{\mathrm{W}}$ is the atomic density of tungsten. The above equation describes the rate of advance of a specific interface on the basis of the oxygen flux differential on both sides of the interface. In essence, an interface must move to maintain mass continuity dynamically. The full derivation of equation (5) is contained in appendix B. To determine $\Delta J_i|_{s_i}$, one must solve a diffusion equation in each oxide sublayer to first determine the free oxygen concentration profile and then obtain the flux using Fick's first law:

$$\begin{align} \Delta J_i|_{s_i}& = J_i|_{s_i+}-J_i|_{s_i-} = D_i \frac{\partial c}{\partial x}\Big|_{s_i+}-D_{i-1}\frac{\partial c}{\partial x}\Big|_{s_i-} \end{align} \tag{ 6 }$$

$$\begin{align} \frac{\partial c}{\partial t} & = D_i \frac{\partial^2 c}{\partial x^2}. \end{align} \tag{ 7 }$$

Equation (7) is subjected to two concentration boundary conditions, one for each interface bounding each oxide sublayer. At the inner edge of the layer (i.e. at $x = s_i$), c = 0 always. At the outer edge of the layer ($x = s_{i-1}$), $c = \rho_{\mathrm{W}}\left(n_{i-1}-n_i\right)$. The exceptions to this are the outermost and innermost layers (which are both fixed, i.e. $\dot{s}_0 = \dot{s}_5 = 0$), where $c = c_0$ and $J_5 = 0$, respectively. The specific values for each sublayer are given in table 1.

Table 1. Details for the formulation of the Stefan model of multilayer interface evolution.

Layer	Compound	n	Outer interface	Inner interface	Outer boundary condition $(\rho_{\mathrm{W}})$	Inner boundary condition $(\rho_{\mathrm{W}})$
1	WO3	3.0	s0	s1	$c = c_0$	c = 0
2	WO2.90	2.90	s1	s2	c = 0.1	c = 0
3	WO2.72	2.72	s2	s3	c = 0.18	c = 0
4	WO2	2.0	s3	s4	c = 0.72	c = 0
5	W–O metal	—	s4	s5	c = 2.0	J = 0

Note that here we are assuming that all the oxygen above the stoichiometric concentration in each layer is insoluble in that layer, while according to the phase diagram part of it will start to precipitate as the next phase richer in oxygen. However, here we treat this as having a minor effect on the kinetics. As well, a delayed formation of a WO_x phase can take place even when the oxygen concentration is lower than $x\rho_{\mathrm{W}}$. In that case, as equation (4) indicates, building up the oxygen concentration to the required value for each stoichiometric oxide phase to form takes an amount of time given by $\Delta t = x/\dot{c}_0$. The assumption made in this work is that this time does not compute towards the total oxide layer formation time. The equations to be solved for each layer, with the appropriate numerical notation and boundary conditions to be used are given below. To indicate the boundary conditions on the right-hand side of a given interface i (with respect to the frame of reference on figure 2), the suffix '$s_i+$' ('$s_i-$' to refer to the left side).

Zoom In Zoom Out Reset image size

Layer 1:

Material	WO3
Governing equations	$\begin{array} {lcl} \dot{s}_1(x)& = &\frac{\Delta J_1|_{s_1}}{3.0\rho_{\mathrm{W}}} \\[6pt] \Delta J_1|_{s_1}& = & D_2 \frac{\partial c}{\partial x}\Big|_{s_1+}-D_1 \frac{\partial c}{\partial x}\Big|_{s_1-} \\[6pt] \frac{\partial c}{\partial t} & = & D_1 \frac{\partial^2 c}{\partial x^2} \\ \\ \end{array}$
Boundary conditions	$\begin{array} {lcl} \\ c(0) = c_0-3\rho_{\mathrm{W}},~c(s_1) = 0 \\[6pt] (\textrm{as~well}:~\dot{c}(0) = \dot{c}_0) \end{array}$

Layer 2:

Material	WO2.9
Governing equations	$\begin{array} {lcl} \dot{s}_2(x)& = &\frac{\Delta J_2|_{s_2}}{2.9\rho_{\mathrm{W}}} \\[6pt] \Delta J_2|_{s_2}& = & D_3 \frac{\partial c}{\partial x}\Big|_{s_2+}-D_2 \frac{\partial c}{\partial x}\Big|_{s_2-} \\[6pt] \frac{\partial c}{\partial t} & = & D_2 \frac{\partial^2 c}{\partial x^2} \\ \\ \end{array}$
Boundary conditions	$\begin{array} {lcl} c(s_1) = 0.1\rho_{\mathrm{W}} ,~c(s_2) = 0 \end{array}$

Layer 3:

Material	WO2.72
Governing equations	$\begin{array} {lcl} \dot{s}_3(x)& = &\frac{\Delta J_3|_{s_2}}{2.72\rho_{\mathrm{W}}} \\[6pt] \Delta J_3|_{s_3}& = & D_4 \frac{\partial c}{\partial x}\Big|_{s_3+}-D_3 \frac{\partial c}{\partial x}\Big|_{s_3-} \\[6pt] \frac{\partial c}{\partial t} & = & D_3 \frac{\partial^2 c}{\partial x^2} \\ \\ \end{array}$
Boundary conditions	$\begin{array} {lcl} c(s_2) = 0.18\rho_{\mathrm{W}} ,~c(s_3) = 0 \end{array}$

Layer 4:

Material	WO2
Governing equations	$\begin{array} {lcl} \dot{s}_4(x)& = &\frac{\Delta J_4|_{s_3}}{2.0\rho_{\mathrm{W}}} \\[6pt] \Delta J_4|_{s_4}& = & D_5 \frac{\partial c}{\partial x}\Big|_{s_4+}-D_4 \frac{\partial c}{\partial x}\Big|_{s_4-} \\[6pt] \frac{\partial c}{\partial t} & = & D_4 \frac{\partial^2 c}{\partial x^2} \\ \\ \end{array}$
Boundary conditions	$\begin{array} {lcl} c(s_3) = 0.72\rho_{\mathrm{W}} ,~c(s_4) = 0 \end{array}$

Layer 5:

Material	W–O (metal)
Governing equations	$\begin{array} {lcl} \frac{\partial c}{\partial t} & = & D_5 \frac{\partial^2 c}{\partial x^2} \\ \\ \end{array}$
Boundary conditions	$\begin{array} {lcl} c(s_4) = 2.0\rho_{\mathrm{W}} ,~J(s_5) = 0 \end{array}$

The above model builds on similar approaches devised to study diffusion in multilaminates.

When a given layer cracks, a modification to the above formulation must be adopted. When a layer i bounded by oxygen concentrations of $c_{i-1}$ and c_i fragments after it reaches a critical width, $w^\ast_i$, it opens up percolation pathways that equalize the free oxygen concentration across it [33], i.e. $c_{i-1} = c_i$. When such a situation arises, the layer in question stops growing and simply acts as an oxygen bridge between its two adjacent layers. The exception to this behavior is WO₃, which grows unabated even after fragmenting. In this work, while a provision is introduced in the model to switch to this mode of oxygen transport once $w^\ast_i$ is reached, the results that will be presented below pertain to growth times shorter than such threshold.

Here we adopt a first-order finite difference scheme to solve the diffusion equations (Fick's 1st and 2nd laws) given in the previous sections. A one-dimensional central finite difference approach results in the following discretized form:

$$\begin{align} c^{n+1}_i & = c^n_i + \delta tD\left(\frac{c^n_{i+1}-2c^n_i+c^n_{i-1}}{\delta x^2}\right) \end{align} \tag{ 8 }$$

$$\begin{align} J^{n}_i& = -D\frac{c^n_{i+1}-c^n_{i-1}}{2\delta x}\approx -D\frac{c^n_{i}-c^n_{i-1}}{\delta x} \end{align} \tag{ 9 }$$

where $\delta t$ and $\delta x$ are the time and space increments, and the subindices n and i refer to temporal and spatial steps. Both equations above are subjected to the Neumann stability condition $\delta t \lt \delta x^2/D$.

Application of the boundary conditions leads to:

$$\begin{align*} c^0_i = 0,~c^n_0 = c_0, ~c^n_{s_5} = c^n_{s_5 -\delta x}. \end{align*}$$

While the total oxygen concentration (stoichiometric plus free) is continuous at each interface, each diffusion equation in each layer is solved only for the free oxygen concentration. This results in jumps in the concentration depending on whether one looks at the interface from the left or from the right in figure 2. With the coordinate x arriving at each interface, one has:

$$\begin{align*} c^n_{s_1} = c^n_{s_2} = c^n_{s_3} = c^n_{s_4} = 0 \end{align*}$$

and with x leaving each interface:

$$\begin{align*} c^n_{s_1} = 0.10,~c^n_{s_2} = 0.18,~c^n_{s_3} = 0.72,~c^n_{s_4} = 2.00 \end{align*}$$

all in units of $\rho_{\mathrm{W}}$.

The above model is fully defined when the external parameters, T and c₀, and the material properties, $\rho_{\mathrm{W}}$ and the diffusivities, are specified. Here, T and c₀ are prescribed, and the value of $\rho_{\mathrm{W}}$ is known. That leaves the diffusivities D_i ($i = 1,5$) as the remaining material constants to determine.

The general expression for the diffusion coefficient is the standard Arrhenius form:

$$\begin{align} D(T) = D_0\exp\left(-\frac{E_\mathrm{a}}{kT}\right) \end{align} \tag{ 10 }$$

where D₀ is the diffusion pre-factor and $E_\mathrm{a}$ is the activation energy. That means that five separate pre-factors and five activation energies must be determined for each oxide layer. In all cases, transmission electron microscopy (TEM) experiments have revealed that the oxide always grows into the tungsten metallic matrix, confirming that oxide growth is controlled by oxygen anion diffusion from the environment side inwards [34].

In any case, the first thing to ascertain is whether oxygen is transported by a vacancy mechanism (through exchanges with empty oxygen sublattice sites) or as interstitial atoms in an otherwise fully-occupied lattice. On first inspection of figure 1(a), the line compound structure of all the oxide sublayers considered here (in figure 2) suggests that oxygen atoms will diffuse as free interstitials in all cases. However, there is sufficient consensus about the abundance of oxygen vacancies, both thermal and chemical (i.e. due to natural stoichiometric deviations) in all oxide phases [35–38], and so here we assume going forward that all oxygen transport within the oxide layers takes place via a vacancy mechanism. Next, we select oxygen diffusion properties for each layer shown in section 2.3, either from literature sources or—if unavailable—by carrying out DFT calculations using the methodology described in appendix C.1.

2.5.1. WO₃ layer.

2.5.2. WO₂ layer.

In contrast to WO₃, no data are available for oxygen diffusion in WO₂ and we proceed to calculate O diffusivity in layer 4. Here too, a vacancy mechanism is assumed. The starting point for the calculations is an expanded expression for the diffusivity:

$$\begin{align} D_4(T) = z c_{\mathrm{V}}(T)\ell^2 \nu(T) \end{align} \tag{ 11 }$$

where z is a geometric factor, $\ell$ is the jump distance,

$$\begin{align*} c_{\mathrm{V}}(T) = \exp\left(-\frac{E^{\mathrm{V}}_\mathrm{f}}{kT}\right) \end{align*}$$

is the thermal vacancy concentration ($E^{\mathrm{V}}_\mathrm{f}$ is the oxygen vacancy formation energy), and

$$\begin{align*} \nu(T) = \nu_0\exp\left(-\frac{E^{\mathrm{V}}_\mathrm{m}}{kT}\right) \end{align*}$$

is the jump frequency with $E^{\mathrm{V}}_\mathrm{m}$ the oxygen vacancy migration energy. Equations (10) and (11) imply that $D_0 = z\ell^2 \nu_0$ and $E_\mathrm{a} = E^{\mathrm{V}}_\mathrm{m}+E^{\mathrm{V}}_\mathrm{f}$. The factor z, which includes the number of jump sites n_p and the dimension of diffusion d, is ${n_p}\big/ {2d} = {5}\big/{3}$, and the jump distance $\ell$ is 2.76 Å for the oxygen sub-lattice in WO₂ [41, 42]. In this work, $E^{\mathrm{V}}_\mathrm{f}$, $E^{\mathrm{V}}_\mathrm{m}$, and ν₀ were obtained using DFT calculations as reported in appendix C.2. The final expression for $D_4(T)$ is given in equation (C.6).

No data whatsoever exist for the diffusivity of oxygen in the WO_2.9 and WO_2.72 layers. In their analysis of oxidation kinetics data, Nagy and Humphry-Baker identified three bands of oxygen partial pressure, low (0.1 atm), medium (0.2 atm), and high (3.4–13.6 atm). They extracted effective activation energies of 1.22, 1.40, and 1.94 eV for each one [26]. By mapping O partial pressures to oxygen/metal ratio, these can loosely be ascribed to diffusion in WO_2.72, WO_2.9, and WO₃, respectively. Indeed, these energies are not inconsistent with the values of 1.3 and 1.5 eV measured and calculated for oxygen diffusion in WO₃ [17, 36]. We thus select 1.22 and 1.40 eV as the activation energies for WO_2.72, WO_2.9. In terms of the prefactor, for lack of a more suitable option, here we take the one in D₁ to be the same as in D₂ and D₃.

2.6.1. Metallic W.

The solubility of oxygen in W is very low [15] and O atoms are found primarily in interstitial positions. Specifically, interstitial oxygen is seen to diffuse in the W lattice throughout tetrahedral positions with a migration energy of 0.15 eV [43]. The diffusion pre-factor has been obtained by Alkhamees et al [44] to be $7.6\times 10^{-8}$ m$^2\,$s⁻¹, i.e.:

$$\begin{align} D_5(T) = 7.6\times 10^{-8}\exp\left(-\frac{0.15}{kT}\right)~~\left({\mathrm{m}}^2\,{\mathrm{s}}^{-1}\right). \end{align} \tag{ 12 }$$

Table 2 lists all the diffusivities used in the multilayer model.

Table 2. Compilation of diffusivity parameters used in the multilayer tracking model developed in this work.

Diffusivity	Material	D0 (m$^2\,$s−1)	$E^{\mathrm{V}}_\mathrm{m}$ (eV)
D1	WO3	$6.8\times10^{-6}$	1.3
D2	WO2.9		1.40
D3	WO2.72		1.22
D4	WO2	$3.1\times 10^{-5}$	1.90
D5	W metal	$7.6\times 10^{-8}$	0.15

As indicated at the end of section 2.3, fragmentation of an oxide scale occurs due to the large epitaxial misfit between the metal substrate and the oxide phase. These creates residual stresses (known as Pilling–Bedworth stresses) that put the oxide under compression, building up as the layer grows. At a critical thickness threshold, which we term $w^\ast$, the layer fragments generally creating a network of cracks through which oxygen can flow unimpeded [33]. At this stage, oxygen transport ceases to be the rate limiting step in oxide layer growth, which is subsequently seen to scale linearly with time. Linear growth can also coexist with parabolic (diffusive) growth. Indeed, the presence of exponents between 0.5 and 1.0 in oxidation experiments is suggestive of the coexistence of cracked layers and protective layers [19, 26].

By its own definition, these stresses set in only between the oxide phase immediately in contact with the metallic base. As such, in a multilayer setting as the present one, one would expect only one single layer to undergo cracking at any given time. A compilation of data relating temperature, critical width $w^\ast$, and oxide phase is provided in table A1 (appendix A.2). These are based on measurements performed by Gulbransen and Andrew [16] and analyzed by Nagy and Humphry-Baker [26]. As the table shows, the WO₂ phase is common to all oxidation conditions, consistent with the expectation that it should be the layer immediately adjacent to the metal substrate ⁶ . On the basis of these data, in the forthcoming simulations, we assume that in impoverished oxygen atmospheres, a WO₂ layer forms and grows up to a temperature-dependent thickness given by the listed values of $w^\ast$.

Next, we present results of the oxide layer growth as a function of temperature. Based on the qualitative observations by Gulbransen et al (table A1), where broadly three oxidation regimes were observed as a function of T, here we consider the following three temperature regions: (a) at or below 600 ^∘C (873 K), (b) between 600 ^∘C and 800 ^∘C (873–1073 K), and (c) at or above 800 ^∘C. We initialize the model with the temperature value and a boundary oxygen concentration, c₀, that favors the formation of the oxide phases listed in table A1. We set $\delta t$ and $\delta x$ to 1 ns and 1 nm, respectively, which are sufficiently small to always guarantee the numerical stability of the finite-difference solution procedure. In all cases, we show results for the first 30 s of evolution, when the thickness of all oxide layers involved is well below the fragmentation threshold (second column of table A1). As such, the results presented below correspond to the regime when the oxide maintains its protective structure, which can thus be assumed until its thickness reaches the critical width at each temperature.

We give the results in each case as a pair of plots, one showing the depth profiles of all the sublayers involved color-coded by the time instant for which they are displayed, and the other showing the position of the different interfaces as a function of time. c is always given in units of $\rho_{\mathrm{W}}$. The depth profiles are always color-coded by the scale used in figure 1(b) to facilitate the identification of each oxide type.

3.1.1. Kinetic model simulations at 600 ^∘C (873 K).

The results at this temperature are shown in figure 3. As indicated in table A1, only WO_2.72 and WO₂ form at this temperature. Figure 3(a) reveals a very thin WO₂ layer behind WO_2.72, also captured in figure 3(a). A power law fit to the s₄-t graph yields $s_4(t) = 3.747t^{0.407}$, indicating a substantial deviation from pure parabolic growth law.

Zoom In Zoom Out Reset image size

3.1.2. Kinetic model simulations at 600 $\lt$ T $\lt$ 800 ^∘C.

The slow growth of the WO₂ layer carries over to the simulations at 650 ^∘C and 700 ^∘C (923 and 973 K). At 650 ^∘C, however, all oxide sublayers are present, figures 4(a) and (b), with WO₃ also displaying slow growth kinetics. At 700 ^∘C, figures 4(c) and (d), the evolution is qualitative similar to that at 600 ^∘C given in the previous subsection. In this case, the growth law exponents at 650 ^∘C and 700 ^∘C are 0.41 and 0.44.

Zoom In Zoom Out Reset image size

3.1.3. Kinetic model simulations at and above 800 ^∘C.

Figure 5 shows the corresponding graphs at 800 ^∘C and 1000 ^∘C. At these temperatures, again according to the observed phases in the experiments listed in table A1, only the WO₂ phase forms. Interestingly, at 800 ^∘C the growth kinetics is quite slow, with the power law fit being $s_4(t) = 1.852t^{0.454}$ (in microns). By contrast, at 1000 ^∘C the temperature is sufficiently high to overcome the diffusion barrier of oxygen in WO₂, giving rise to faster kinetics compared to the previous scenarios. As such, the fitted expression in this case is $s_4(t) = 7.812t^{0.490}$. As revealed by the value of the power exponents, having one layer only forming part of the oxide significantly approximates the kinetics to the pure parabolic case.

Zoom In Zoom Out Reset image size

All the fitting coefficients are compiled for each temperature in table 3. It is worth mentioning that the power law exponents obtained here are similar to those measured in 'smart' W-based alloys [10, 45].

Table 3. Power law fitting coefficients $s_i(t) = at^b$ for all the temperatures considered in this work. Sublayer thickness given in units of microns.

T (K)	Number of sublayers	a	b	$k_\mathrm{p}$ (mg2 cm−4 min−1)
873	2	3.747	0.407	4.28
923	4	7.930	0.413	20.31
973	2	6.435	0.441	16.69
1073	1	1.852	0.454	3.85
1273	1	7.812	0.490	34.20

The model presented here is formulated as a Stefan interface tracking approach where different oxide/oxide and oxide/metal interfaces evolve in time driven by discontinuities in oxygen fluxes across each one of them. These fluxes emerge in each oxide sublayer from the diffusion of free oxygen (i.e. not chemically immobilized) as a function of temperature. As such, it is crucial to have accurate oxygen diffusivities for the different oxide layers. These can be inferred from mass gain experiments, by assuming parabolic growth laws governed by pure diffusion. Unfortunately, most experimental works do not discriminate among the different oxide subphases, which has prompted the use of computational techniques to calculate the diffusivity in specific WO compounds [39, 46]. Even then, the computational data is scarce and does not address many of the gaps in our understanding of oxygen diffusion in tungsten oxide. For example, it is not clear whether oxygen diffusion primarily takes place by an interstitial or a vacancy mechanism in each of the different oxide phases. Here, we have applied electronic structure calculations to quantify the migration parameters of oxygen atoms in WO₂ by a vacancy mechanism, which was lacking in the literature.

From this point of view, this work is a first attempt at embedding the thermodynamic and kinetic complexity of tungsten oxide growth into a comprehensive mesoscale kinetic model that captures the more important features of the oxidation process. Below we provide further discussion about how the model has been designed, implemented, and validated.

In this work, the simulated oxide structures are inspired by empirical observations of the outermost oxide layer during environmental exposure of pure W [16, 17, 19, 47] (see appendix A.2). Of note is the apparent absence of tungsten trioxide in oxidized samples. Indeed, WO₃ is not customarily seen to form as a standalone layer in the majority of experimental works scrutinized in this paper. Instead, WO₃ forms only after lower order oxides have grown and cracked, allowing for fast oxygen transport pathways that stimulate the formation of the trioxide. For example, in their recent work, Cifuentes et al [19] note that at 600 ^∘C the oxidation was initially controlled by the development of a continuous WO_2.72 layer, which grows until it starts cracking leading to linear kinetics. At 700 ^∘C and 800 ^∘C, after cracking of the WO_2.72 layer, WO_2.92 is progressively seen to form, eventually turning into WO₃ in the course of the oxidation. Other times, WO_2.92 volatilization after cracking exposes tungsten to the exposed surface, leading to the formation of a yellow external WO₃ layer [19]. This evidence points to a mechanism of formation of WO₃ by which its formation takes place only after the WO_2.72 (and WO₂, discussed further below) layer has fragmented and lost its protective abilities. Another possibility worth mentioning is based on experiments performed in synthetic air (O₂ partial pressures above 20 kPa) where WO₃ was seen to form as a standalone layer [9, 48].

The lack of stated observation of WO₂ in many of these experiments may be due to its intrinsic small thickness. Indeed, our calculations suggest that oxygen transport through tungsten dioxide is slow (mostly owing to a migration energy of almost 2.0  eV, table 2), making WO₂ a true protective layer. Only above 1000 ^∘C is WO₂ seen to form at high rates, although at such temperatures all oxide layers start becoming structurally unstable due to sublimation. On this point, it must be noted that sublimation is not a feature of the present model, although studies have consistently shown that it becomes the dominant effect at high temperatures [49, 50].

Based on this presumed behavior, our interpretation of the kinetics of oxide growth is that: (a) the free oxygen concentration under nominal conditions is generally never sufficiently high to form WO₃ directly; (b) the most common oxide sublayers are WO_2.72 and WO₂; (c) these layers will grow in a protective fashion until a critical thickness is reached, at which point cracking occurs and the oxygen-rich phases WO_2.92 and WO₃ begin to form; (d) whenever multiple oxide sublayers form in succession, deviations from the standard parabolic growth law are clearly seen. In our case, fracture and fragmentation is not a part of the current model, and we control the formation of the different oxides simply by regulating the amount of available oxygen, c₀, on the environment side of the exposed metal. c₀ is in fact the only parameter of our model that is 'adjusted' to reproduce certain observations.

A useful exercise is to compare the growth kinetics simulated in this work to available experimental data. As is customary in most works, oxide scale growth is measured in terms of tungsten mass gain assuming parabolic growth. Consequently, the mass gain constant is obtained from first to the data expressed linearly in time

$$\begin{align*} (\Delta m)^2 = \left[(n\rho_{\mathrm{W}}) s_4(t)\right]^2 = k_\mathrm{p}t \end{align*}$$

where $\Delta m$ is the mass gain, $n\rho_{\mathrm{W}}$ is the amount of oxygen gained when a WO_n oxide forms (expressed as a mass density), and $k_\mathrm{p}$ is the parabolic kinetic constant. Next, we fit all the sublayer evolutions shown in figures 3–5 to ideal parabolic growth expressions and plot the simulated values of $k_\mathrm{p}$ in figure 6 along with a set of experimental measurements. The numerical values of $k_\mathrm{p}$ extracted from the simulations are also given in table 3.

Zoom In Zoom Out Reset image size

The results show that, while the simulated values generally overestimate the experimental measurements, they lie within the broader range of the data at some temperatures. However, it is worth noting that the experimental results are not all representative of the same conditions. Indeed, broadly speaking, they can be divided into two categories: short and long exposures to the oxidizing medium. For example, in the experiments by Sikka et al [17], the samples were exposed only for several hundred seconds, while in all the other cases [16, 19, 47], they were exposed to oxygen for several hours. Given the slow diffusivity of oxygen across most of the oxide phases, hundreds of seconds may not be a sufficiently long time to reach visible thickness. Using one micron as the minimum observable oxide thickness (reasonable based on the data in table A1), the diffusive time scales for oxygen transport in the outer layers (WO₃, WO_2.9, and WO_2.72), and in WO₂ at 850 K are respectively 7.5 s and 1.7 h, approximately (estimated using the diffusivities reported in table 2). Therefore, within the short timescales simulated here (≈30 s) and those in the experiments by Sikka et al (≈200 s), the oxide layers may not have fully reached a steady growth rate. The exact diffusive timescale will vary depending on the temperature and therefore not all the experimental data may correspond to the same kinetic state. In any case, given the figurative 'distance' between how the model was developed and how experiments are performed, we consider these differences reasonable.

As discussed above, during oxidation of W, several different oxides may form. WO₃ is the most stable oxide, and at ambient temperature and pressure it displays a monoclinic structure (γ-WO₃) consisting of a network of WO₆ octahedra. Above 330 ^∘C, WO₃ is seen to transform into an orthorhombic structure (β-WO₃), and tetragonal above 740 ^∘C (α-WO₃) [51]. WO₂ has a monoclinic structure in ambient conditions, and an orthorhombic structure at high temperatures and pressures. The crystal structure of monoclinic WO₂ also consists of WO₆ octahedra but with the difference that the octahedra in WO₂ are edge-shared while they are corner-shared in WO₃ [46]. Between WO₃ and WO₂, the stoichiometric oxides W₁₈O₄₉ (WO_2.9), W₂₄O₆₈ (WO_2.72) have been reported at temperatures of 485 ^∘C and above [27, 52].

These intermediate phases W₁₈O₄₉ and W₂₄O₆₈ form due to the thermal stabilization of distorted arrangements of the WO₆ octahedra [53–55]. As the amount of chemical oxygen decreases, the usual corner-sharing arrangement of octahedra is partially replaced by groups of edge-sharing octahedra, which form pockets of shear planes. This is due to the tendency of the lattice to eliminate single oxygen vacancies by a crystallographic shear process [54, 55]. When these shears become coordinated, new stoichiometric lattices can form, e.g. WO_2.9 and WO_2.72. However, as revealed by the W–O phase diagram shown in figure 1(a) [27], these phases have very narrow stability ranges close to the theoretical stoichiometry values. Interestingly, each phase is characterized by a distinctive color, as shown in figure 1(b), associated with the optical diffraction of the different octahedral arrangements.

To extract the critical scale thicknesses at which a transition from parabolic (protective layer) to linear (loss of protective layer) growth occurs, we resort to the tabulated data given by Gulbransen et al [16] in a series of oxidation experiments at several temperatures and different oxygen atmospheres. The data are listed in table A1 (units from the original paper have been updated to more current usage). The presumed phase is inferred from the observed color and the oxygen concentration (high O concentrations favoring WO₃ and low ones favoring WO₂). The oxide scale thickness at which a deviation from parabolic growth is observed is marked as the value of $w^\ast$ to be used in the simulations.

Table A1. Transition oxide layer thickness calculated as a function of temperature at an O₂ partial pressure of 0.1 atm (adapted from [16]). The observed color and presumed oxide phase are also included. The presumed phase is inferred from the observed color and the oxygen concentration (high O concentrations favoring WO₃ and low ones favoring WO₂). The mass density of metallic W is 19.3 g cm⁻³, which amounts to an atomic concentration of $6.3\times10^{28}$ m⁻³.

T (∘C)	$w^\ast$ (µm)	Observed surface color	Presumed oxide
600	52	Blue black	WO2.72/WO2
625	77	Black with yellow spots	WO3/WO2
650	100	Blue black with yellow covering	WO3/WO2.72/WO2
700	155	Blue black	WO2.72/WO2
750	181	Black	WO2
800	207	Black	WO2
850	310	Black	WO2
900	310	Blue black	WO2.72/WO2
950	340	Black	WO2
1000	360	Black	WO2
1050	360	Black	WO2
1100	390	Yellow black	WO3/WO2

It is worth mentioning that the present model disregards the reaction time between W and O. This is justified by the reaction rate being much faster than the diffusion rates. This can be quantitatively demonstrated by considering typical rates for the W–O reaction, e.g., as given by the expression:

$$\begin{align*} K_{\textrm{W-O}} = 6.1\times10^{23} \exp\left(-2.1/kT\right) \end{align*}$$

which represents the number of reactions per unit area per unit time in (cm$^{-2}\,$s⁻¹) [56] for an oxygen partial pressure of 49 Torr. At a temperature of 1500 K, this gives $5.5\times10^{16}$ reactions per squared centimeter per second. When multiplied times $\rho_{\mathrm{W}}$ (assuming an O-rich environment), one gets approximately ${\sim}10^{38}$ molecules of oxide (cm⁻² s⁻¹). Compare this with diffusion-controlled oxygen fluxes obtained as $J_{\mathrm{O}} = -D (\partial c/\partial x)$ with a diffusivity of $6.1\times10^{-5}$ (cm$^{2}\,$s⁻¹) at 1500 K in WO₂ (table 2), which for a generic oxygen atomic concentration gradient of $\partial c/\partial x = 3.1\times10^{26}$ cm⁻⁴ (e.g. figure 5) yields fluxes of ${\sim}10^{21}$ atoms (cm⁻² s⁻¹). As the numbers show, the rate of arrival of oxygen to an interface is orders of magnitude slower than the reaction rate of oxygen with tungsten atoms. This substantiates the notion that oxygen transport is convincingly the rate limiting step of the oxidation process.

The use of the DFT+U technique is necessary to reach a correct description of the electronic configuration of tungsten, which is critical to calculate formation energies and migration energy barriers [57]. Tungsten is a transition metal with electron configuration [Xe]4f¹⁴5d⁴6s², i.e. a partially filled d orbital that requires special DFT treatments [58].

When applying a one-electron method with an orbital-independent potential to transition metal compounds, one has as a result a partially filled d band with metallic type electronic structure and itinerant d electrons. This is definitely the wrong model for late-transition-metal oxides and rare-earth metal compounds where d (f) electrons are well localized and there is a sizable energy separation between occupied and unoccupied sub-bands (the lower Hubbard band and upper Hubbard band in a model Hamiltonian approach [59]). Local density approximation (LDA) and generalized gradient approximation (GGA) functionals describe poorly the electronic properties of these localized orbitals and that is where DFT+U technique comes in. DFT+U accounts for electron–electron interactions in strongly correlated materials because the self-interaction error is dominant in d and f states.

DFT+U is based on Hubbard model, which solves the energy of electron correlation introduced or Coulomb repulsion between electrons at the same atomic orbitals [60]. The Hubbard model explains the transition between the insulation and conduction features of the systems with strong on-site repulsion using the following Hamiltonian:

$$\begin{align} H_\textrm{Hub} = t\sum_{\langle i,j \rangle,~\sigma} c^{\dagger}_{i, \sigma} c_{j, \sigma} + U \sum_{i} {n}_{i \uparrow} {n}_{i\downarrow}. \end{align} \tag{ C.1 }$$

The Hubbard Hamiltonian has two components, namely, the hopping integral t and the electron repulsion strength U. $c^{\dagger}_{i, \sigma}$ and $c_{j, \sigma}$ are the creation and annihilation operators of electron state at sites i and j and spin σ, respectively. The sum runs over all $\langle i, j \rangle$ nearest electron pairs, i and j. ${n}_{i} = {n}_{i \uparrow}+{n}_{i\downarrow}$ is the total density operator at site i.

The Hubbard model allows us to write the total energy of the system as a standard DFT energy plus a term representing the strong correlation of the electronic states in d (and f) orbitals

$$\begin{align} E_{\mathrm{DFT}} + U[n] = E_{\mathrm{DFT}}[n] + E_{\mathrm{Hub}}[n_{i}^{\sigma}] - E_{\mathrm{dc}}[n_{i}^{\sigma}] \end{align} \tag{ C.2 }$$

where $U[n]$ is the correction induced by the DFT+U method, $E_{\mathrm{Hub}}[n_{i}^{\sigma}]$ is the Hubbard functional that describes the correlated states, and $n_{i}^{\sigma}$ is the occupation number of the localized orbitals with a spin of σ at site i. The term of $E_{\mathrm{dc}}[n_{i}^{\sigma}]$ is a 'double counting' term, which is included because when we additively append the Hubbard term $E_{\mathrm{Hub}}[n_{i}^{\sigma}]$, the energy contribution of the related orbitals has already been counted in the DFT term. In our calculations we use a Coulomb interaction parameter of LDAUU = 1.82 eV for pure metal W [61] because it does not affect the distribution of electrons for the attempt frequency calculations. However, we have also repeated the calculations with a value of LDAUU = 6.2 eV used for WO₃ [62] and found only a change of 5% in the migration energy barriers.

Tungsten is a transition metal with partially-filled d orbitals that requires special DFT treatments. Here we use the DFT+U technique with the parameterization for W used by Feng [61]. The GGA and Perdew–Burke–Ernzerhof exchange-correlation potentials were applied to build the pseudopotential of the system and the frozen-core convention was implemented with the projector-augmented wave method. The DFT simulations are performed with a cutoff energy of 600 eV and a k-point sampling of 2 × 2 × 4. The convergence criterion for the energy relaxation of the atomic positions was set to be 10⁻⁵ eV. Defect energies were calculated in a WO₂ supercell containing 36 atoms under constant volume.

C.2.1. Formation energy, $E^{\mathrm{V}}_\mathrm{f}.$

The calculations are based on supercells of W_N O$_{2N}$ lattices, where N is the number of W atoms. A vacancy is created by removing an O atom from a lattice position and placing it at infinity, i.e.:

$$\begin{align} E^{\mathrm{V}}_\mathrm{f} = E_{{\mathrm{W}}_N{\mathrm{O}}_{2N -1}} + \mu _{\mathrm{O}} - E_{{\mathrm{W}}_N{\mathrm{O}}_{2N}} , \end{align} \tag{ C.3 }$$

where $E_{{\mathrm{W}}_N{\mathrm{O}}_{2N-1}}$ and $E_{{\mathrm{W}}_N{\mathrm{O}}_{2N}}$ are the energy of the supercells with and without a vacancy, and $\mu _{\mathrm{O}}$ is the chemical potential of an isolated oxygen atom. Under low oxygen partial pressures, $\mu _{\mathrm{O}}$ satisfies:

$$\begin{align} \mu_{{\mathrm{O}}} = \mu_{{\mathrm{O}}_2} + \frac{1}{2} E^{\,f}_{{\mathrm{WO}}_2} \end{align} \tag{ C.4 }$$

where $\mu_{O_2}$ is the chemical potential of bimolecular oxygen and $E^{f}_{{\mathrm{WO}}_2}$ is the formation energy of tungsten dioxide. For its part, $E^{f}_{{\mathrm{WO}}_2}$ is given by:

$$\begin{align*} E^{\,f}_{{\mathrm{WO}}_2} = \mu_{{\mathrm{WO}}_2} - \mu_{{\mathrm{O}}_2} - \mu_{\mathrm{W}} \end{align*}$$

where the r.h.s. of the equation contains the chemical potentials (energy per molecule) of the perfect WO₂ crystal, bimolecular oxygen, and pure (metal) W. After carrying out the pertinent DFT calculations, we arrive at values of $\mu_{\mathrm{O}} = -7.65$ eV and $E^{\mathrm{V}}_\mathrm{f} = 2.49$ eV. It must be kept in mind, however, that measurable concentrations of vacancies have been detected in tungsten oxide under ambient conditions [63, 64], which suggests that they exist in natural concentrations so as to naturally enable oxygen transport by a vacancy mechanism. This will the assumed mode of oxygen migration in WO₂.

C.2.2. Migration energy: $E^{\mathrm{V}}_\mathrm{m}.$

Migration energies are obtained using the nudged elastic band (NEB) method [43]. In the WO₂ crystal structure (monoclinic), there are three non-equivalent oxygen lattice positions, shown in figure C1(a) and labeled as V$_{{\mathrm{O}}_1}$, V$_{{\mathrm{O}}_2}$, and V$_{{\mathrm{O}}_3}$. The NEB paths for the V$_{{\mathrm{O}}_1}{\longleftrightarrow}$V$_{{\mathrm{O}}_2}$ and V$_{{\mathrm{O}}_2}{\longleftrightarrow}$V$_{{\mathrm{O}}_3}$ migration trajectories are shown in figure C1. The resulting energy barriers are summarized in the following table.

Path	$E^{\mathrm{V}}_\mathrm{m}$ (eV)
V$_{{\mathrm{O}}_1}{\longrightarrow}$V$_{{\mathrm{O}}_2}$	2.19
V$_{{\mathrm{O}}_2}{\longrightarrow}$V$_{{\mathrm{O}}_1}$	1.86
V$_{{\mathrm{O}}_2}{\longleftrightarrow}$V$_{{\mathrm{O}}_3}$	1.96

Zoom In Zoom Out Reset image size

By way of comparison, an energy of 2.46 eV has been calculated for the V$_{{\mathrm{O}}_1}{\longrightarrow}$V$_{{\mathrm{O}}_2}$ path [46]. While all of these oxygen vacancy exchanges are viable diffusive transitions, here we take the V$_{{\mathrm{O}}_1}{\longrightarrow}$V$_{{\mathrm{O}}_2}$ as representative of all vacancy jumps and proceed to calculate the attempt frequency for that path.

C.2.3. Attempt frequency: $\nu_0$ $.$

$\nu_0$ can be obtained by resorting to harmonic transition state theory as [65, 66]:

$$\begin{align} \nu_0 = \frac{ \prod_i^{3N}\omega_i^{\mathrm{eq}} }{ \prod_i^{3N-1}\omega_i^{\mathrm{sp}} } \end{align} \tag{ C.5 }$$

where $\omega_i^{\mathrm{eq}}$ and $\omega_i^{\mathrm{sp}}$ are, respectively, the normal frequencies of vibration of the system at the equilibrium position (e.g. V$_{{\mathrm{O}}_1}$) and saddle point of the migration path. In each case, the normal frequencies can be obtained by diagonalizing the dynamical matrix, D , generated at the equilibrium and saddle point positions, which amounts to solving the eigenvalue problem defined by:

$$\begin{align*} \boldsymbol{D}\boldsymbol{\epsilon} = \boldsymbol{\lambda}\boldsymbol{\epsilon} \end{align*}$$

where ε is the eigenvector matrix and λ is a diagonal matrix containing all the eigenvalues of the system from which the eigenfrequencies can be calculated as $\omega_k = \sqrt{\lambda_{kk}}$. At the saddle point, one of the normal modes of vibration is undefined (imaginary eigenvalue) and eliminated from the product in the denominator of equation (C.5), hence the $3N-1$ limit of the product.

The results for the V$_{{\mathrm{O}}_1}{\longrightarrow}$V$_{{\mathrm{O}}_2}$ path give a value of $\nu_0 = 243.2$ THz, which, when inserted in equation (11), results in $D_0 = 3.087\times 10^{-5}$ m$^2\,$s$^{-1}.$ With this, the final expression for the oxygen diffusivity in WO₂ becomes:

$$\begin{align} D_4(T) = 3.087\times 10^{-5}\exp\left(-\frac{1.9}{kT}\right)~~\left({\mathrm{m}}^2\,{\mathrm{s}}^{-1}\right). \end{align} \tag{ C.6 }$$

These results are added to table 2 to be used in the kinetic model simulations. To close this section, it is worth returning to the point mentioned in section 2.2 about neglecting the differences in specific volume between the oxide and the metal. All the oxide phases considered in this work form and grow under a state of compressive stress. As such, the oxygen diffusivities should in principle be calculated under uniaxial compressive deformation. While experimentally-measured diffusivities are intrinsically obtained under such conditions, this is not the case for the DFT calculations in WO₂ presented here.