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.

Introduction

Ordering and self-organization in reaction-diffusion systems are of great significance in determining the ordered patterns in chemistry, physics and biology[1-7]. Pioneered by Alan Turing, the mathematical description of reaction-diffusion systems are well established, and it has been recognized that the Turing instability originates from the reaction kinetics and the ordering is dictated by the breaking of continuous translational symmetry (i.e., the loss of homogeneity)[8,9]. However, the studies on another type of patterning mechanism in reaction-diffusion systems are limited, which originates from thermodynamic instability[10,11] and the breaking of point symmetry (e.g., rotational or mirror symmetry)[12]. In physics, the self-organizations induced by thermodynamic instability and the breaking of point symmetry are widely observed in systems without reactions, e.g., multi-domain patterning in second-order ferroelectric and ferromagnetic phase transitions, which is dominated by long-range electric/magnetic interactions[13-16]. In those cases, long-range interactions play a critical role in pattern formation. In parallel, if local reactions are considered instead of non-local interactions, one should expect a new type of patterning mechanism in reaction-diffusion systems. Here we report a unique self-organization mechanism to understand the formation of void/gas bubble superlattices in crystals, which originates from the interplay of thermodynamic instability, diffusion anisotropy and reaction kinetics.

Literally, the dynamics of reaction-diffusion systems is dictated by the coupling of reaction and diffusion. Most of the previous studies focus on different kinds of reactions, while diffusivities are usually taken as positive (i.e., down-hill diffusion) and isotropic for simplicity[17,18]. In such cases, local reactions dominate the instability and the breaking of translational symmetry, while diffusion plays a secondary role on the ordering process. However, we cannot ignore another possibility that diffusion plays a dominate role over local reactions, when the diffusion flux is against concentration gradient (i.e., up-hill diffusion) and/or the diffusivity is anisotropic. In physics, up-hill diffusion could arise from thermodynamic instability[19,20], while anisotropy suggests the breaking of point symmetry[21]. Given the above two fundamental pieces, instability and ordering could be realized even with a very simple reaction. Note that the reaction term in this case should not follow the constraint in Turing instability[18], since it does not directly contribute to the instability and ordering. However, additional phenomena associated with local reactions could occur, e.g., bifurcation[17,22], which might affect the symmetry selection of ordered patterns. In general, the bifurcation associated symmetry breaking could lead to diversified types of self-organized patterns. In real material systems, reaction-diffusion with the coupling of thermodynamic instability and breaking of point of symmetry can be found in crystals under irradiation. In particular, self-organized void/gas bubble superlattices have been widely observed in a large number of metals and alloys under irradiation[23-27]. In those systems, the reaction and diffusion of point defects are coupled with void formation and 1-dimensional interstitial diffusion, which result in diversified types of void/gas bubble superlattices as reported in the literature[28-30].

In this paper, we investigate the ordering and self-organization in a reaction-diffusion system with thermodynamic instability and kinetic anisotropy, through a combination of analytical study and phase field simulations. We establish a general theoretical framework to predict the symmetry of superlattices associated with anisotropic diffusion. In particular, we demonstrate a unique formation mechanism of superlattice structures dictated by the interplay of diffusion anisotropy and local reactions. For a fixed type of anisotropic diffusivity, reaction rate works as a bifurcation parameter that could lead to different superlattice symmetries. Our discovery suggests a new way to design and control the symmetry of void/gas bubble superlattices in solid crystals under irradiation.

Analytic Study of Reaction-Diffusion Systems with Kinetic Anisotropy

Without the loss of generality, our analytical study starts with a generic description of reaction-diffusion systems, followed by specific anisotropies incorporated into phase field modeling and simulations. Mathematically, the dynamics of reaction-diffusion systems is described by partial differential equations including two kinds of terms, i.e., diffusion terms and reaction terms. In this study, we focus on the instability and ordering caused by different diffusion terms. Here we consider two ways to modify the diffusion terms, by adopting Cahn-Hilliard type diffusion[31] and anisotropic diffusion. The former could introduce thermodynamic instability and up-hill diffusion, while the latter suggests a breaking of point symmetry. In a reaction-diffusion system, we consider two components, the diffusion kinetics of which are dominated by Cahn-Hilliard equation and anisotropic diffusion, respectively. A simple reaction term, i.e., annihilation, is considered to couple the evolutions of the two components. Source terms are also considered to balance the annihilation. The following equations are employed,

X and Y are the concentrations of the two components. The diffusion of X is Cahn-Hilliard type, with a mobility of M. The diffusion of Y is anisotropy, and there could be n-th types of Y with different diffusivity D. P and P are the source terms for X and Y, respectively. F is the total free energy of the system. K is the reaction rate for the annihilation between X and Y.

F is the total free energy of a non-uniform system, which 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.

Here we consider an anisotropic diffusivity along a 1-dimensional (1D) direction, which can be represented in a tensor form.where m is a unit vector describing the direction of the 1D diffusion for Y, and D0 is the 1D diffusivity along m direction. ⊗ is the diadic product operator.

For an analytic study, we consider a linear approximation of the reaction term.

and are the averaged values of X and Y, respectively, which are spatially-independent. As a result, we can express and , where δX and δY are perturbations. The above approximation is valid when the pertubations are relatively small (comparing with the averaged values), which corresponds to the initial stage of modulation. As a result, Eqs 1 and 2 can be represented in Fourier space.

and are the Fourier transforms of X and Y, respectively. Excluding k = 0, and are also the Fourier transforms of δX and δY. G is the coefficient matrix (a square matrix of n + 1 order) of the above partial differential equations. In the index of G, subscript 1 indicates X, and indicate n types of Y. G is a function of the wave vector k, as well as all the above thermodynamic and kinetic parameters (i.e., F, , , M, D, K).

Assuming the first developed wave is k. The critical conditions in determining k should include the following two equations,

Note that k is a wave vector having both magnitude k and direction , which determine the characteristic wavelength and symmetry of a superlattice, respectively. In order to simplify the above equation, we consider n equivalent types of Y (in terms of their diffusion anisotropy), with . As a result, we can separate the length and direction of k in det[G(k)].wherea1, a2 and a3 are parameters depending on the magnitude of  k only, while other part of Eq. 10 depends on the direction only. The solution of k corresponds to a minimum of det(G) if n is odd, while it corresponds to a maximum if n is even. In the following discussion, we take an example of the 1D diffusion along equivalent directions in a Cartesian coordinate system, i.e., , , , , with n = 4. For any given values of a1, a2 and a3, we can determine the critical direction that maximizes det[G(k)]. In the above equation, a1 and a2 are parameters being coupled with such a maximization process, while a3 is decoupled. As a result, a numerical calculation to maximize det(G) (with respect to ) can be performed for a given set of a1 and a2.

For 1D diffusion of Y along equivalent directions, since a1 is negative near the critical point, we use −a1 in our following discussions for convenience. The symmetry selection of superlattice associated with type 1D diffusion is shown in Fig. 1. The horizontal axis is −a1 with logarithmic scale, while the vertical axis is a2. Depending on the choice of −a1 and a2, there are four distinctive regions, in which different k directions are dominant. When both −a1 and a2 are large, is the preferred k direction (yellow region). When both −a1 and a2 are small, (100) is the preferred k direction (green region). Between the yellow and green regions, there is a region preferring the direction (blue region). There is also a slim transition region (red region noted by T) between the blue and green regions, with the preferred k direction changing gradually from to (100). Note that the four regions in Fig. 1 only illustrate the effects of a1 and a2 on the symmetry selection based on Eq. 9. The critical condition described by Eq. 8 should also be taken into account, which leads to another relation between a1 and a2.

Figure 1

Symmetry selection of superlattices associated with type 1D diffusion. The dominant concentration waves in different colored regions are distinctive. The critical condition for the development of a concentration wave is described by solid lines (different types of wave are suggested by different colors). The interplay of the two conditions implies the formation of the superlattices, i.e., the blue line in the blue region suggests BCC superlattice formation, and the yellow line in the yellow region suggests FCC superlattice formation. Numerical simulations are performed at the star points.

Here a1 is a function of a2 for a given direction of k, which suggests a curve including all critical points. In Fig. 1, we plot three solid lines for along (yellow curve), (blue curve) and (100) (green curve), respectively. As a result, when a curve is located in a region with the same color, it suggests the formation of superlattice. In Fig. 1, there are two possible regions of superlattice formation. In the region circled by magenta dashed lines, the preferred k direction is , which results in a face-centered cubic (FCC) reciprocal lattice and a body-centered cubic (BCC) real superlattice. In the region circled by brown dashed lines, the preferred k direction is , which results in a BCC reciprocal lattice and an FCC real superlattice. Note that superlattice could form when a1 is slightly larger than its critical value described by Eq. 14, which suggests that superlattice formation region should attach to the left of the critical curve.

Similar analysis can be applied to the systems with 1D diffusion along directions, i.e., , , , , , , with n = 6. The preferred k direction is always type (for any given a1 and a2 constrained by Eq. 14). As a result, the k vectors suggest a BCC reciprocal lattice and an FCC real superlattice. The mathematical expressions of G as well as detailed predictions of characteristic length and symmetry of superlattices are presented in Supporting Information, which are consistent with previous theoretical studies[32,33].

Phase Field Modeling and Simulations of the Formation of Self-Organized Superlattices in Crystals

We perform phase field modeling and simulations (in a 100 × 100 × 100 simulation cell)[34-36] to validate our analytical predictions. Modeling details are presented in Supporting Information. By choosing different sets of a1 and a2, we obtain an FCC superlattice associated with type 1D diffusion (Fig. 2(a)), BCC and FCC superlattices associated with type 1D diffusion (Fig. 2(b,c), corresponding to the magenta and brown stars in Fig. 1, respectively). The related parameters are listed in Table 1. As expected, the simulation results perfectly agree with our analytical predictions shown in Fig. 1.

Figure 2

Phase field simulation results of superlattices in reaction-diffusion systems with kinetic anisotropy (white: X-rich, black: X-lean). (a) an FCC superlattice associated with type 1D diffusion; (b) a BCC superlattice associated with type 1D diffusion; (c) an FCC superlattice associated with type 1D diffusion.

Table 1

Parameters for phase field simulations.

Example	a 1	a 2	1D diffusion	Superlattice	Lattice constant
Figs 2(a), 3 and 4	−157.3	0.0101	〈110〉	FCC	16.0
Figs 2(b) and 5	−101.8	0.0105	〈111〉	BCC	13.3
Figs 2(c) and 6	−2.37	1.26	〈111〉	FCC	26.7

To further understand the coupling between different components, we plot the concentrations of X in Fig. 3, for type of 1D diffusion. At the initial stage, a chessboard-like modulation of X is developed (Fig. 3(a)), which finally evolves to an FCC superlattice (Fig. 3(b)). When the superlattice forms, the concentrations of six types of Y are plotted in Fig. 4, which correspond to the six equivalent directions of type. It can be found that a 1D diffusion along [110] weakens the modulation along this [110]. The interplay of six types of Y modulations are coupled with concentration wave of X, which finally results in an FCC superlattice.

Figure 3

Simulation results of X concentration for an FCC superlattice associated with type 1D diffusion, corresponding to Fig. 2(a). (a) Initial stage: chessboard-like modulation; (b) final stage: FCC superlattice.

Figure 4

Simulation results of six types of Y concentrations for an FCC superlattice associated with type 1D diffusion, corresponding to Fig. 2(a).

Similarly, the concentration modulations of four different types of Y are shown in Figs 5 and 6, for type of 1D diffusion. According to our previous calculations, either BCC or FCC superlattice can form. The concentration modulations of associated with a BCC superlattice is shown in Fig. 5, while those associated with an FCC superlattice is shown in Fig. 6. Comparing the figures in Fig. 5, we can clearly identify the difference in the concentration modulations of Y, i.e., different symmetry. Note that the symmetry of modulations originates from the anisotropic diffusivities of Y. Without changing the anisotropy, i.e., type of diffusion, we can get completely different symmetries, which is caused by the change of scalar parameters (i.e., a1 and a2). For simplicity, we only consider one independent variable, a2, since a1 is a function of a2 at the critical condition (Eq. 14).

Figure 5

Simulation results of four types of Y concentrations for an BCC superlattice associated with type 1D diffusion, corresponding to Fig. 2(b).

Figure 6

Simulation results of four types of Y concentrations for an FCC superlattice superlattice associated with type 1D diffusion, corresponding to Fig. 2(c).

In fact, a2 reflects a competition between reaction (K) and diffusion (D0), and it is also influenced by the critical wavelength (λ = 2π/k) and the critical concentration when the homogeneous system starts to lose stability. In reaction-diffusion systems, it is well known that the change of reaction/diffusion rate could lead to bifurcation phenomena. Here a2 can be taken as a bifurcation parameter. Note that a2 is a scalar without any symmetry information. However, the change of a2 switches the symmetry of a superlattice. Such a unique bifurcation phenomenon could provide a new insight into the formation of void and gas bubble superlattices in crystals induced by irradiation.

The reaction-diffusion system we discussed above could correspond to a solid crystal under irradiation. For example, a lattice atom can be knocked out by implanted particles (e.g., ions or neutrons), which generates a frenkel pair, i.e., a vacancy and an interstitial atom[37]. In crystals, vacancies and interstitials (also called self-interstitial atoms, SIAs) are point defects of opposite nature (e.g., defect and anti-defect), which can disappear through recombination, i.e., annihilation. The concentrations of vacancies and interstitials can be described by X and Y. In the literature, it has been reported that SIAs and their clusters (e.g., interstitial loops) usually diffuse along a 1D crystallographic direction because of the lattice discreteness of crystals[38,39]. As a result, Eqs 1 and 2 capture the production, reaction and evolutions of defects in crystals under irradiation, in which the rate theory for production and reaction kinetics[40] and the Cahn-Hillard approach for the phase separation description of void formation[19,41] are coupled together. Thermodynamic instability are described through the formulation of F in Eq. 1, which leads the accumulation of vacancies and the formation of voids[28-30,34,42,43]. The types of SIAs depend on the symmetry of host crystals. For example, in a BCC crystal, if the SIAs diffuse in 1D along the close-packed directions 〈111〉, there are four types of SIAs (n = 4), which diffuse along , respectively. In an FCC crystal, if the SIAs diffuse in 1D along 〈110〉, there are six types of SIAs or clusters (n = 6), which diffuse along , respectively. Note that we treat the SIAs diffusing along different directions as distinct types of Y, rather than treating them as a single type with equal diffusivity along several crystallographically equvalent directions. Theoretically, diffusivity is a second-rank tensor, which has to be isotropic in cubic crystal[21]. For example, if an interstitial atom jumps along four 〈111〉 directions with equal probability in one step, it essentially diffuses isotropically after several steps. Furthermore, it is possible that an interstitial atom does not jump strictly along one direction in reality, i.e., it may change its direction. In such a case, reaction terms among different Y should be included. However, we do not consider those terms due to analytical complexity.

As suggested by experimental observations, only FCC superlattice can form in FCC host crystals[27,44,45], while either BCC or FCC superlattice can form in BCC host crystals[46-48], which agree with our theoretical predictions. As suggested by our analyses, the symmetry selection of superlattices is dictated by a1 and a2, which are determined by the interplay of thermodynamic/kinetic properties of the material systems and radiation conditions (details can be found in Supplemental Information). Void/gas bubble superlattices widely observed in experiments also suggest the stability of superlattice under irradiation. In our phase field simulations, we obtain stable superlattices upon further relaxation, which is caused by the dynamic equilibrium between defect generation and recombination. Without irradiation (e.g., the reaction terms), the superlattice is not thermodynamically stable. Coarsening can occur driven by the minimization of surface energy. Such coarsening will be very slow for order superlattices. The study of superlattice stability is beyond our symmetry analyses in this paper, which will be conducted in future work.

Conclusion

Through analytical study and phase field simulations, we establish a general theoretical framework to predict the ordering and self-organization in reaction-diffusion systems with thermodynamic instability and kinetic anisotropy. It is found that the pattern symmetry is determined by the interplay of anisotropic diffusions and local reactions. A new bifurcation phenomenon is demonstrated, which provides a new insight into the formation mechanism of irradiation-induced void/gas bubble superlattices in crystals.

Supplemental Info

SI LaTeX File