A Turing instability in the solid state: void lattices in irradiated metals

TL;DR

This paper demonstrates that Turing instabilities can explain the formation of void superlattices in irradiated metals, using a minimal model and phase field simulations, revealing a new pattern formation mechanism in solids.

Contribution

It introduces a simple, generic model showing Turing patterns can form in solids, specifically explaining void lattice structures in irradiated metals, expanding the understanding of pattern formation in solid materials.

Findings

Void superlattices can be explained by Turing mechanism.

Analytical results are supported by phase field simulations.

The model applies broadly to pattern formation in solid state systems.

Abstract

Turing (or double-diffusive) instabilities describe pattern formation in reaction-diffusion systems, and were proposed in 1952 as a potential mechanism behind pattern formation in nature, such as leopard spots and zebra stripes. Because the mechanism requires the reacting species to have significantly different diffusion rates, only a few liquid phase chemical reaction systems exhibiting the phenomenon have been discovered. In solids the situation is markedly different, since species such as impurities or other defects typically have diffusivities $\propto\!\exp\left( -E/k_{\rm B} T\right)$, where $E$ is the migration barrier and $T$ is the temperature. This often leads to diffusion rates differing by several orders of magnitude. Here we use a simple, minimal model to show that an important class of emergent patterns in solids, namely void superlattices in irradiated metals, could also be explained by the Turing mechanism. Analytical results are confirmed by phase field simulations. The model (Cahn-Hilliard equations for interstitial and vacancy concentrations, coupled by creation and annihilation terms) is generic, and the mechanism could also be responsible for the patterns and structure observed in many solid state systems.

Figures

Equations

$$\dot{s}$$

$$\dot{v}$$

$$F={\int }_V\left[s^2{\left(1-s\right)}^2+\frac{{\gamma }_s}{2}|\nabla s{|}^2+v^2{\left(1-v\right)}^2+\frac{{\gamma }_v}{2}|\nabla v{|}^2\right]\mathrm{d}V.$$

$$\dot{s}$$

$$\dot{v}$$

$$s(\mathbf{x})=\bar{s}+S(\mathbf{x});v(\mathbf{x})=\bar{v}+V(\mathbf{x});c-a\bar{s}\bar{v}=0.$$

$$\bm{A_{q}}=\left(\begin{array}[]{cc}-D_{s}(q^{2}g_{s}+\gamma_{s}q^{4})-a\bar{v}&-a\bar{s}\\ -a\bar{v}&-D_{v}(q^{2}g_{v}+\gamma_{v}q^{4})-a\bar{s}\\ \end{array}\right),$$

$$\det {\mathbf{A}}_{\mathbf{q}}=D_v{D}_s{q}^8+D_v{D}_s(g_s+g_v)q^6+(D_v{D}_s{g}_s{g}_v+a(D_v\bar{v}+D_s\bar{s}))q^4+a(D_v\bar{v}{g}_v+D_s\bar{s}{g}_s)q^2,$$

Full text

Introduction — Patterns formed by Turing instabilities Turing (1952) arise in reaction-diffusion systems due to the competition between diffusion and nonlinear reaction terms. Counterintuitively, a uniform solution for reactant concentrations (known as a base state), stable in the absence of diffusion, can become unstable to the emergence of patterns and ordering once diffusion is switched on. This runs counter to the standard picture of diffusion as a smoothing influence, and is interesting to study from a non-equilibrium physics point of view. Some time after Turing’s original prediction, chemical systems were discovered that exhibited the effect, though they remain rare since the Turing model typically requires the reacting species to diffuse at significantly different rates – unusual in liquid phase chemical systems Cross and Greenside (2009). In the solid state, however, different species’ diffusion rates generically differ by many orders of magnitude, since they are usually governed by nonlinear Arrhenius escape rates $\propto\!\exp\left(-E/k_{\rm B}T\right)$, where the migration barrier $E$ can vary from fractions-of to several eV. We note that crowdion defects in body-centred-cubic (bcc) metals have migration barriers too low for the Arrhenius formula to apply, and their diffusion rates are linear in temperature Fitzgerald and Nguyen-Manh (2008); Swinburne et al. (2013).

An intriguing and technologically important example of solid state pattern formation is void and gas bubble superlattice formation in irradiated metals. First observed in the 1970s Evans (1971); Johnson and Mazey (1978); Sikka and Moteff (1972), the voids generated by the agglomeration of the radiation-induced vacancies can form an ordered superlattice under certain conditions. This runs counter to the more intuitive picture of Ostwald ripening, where large voids grow at the expense of smaller ones. Also, noble gases formed in fission reactors (e.g. Kr, Xe) generally have very low solubility in metals, and hence segregate to regions of high tensile strain. At grain boundaries, this leads to embrittlement, and accelerated mechanical failure. Engineering a stable bubble lattice (formed of voids filled with gas atoms) potentially offers a way to sequester this gas atoms safely away from grain boundaries and extend the life of reactor materials Harrison et al. (2017). Superlattices are most often observed within a temperature window of 0.2-0.4 of the melting point Robinson et al. (2017), and often mimic the lattice symmetry of the underlying crystal, though with a spacing tens or hundreds of times larger; see Ghoniem et al. (2001) for a thorough review. These lattices form over minutes and hours, meaning molecular dynamics simulations cannot hope to directly capture the processes at work.

Various competing mechanisms for superlattice formation have been proposed, including elastic interactions between voids, isomorphic decomposition, phase instability, interstitial dislocation loop punching and anisotropic interstitial diffusion Krishan (1982); Tewary and Bullough (1972); Woo and Frank (1985); Khachaturyan and Airapetyan (1974); Evans (1985); Dubinko et al. (1986). Here we propose an alternative mechanism, and argue that void lattices could emerge as a Turing instability, where diffusion itself destabilizes the uniform base states which solve the steady-state, diffusionless equations of motion. Whilst some or all of the mechanisms above may play a role in the details of the superlattice formation, we show all that is actually required is a region in which local vacancy and interstitial concentration, generation, and annihilation rates satisfy a specific relation, and vacancy and interstitial diffusion rates that are sufficiently different. Ours is the simplest possible model that can capture the diffusion of two reacting species, with like species tending to cluster. It is a gross idealizeation, and neglects many important features of real crystal systems, in particular the anisotropic nature of self-interstitial diffusion and the elastic interactions between species. Nevertheless, it is sufficient to predict the formation and lengthscale of ordered patterns, as we show below. Our purpose here is to present a minimal and general model, which may be applied to many different systems, rather than to focus on the details of specific materials. A systematic study dealing with particular metals and radiation conditions will be published elsewhere.

In the next section, we apply Turing’s linearized analysis to the pair of coupled equations governing the diffusing defects, and extract analytical conditions for the system to support a superlattice of a given wavenumber. We then perform fully non-linear phase field simulations to investigate the system behaviour at longer times, confirming that the superlattice wavenumber predicted by the linear analysis is indeed realised in the full system.

The model — In what follows, $v(\bm{x},t)$ and $s(\bm{x},t)$ denote the concentrations of vacancies and self-interstitials respectively. A phase field model Moelans et al. (2008); Chen (2002) for their evolution leads to Cahn-Hilliard equations Cahn and Hilliard (1958), with additional terms corresponding to creation ($c$) and annihilation ($-asv$, according to the law of mass action):

[TABLE]

The terms in brackets are functional derivatives of the following simple double-well free energy $F[s,v]$ with respect to $s$ and $v$:

[TABLE]

The quartic bulk free energy terms have minima when the concentrations $s$ and $v$ are 0 or 1, encouraging the formation of voids and clusters. The $D$s are the diffusivities, with $D_{s}\gg D_{v}$ in metals, and the $\gamma$s are proportional to the square of the effective interface size between solid and void/cluster regions. We stress that all these parameters take effective values. Since superlattice formation takes place on a timescale of hours, the underlying atomic processes will be averaged over many realizations. For example, the annihilation rate $a$ does not represent the probability of mutual annihilation when a vacancy and self-interstitial atom meet, but rather the fraction of defects which annihilate over a representative region in a representative time interval.

The explicit form of the equations is

[TABLE]

These equations conserve the number of defects during evolution (apart from the explicit creation and annihilation terms), in contrast with the coupled rate equation model Bullough et al. (1975) explored in ref.Ghoniem et al. (2001), which involves only two spatial derivatives. Note that the defects do not interact until they meet and react: this is not a Fokker-Planck model of diffusion in a position-dependent potential, but rather a reaction-diffusion one.

Analytical results — We now follow the analysis due to Turing, and linearize the system about a so-called base state $\bar{s},\,\bar{v}$ which satisfies the static equations, Eqs.(3) with all spatial and temporal derivatives set to zero:

[TABLE]

This imposes a relation between the uniform base states and the creation and annihilation rates. Seeking solutions of the form $(S,V)\equiv\bm{S}=\bm{S_{0}}\exp\left[\lambda t+i\bm{q}\cdot\bm{x}\right]$ leads to the eigenvalue equation $\lambda\bm{S_{0}}=\bm{A_{q}S_{0}}$, with

[TABLE]

where $q=|\bm{q}|=(q_{x}^{2}+q_{y}^{2})^{1/2}$ in 2D, $g_{s}=2(6\bar{s}(\bar{s}-1)+1),$ and $g_{v}=2(6\bar{v}(\bar{v}-1)+1)$. The eigenvalues $\lambda$ are given by the two solutions to det$(\bm{A_{q}}-\lambda\bm{I})=0$. If both solutions for $\lambda(q)$ are negative, the solution decays in time, and hence the base state is stable to perturbations of wavenumber $\bm{q}$. A Turing instability arises when a base state is stable for $D_{s}=0=D_{v}$ (equivalently $q=0$), but becomes unstable when it is perturbed by a certain wavenumber $\bm{q}$. The growing solution then leads to periodic patterns with wavenumber $\bm{q}$.

When $q=0$, $\lambda=0$ or $-a(\bar{s}+\bar{v}),$ so for all base states, no unstable ($\lambda>0$) pattern-forming mode is possible without diffusion. When diffusion is switched on, one or both eigenvalues are pushed above zero when either ${\rm tr}\bm{A_{q}}>0$ and $({\rm tr}\bm{A_{q}})^{2}-4\,{\rm det}\bm{A_{q}}>0$, or ${\rm tr}\bm{A_{q}}<0$ and ${\rm det}\bm{A_{q}}<0.$ A sufficient (but not necessary) condition for the Turing instability is hence ${\rm det}\bm{A_{q}}<0.$ Assuming $\gamma_{s}=\gamma_{v}$, and working in units where $\gamma_{s}=\gamma_{v}=1$ leads to

[TABLE]

a quartic in $q^{2}$, passing through $q^{2}=0$ (reflecting the conservation of vacancies and interstitials). Positive values of $q^{2}$ that lead to a negative value of det$\bm{A_{q}}$ correspond to a pattern with wavenumber $\bm{q}$. $q>2\pi$ is not physically interesting, since it corresponds to patterns of wavelength less than the interface width. Also, $q\to 0$ corresponds to complete decomposition into void and undefected crystal, thus the most predictive, and hence physically interesting, case is the third in Fig.1 (inset), where only a certain range of wavenumbers lead to instability.

Since the equation for the determinant is effectively a cubic, it can be solved analytically, and the value of the superlattice spacing $\Lambda$ can be extracted as a function of the input parameters. For case 3, this is given by $\Lambda=2\pi/\sqrt{Q_{*}/2}$, where $Q_{*}$ is the largest root of ${\rm d}({\rm det}\bm{A_{q}})/{\rm d}Q=0$ (see Fig. 1).

$D_{s}\gg D_{v}$ means that the interstitials generated during a cascade diffuse away faster than the vacancies, typically leading to a “halo” of interstitials surrounding a region of high vacancy density. Setting $\bar{v}=0.25$ and $\bar{s}=0.01$ to reflect an example of this results in the third scenario described above. The determinant is shown in Fig. 1 for several values of the diffusivity ratio ($a=1$ in this plot. The values of $a$ and $c$ are constrained by Eq.(4).

When $D_{s}=D_{v}$, the determinant barely dips below zero, but as the ratio $D_{s}/D_{v}$ increases up to the value of 1000 typical for bcc metals, the instability deepens. The minimum, most unstable, wavenumber $q$ for these parameters is approximately $\sqrt{0.15}$, corresponding to a pattern period $2\pi/(q/\sqrt{2})$ of about 23 times the interface width, or around 100 spacings of the underlying crystal lattice, if we take the interface to be 4 crystal lattice spacings in width (again, this is an effective quantity, chosen to appropriately balance the bulk and interface terms in the free energy, and need not correspond precisely to the size of the physical interface at the void surface). This is consistent with experimentally observed void lattices.

The above values for $\bar{v}$ and $\bar{s}$ represent a reasonable example, but in any irradiated crystal, different regions will have different values. The conditions for instability are not particularly restrictive, however. According to Descartes’ rule of signs, a cubic has two positive roots (i.e. case 3 discussed above) when there are two sign changes between the successive terms in Eq. (6), and the discriminant is positive. Since the first coefficient is always positive, this means the last coefficient must be positive, and at least one of the second and third coefficients must be negative. Fig.2 shows the fraction of the region in parameter space defined by $(\bar{s},\bar{v})\in[0,0.5]\times[0,0.5]$ that satisfies these conditions, and hence supports a pattern-forming instability, as a function of the diffusivity ratio (we restricted the full $(u,v)\in[0,1]$ range to exclude unrealistic base states with $>50\%$ vacancies or interstitials). Several values of the annihilation parameter $a$ are shown. For each value of $a$, the unstable region reaches a plateau when the diffusivity ratio exceeds approximately 100. The lattice-forming region also grows as the effective annihilation parameter falls. For $a<0.1$, around a quarter of the possible values for $\bar{s}$ and $\bar{v}$ lead to case 3 and hence an instability. This condition is again sufficient, but not necessary.

Phase field simulations — The Turing analysis is based on linearization, and it is reasonable to ask whether the patterns remain once the nonlinearity becomes important, and the nascent regions of high vacancy concentration grow into voids. We used the open source Multiphysics Object Oriented Simulation Environment Tonks et al. (2012); Schwen et al. (2017) to integrate Eqs.(3) numerically on a 2D domain, using the finite element method with implicit time integration, starting from an initial condition randomized about $\bar{s}=0.25,\bar{v}=0.01$ with $a=0.5$. The results are shown in Fig. 3. Ordering is absent when $D_{s}/D_{v}<5$ and clearly emerges when $D_{s}/D_{v}\geq 10$. The superlattice spacing is approximately 25 units, confirming that the system selects the fastest-growing unstable mode as predicted by the analytical model. The lattice is hexagonal, which is the expected symmetry that minimizes the free energy for a given wavenumber; the equivalent in three dimensions is body-centred-cubic (bcc) Shoji et al. (2007).

Fig.4 shows the average void area and number of voids against time. Initially, voids nucleate in the regions where the fluctuating initial vacancy concentration is high. For $D_{s}/D_{v}<5$, the standard picture of Ostwald ripening emerges, with large voids growing at the expense of smaller ones. As $D_{s}/D_{v}$ is increased however, the number of voids stabilizes, and an ordered lattice emerges, as is clear from Fig.3. We simulated the system under a variety of different initial conditions, including pre-existing populations of voids of different sizes and distributions, and several different values for the creation term. In all cases with $D_{s}/D_{v}\gg 1$, we found a stable void lattice (see Supplementary Material). The voids do not nucleate in an ordered pattern, and the lattice begins to form after nucleation. Smaller voids on the lattice grow and larger voids shrink, and those not on lattice sites shrink until they disappear. Intriguingly, we observed diffusion-driven migration of established voids to lattice locations, consistent with experimental observations Ghoniem et al. (2001). This occurs in the absence of any advective term in the governing equations (1, 3), and is purely due to the preferential diffusion of vacancies and interstitials so as to form the superlattice. This provides a mechanism for the fast migration of fairly large voids, which might intuitively be expected to be immobile.

Discussion — We have shown that the simplest possible model for diffusing populations of vacancies and interstitials, subject to uniform creation and annihilation, supports void superlattice formation, even in the absence of refinements such as anisotropic interstitial diffusion and elastic interactions. The mechanism responsible for the ordering is the well-known Turing instability. This also offers a possible explanation for the observed temperature window for superlattice formation: the mechanism requires the diffusivities of the vacancies and interstitials to differ significantly. The ratio $D_{s}/D_{v}\propto\exp(-(E_{\rm mig}^{\rm int}-E_{\rm mig}^{\rm vac})/k_{\rm B}T)$, and since $E_{\rm mig}^{\rm int}<E_{\rm mig}^{\rm vac}$, it decreases at high temperature. At low temperatures, the vacancy diffusion rate is simply too slow for sufficient vacancies to cluster and form voids on experimental timescales.

This simple model is sufficient to qualitatively account for most of the phenomena observed in void lattice formation: the temperature window for formation, bcc superlattices appearing in bcc crystals, and hexagonal superlattices in hexagonal crystals (where diffusion within the basal plane is sufficiently faster than diffusion normal to it to make the superlattices effectively 2D Mazey and Evans (1986)). Our model cannot predict fcc lattices (which have more than one inherent lengthscale). We have also observed the unexpected purely diffusion-driven migration of established voids to superlattice sites. The lineararized Turing analysis predicts analytically the superlattice parameter in excellent agreement with fully nonlinear phase field simulations, even when the simulations are initialized with a pre-existing population of randomly distributed voids. The remarkable robustness of stable superlattice formation, together with the simple and general nature of the model, suggests that Turing instabilities and their associated patterns could be generic in many solid state systems, where widely differing diffusivities of different species are ubiquitous.

*Acknowledgments — * SPF acknowledges useful discussions with Dr P. Edmondson and Prof S. Donnelly, and financial support from the UK EPRSC under grant number EP/R005974/1.

I Supplementary material

In order to investigate the robustness of the pattern formation, we investigate the impact of various parameters on the pattern formation. In each of these analyses, we generate an initial set of randomly distributed voids and then let them evolve over time. The interstitial concentration throughout the domain is initialized at a value of $\bar{s}=0.007$.

(I) Impact of the initial average vacancy concentration on pattern formation— Three simulations were conducted, starting with 90, 120, and 150 voids. The initial void radius was 5.8. This results in an initial average vacancy concentration of 0.137, 0.181, 0.225, respectively. In each case, a stable lattice of voids formed in the material. The production term in each was 0.00125.

(II) Impact of the magnitude of the source term— Four simulations were conducted with different values for the defect production term, each starting with the same 90 voids. The four production terms were 0.0005, 0.00075, 0.0001, and 0.00125.

(III) Impact of Variation in the Initial Void Size— In the previous simulations, each void started with the same size. Now, we compare the impact of randomly varying the initial size of the voids. In each case, the average void size is 5.8 and we start with 120 voids. We run one simulation with no variation, one in which the void size uniformly varies by $\pm$20% of the void radius, and one which varies by $\pm$40% of the void radius.