Bifurcation and Pattern Symmetry Selection in Reaction-Diffusion Systems with Kinetic Anisotropy

Ordering and self-organization are critical in determining the dynamics of reaction-diffusion systems. Here we show a unique pattern formation mechanism, dictated by the coupling of thermodynamic instability and kinetic anisotropy. Intrinsically different from the physical origin of Turing instability and patterning, the ordered patterns we obtained are caused by the interplay of the instability from uphill diffusion, the symmetry breaking from anisotropic diffusion, and the reactions. To understand the formation of the void/gas bubble superlattices in crystals under irradiation, we establish a general theoretical framework to predict the symmetry selection of superlattice structures associated with anisotropic diffusion. Through analytical study and phase field simulations, we found that the symmetry of a superlattice is determined by the coupling of diffusion anisotropy and the reaction rate, which indicates a new type of bifurcation phenomenon. Our discovery suggests a means for designing target experiments to tailor different microstructural patterns.

The total free energy of a non-uniform system can be described as below, with the gradient term incorporated.
f is the bulk free energy density, and κ is the coefficient of gradient energy, which captures the energetic penalty of inhomogeneity.
The general way to describe a 1D diffusion is through a diffusivity tensor (Eq. 3), and its component form can be written as Considering Eqs. S3 and S5, we can express the Fourier transform of Eqs. S1 and S2 as f is the second-order derivative of f , with respect to X. In the index of G, subscript 1 indicates X, and 2 ∼ (n + 1) indicate n types of Y i . Only the diagonal terms G ii depends on k. The non-zero off-diagonal terms are G 1,i+1 = −KX and G i+1,1 = −KȲ i . The solution of the above partial differential equations relies on the eigenvalues of the coefficient matrix (G), which depends on k. The eigenvalues of G, R(k), can be determined by For a given k, if there is one eigenvalue larger than zero, the system will lose stability to the X wave with wave vector k, i.e., such a wave will develop. At the critical point that the system starts to lose stability, the largest eigenvalue is a maximum at zero, with respect to k.
For the largest eigenvalue, the critical condition is described by The above problem for the extremum of eigenvalue (of coefficient matrix G) can be converted to a problem for the extremum of determinant, so that the conditions in determining the critical wave vector k c should include the following two equations.
which can be expressed as below, with all the non-zero terms in G considered.
Note that the solution k c corresponds to a minimum of det(G) if n is odd, while it corresponds to a maximum if n is even. For simplicity, we assume all types of Y i are equivalent (could be due to crystallographic degeneracy), so thatȲ i =Ȳ . Here we consider the following cases.
Case (i): There is only one type of Y i . Eqs. S13 and S14 are reduced to Since G 12 and G 21 are independent of k, Eq. S16 can be written as It is clear that the solution of the above equation is With Eq. S15 considered, we can obtain k c .
In this case, there are infinite number of the critical wave vectors perpendicular to m 1 (i.e., if m 1 is a plane normal vector, all the critical wave vectors are in this plane), with the above magnitude.
Case (ii): There are two types of Y i . Eqs. S13 and S14 are reduced to Similar to Case (i), we can get m 1 and m 2 are the 1D diffusion directions for two different types of Y i . In this case, there is a unique critical wave vector perpendicular to both m 1 and m 2 , with the above magnitude.
Case (iii): There are n types of Y i (n 3). Eqs. S13 and S14 can be written as Because of the complexity of the above equations, we reduce Eq. S25 to a dimensionless form, which is used in numerical calculations.
a 1 , a 2 and a 3 only depend on the magnitude of k, i.e., they are independent of the direction of k. The effects of a 1 , a 2 , a 3 on the symmetry selection of superlattices are presented in the paper.
When the above equations are applied to understand the formation of irradiation-induced void/gas bubble superlattice, a 1 , a 2 and a 3 are determined by the interplay of thermodynamic/kinetic properties of the material systems and radiation conditions. a 1 depends on the diffusivity of interstitials (D 0 ), the mobility of vacancies (M ), average concentrations of vacancies and interstitials before modulation occurs (X andȲ ), chemical free energy (f ), surface energy (κ), recombination rate between vacancies and interstitials (K) and the critical wave vector of superlattices (k). a 2 depends on K,X, D 0 and k. a 3 depends on K,X,Ȳ , D 0 and k. According to unit analyses, a 1 and a 2 are determined to be dimensionless, which suggests that they are universal parameters to determine the symmetry selection in such reaction-diffusion systems.
In Case (i) and (ii), the critical vector k c can be analytically determined because the direction and magnitude of k c could be easily separated in the equations. However, in Case (iii), when the number of Y i is larger than (or equal to) 3, the direction and magnitude of k c are coupled together, so that analytical solution cannot be obtained. Instead, we have to adapt numerical calculations based on the values of a 1 , a 2 and a 3 . It can be found that the direction of k c is determined by a 1 /a 2 /a 3 , while a 1 /a 2 /a 3 are functions of k c magnitude. The coupling between the direction and magnitude of k c is also the mathematical origin of the bifurcation phenomenon presented in our paper.

II. PHASE FIELD MODELING OF SUPERLATTICE FORMATION
A phase field model to investigate the superlattice formation is developed based on the open source MOOSE finite element framework. In MOOSE, the phase field kinetic equations (e.g., Cahn-Hilliard equation) are implemented in a general form which is separated from the thermodynamic information. In addition, the spatial discretization is finite-element-based, which can accommodate various kinds of geometries and boundary conditions. Implicit time integration is utilized with adaptive time stepping, employing the preconditioned Jacobian-free Newton Krylov method and utilizing capabilities provided by the libMesh and PETSc libraries. The MOOSE code is open access and guidance is available for downloading, compiling and using the code: https://www.mooseframework.org/. The results in the present paper can be fully reproduced with this version of MOOSE: commit f6fb066963c59cb24b2f1aea175ee338dc3e40a0 (HEAD, origin/master) Merge: 20f6c3043f7bb602ed94f219c4a42a708ee1ecc7 Date: Wed Apr 24 2019 The MARMOT code, which is based on MOOSE and is free to access, requires a license. Request can be submitted at: https://moose.inl.gov/marmot/SitePages/Access.aspx.
Without loss of generality, a genetic formulation of free energy is considered, which suggests an intrinsic thermodynamic instability of the system.
We use a 100 × 100 × 100 simulation cell, which can capture the formation of FCC and BCC superlattices with given input parameters. Note that the critical wavelength of the superlattice is determined by the input parameters. A smaller simulation cell is also acceptable, if the critical wave length is smaller. Periodic boundary conditions are applied to field variables, X and Y i . Initially, it is a uniform system with X = Y i = 0. All simulations are performed using the Falcon supercomputing system at Idaho National Laboratory, with 36 cores, 128G Random Access Memory and 60 hours. The dimensionless parameters used in phase field simulations are listed in Table S1. Note that phase field is a deterministic method, and there is no uncertainty in calculations beside the intrinsic rounding error when a double precision float is stored in the computer memory. The critical condition for a perfect superlattice formation is described by the solid lines (e.g., blue line for BCC, yellow line for FCC) in Fig. 1. If a 1 and a 2 gradually change their values deviating from this line, the ordering of the superlattices will become worse and worse, so there is no sharp boundary between lattice and non-lattice regions.
To further confirm the symmetry of superlattices, we plot the superlattice configurations in different viewing directions in Figs. S1-S3, corresponding to Figs. 2(a)-(c) in the paper, respectively.