Symmetry breaking during defect self-organization under irradiation

One of the most intriguing phenomena under radiation is the self-organization of defects, such as the void superlattices, which have been observed in a list of bcc and fcc metals and alloys when the irradiation conditions fall into certain windows defined by temperature and dose rate. A superlattice features a lattice parameter and a crystal structure. Previously, it has been shown that the superlattice parameter is given by the wavelength of vacancy concentration waves that develop when the uniform concentration field becomes unstable. This instability is driven thermodynamically by vacancy concentration supersaturation and affected by the irradiation condition. However, a theory that predicts the superlattice symmetry, i.e., the selection of superlattice structure, has remained missing decades after the first report of superlattices. By analyzing the nonlinear recombination between vacancies and self-interstitial-atoms (SIAs) in the discrete lattice space, this work establishes the physical connection between symmetry breaking and anisotropic SIA diffusion, allowing for predictions of void ordering during defect self-organization. The results suggest that while the instability is driven thermodynamically by vacancy supersaturation, the symmetry development is kinetically rather than thermodynamically driven. The significance of SIA diffusion anisotropy in affecting superlattice formation under irradiation is also indicated. Various superlattice structures can be predicted based on different SIA diffusion modes, and the predictions are in good agreement with atomistic simulations and previous experimental observations.


Introduction
Self-organization and pattern formation have been widely observed in various far-fromequilibrium systems. As an effective tool to create far-from-equilibrium environments, radiation has been commonly known to produce point and extended lattice defects in crystals, usually with disordered distributions. The surprising observations of selforganization phenomena under radiation such as compositional patterning (Enrique and Bellon 2000), nanodroplet (Wei et al. 2008), periodic walls of vacancy loops and stacking fault tetrahedrons (Jager and Trinkaus 1993), and void superlattice (Evans 1971), have caught intensive research interests. In particular, the formation of void or bubble superlattices has been reported in both metallic and ionic systems under different radiation conditions including neutron, ion, and electron irradiation, as summarized in previous reviews (Krishan 1982;Ghoniem et al. 2001). The self-organization of void superlattices is not only scientifically fascinating, but also a promising way to mitigate swelling and He embrittlement, which are critical concerns on the structural materials used in nuclear reactors, by forming high density of nanoscale voids/bubbles instead of large ones. Similar mitigation may also be achieved by tailoring the spatial distribution of voids/bubbles using nano-layered alloys (Chen et al. 2017). Moreover, it suggests a promising way of tailoring material microstructure using radiation for novel properties. However, albeit their strikingly interesting nature and numerous observations, the governing physics that leads to superlattices formation is yet to be fully discerned.
Similar to atomic lattices, void superlattices adopt certain crystal structures. The crystal structure of an atomic lattice is usually defined by long-range and short-range thermodynamic interactions between atoms. Similar interactions exist between nanoparticles forming superlattices (Shevchenko et al. 2006). As to void superlattices, a possible such interaction is the elastic interaction between voids, which can be anisotropic and leads to void ordering (Malen and Bullough 1971). However, this ignores the influence of radiation conditions, against experimental observations, and is unable to explain the formation of superlattices in elastically isotropic W (Sikka and Moteff 1972). Without such thermodynamic interactions, what drives the ordering of voids?
It has been long realized that during defect accumulation under continuous irradiation, the mean defect concentration fields can lose stability to modulated waves at a certain critical point (Krishan 1980;Walgraef et al. 1996). Such a dynamic instability breaks the translational symmetry, thereby giving a characteristic length which is related to the superlattice parameter. In particular, by introducing thermodynamic description of vacancies, an analytic solution is derived for the characteristic length (Gao et al. 2018a). Similar to the so-called chemical freeing of phase separation in reaction-diffusion systems (Carati and Lefever 1997), such a characteristic length is stabilized by the competition between phase separation kinetics and defect dynamics. However, the appearance of a single characteristic length is insufficient to define a lattice. Therefore, further symmetry breaking process is required to explain the selection of superlattice structure. Two main hypotheses have been proposed in the literature; both have been successful in explaining superlattice formation in different material systems without necessarily excluding each other. The first one states that, based on experimental observations, the superlattices are isomorphic with the hosting matrices for unknown reasons. This theory, albeit the unknown physics behind, has remained to be consistent with most experimental results except for recent observations of fcc superlattices in bcc UMo (Gan et al. 2015) and bct superlattices in bcc Mo (Sun et al. 2018). The second is the so-called shadowing effect. Assuming that SIAs perform one-dimensional (1D) diffusion along a family of crystal axes, voids aligned along the SIA diffusion directions can shield each other from being annihilated by incoming SIAs (Forreman 1972;Woo and Frank 1985). This hypothesis is appealing because it is consistent with most experimentally observed superlattice structures. Considering 1D SIA diffusion, superlattice formation has been reproduced by kinetic Monte Carlo (KMC) (Heinisch and Singh 2003) and phase field (Hu and Henager 2009) simulations. The 1D SIA diffusion assumption is also supported by atomistic calculations in bcc metals (Nguyen-Manh et al. 2006). However, it is regarded as against the common observations of planar void ordering at early stage of superlattice formation. Alternatively, the planar ordering can be explained by 2D SIA diffusion (Evans 2006). Such controversies call for a fundamental understanding of the physical connection between anisotropic defect diffusion and symmetry breaking.
Following the rate theory based instability analysis (Gao et al. 2018a), this study addresses possible symmetry breaking induced by anisotropic defect diffusion. By analyzing the nonlinear recombination between vacancies and SIAs in the discrete lattice, it is found that, when instability in the vacancy concentration field occurs, the selection of wave directions depends on SIA diffusion anisotropy. The recombination rate becomes dependent on the directions of perturbation waves. Minimums in the recombination rate are obtained for waves that are along certain crystal directions. As a result, these waves are favored to grow over others , because the recombination rate contribute negatively to the growth rates of perturbations. This suggests a theoretical connection between anisotropic defect diffusion and symmetry development during defect self-organization under irradiation. The corresponding superlattice structures that develop after instability can be predicted by analogizing the wave planes to the close-packed planes in crystals. With various SIA diffusion modes, the theoretically predicted superlattice structures are consistent with those from atom kinetic Monte Carlo (AKMC) simulations and experimental observations.

Methodology and computation details
The rate theory framework for defect evolution under irradiation Although the current work focuses only on the nonlinear recombination of vacancy and SIAs, the rate theory framework for instability analysis (Gao et al. 2018a) is still needed to show how the recombination rate affects the growth rates of perturbation waves, and it is presented in the following. Theoretical predictions of the formation condition (temperature and dose rate) and the characteristic wavelength have been done elsewhere (Gao et al. 2018a).
In the rate theory, the lattice defects produced by irradiation are represented by the concentrations of vacancies (c v ) and SIAs (c i ) (Bullough et al. 1975). Replacing the Fickian diffusion by free energy (F) driven Cahn-Hilliard dynamics (Cahn 1961) leads to the below evolution equations for c v and c i (Gao et al. 2018a): In the above equations subscripts i, v and s represent SIA, vacancy, and sink, respectively. P is the production rate in unit of displacement-per-atom per second, dpa/s. It is multiplied by (1 − c v ) to ensure defect production in the bulk only but not in the voids where c v = 1. M is the atomic mobility and is given by M = D/2K B T following the constant mobility version of the Cahn-Hilliard equation for void superlattice formation (Hu and Henager 2009). Here D is the diffusivity, K B the Boltzmann constant and T temperature. K iv is the recombination constant between vacancy and SIA, usually expressed as K iv = 4πR iv (D i +D v ) , with R iv being the recombination radius between vacancy and SIA, and being the atomic volume. k 2 Xs is the sink strength, with X = i for SIA and X = v for vacancy, respectively.
The terms on the right-hand side of Eq. 1 represent in turn production, transport, mutual recombination between vacancies and SIAs, and their absorption by sinks. Compared to the often used mean-field rate theory framework (Bullough et al. 1975), the above description considers the thermodynamics of vacancies with spatially dependent vacancy concentration. Such a consideration is of particular importance as it allows the formation and migration of voids driven by the free energy. The voids, effectively represented by c v =1 in local regions, can precipitate out as void phase precipitates from a over-saturated matrix phase, in a way similar to phase separation in immiscible regular solid solution. Mathematically, they can be regarded as perturbations in the concentration field, and can either grow or shrink depending on the growth rate. Such an approach has widely been used for void formation under irradiation using the phase field method (Hu and Henager 2009). Focusing on the vacancy concentration, the evolution rate can be re-written as: The three terms in the right hand side of Eq. 2 describe Cahn-Hilliard phase separation dynamics, chemical reaction and a source term, respectively. Consider a small perturbation wavec v (k) with non-zero k, whose evolution is given by: Here f is the second-order derivative of the bulk free energy density f with respect to c v . κ is associated with surface energy. When c v is low, the growth factor R(k) is negative for any non-zero k, implying that the uniform concentration field is stable. As irradiation dose increases, c v increases. Once c v is above a critical value inside the spinodal zone, the uniform concentration loses stability to perturbation waves with a critical wavelength, whose growth rates transitions the first from negative to positive. The instability is actually a spinodal phase separation process driven thermodynamically by vacancy supersaturation and also affected by the reaction term for defect annihilation. This is similar to the case in a nonequilibrium superconducting film (Scalapino and Huberman 1977). The critical wavelength has been derived elsewhere as λ c = 2π k c = 2π( κM v (P+K iv c i +k 2 vs D v ) ) 1/4 ; here k c is the critical wavenumber (Gao et al. 2018a). Each of the concentration waves represents a close-packed plane of the superlattice, with the wavelength corresponding to the interplane distance. As such, the superlattice parameter is associated with the critical wavelength and it shares the dependence on the surface energy via κ, the vacancy mobility M v , the defect production rate P, and the vacancy annihilation rate (K iv c i + k 2 vs D v ). As the analytical solution shows, the superlattice parameter increases with increasing temperature and decreasing dose rate. It should be noted that such instability takes place only in certain irradiation conditions defined by temperature T and dose rate P, as has been shown in the literature (Ghoniem et al. 2001). At very high dose rate P or low temperature T, the system is in the recombination-controlled regime and the critical vacancy concentration will never be reached (Gao et al. 2018a). At very low dose rate P or high temperature T, sink absorption dominates, and voids may lose thermal stability due to vacancy emission (Ghoniem et al. 2001).
The previous analysis (Scalapino and Huberman 1977;Gao et al. 2018a) concerned only the development of the critical wavelength, which is related to the superlattice parameter but has no contribution to defining the superlattice structure. Note that the wave vector k comprises a wavelength (or wave number, k 0 ) and a directionk. Only the wavelength was considered in the previous analysis. As will be shown below, when SIAs diffuse anisotropically, the growth rate R(k) becomes anisotropic with respect to the wave directionk. Such an anisotropy in R(k) comes from the recombination term. The other terms, such as the even-order terms, the source and the sink terms, are apparently independent ofk. As shown in Eq. 3, the recombination rate correlates negatively with the growth rate R. Therefore, with the same wavelength, waves along directions associated with minimal defect recombination will be favored over others. By numerically assessing the dependence of the recombination rate on the directions of perturbation waves, the most favorable wave directions can be obtained. These wave directions define the normals of close-packed planes and thereby the superlattice structure.

Defect recombination with perturbations
For their opposite nature, vacancies and SIAs recombine immediately when they are within a certain distance -the recombination radius, R rec -from each other. The recombination rate with uniform concentrations c 0 v and c 0 i is given by: In Eq. 4, K i iv and K v iv represent the fractions of recombination induced by SIA diffusion and vacancy diffusion, respectively. N new is the number of lattice sites entering the recombination radius of a point defect after each atomic hop; ω i and ω v are the hopping rates of SIA and vacancy, respectively. In typical materials at the conditions for superlattice formation, ω i >> ω v , so that K iv ∼ = K i iv . In the following, we will focus on SIA diffusion only. In deriving Eq. 4, it is assumed that the defect concentrations are uniform before instability, but the defect fields are dynamically evolving due to diffusion and defect production. With such an assumption, Eq. 4 is applicable for both 1D and 3D SIA diffusion. We note that in the literature a different expression has been used for the recombination rate when SIAs perform 1D diffusion (Amino et al. 2011), which followed the analysis of annihilation of 1D diffusing SIA clusters at static traps such as voids (Barashev et al. 2001). In bcc metals such as Mo and W, SIAs diffuse much faster than SIAs. Accordingly, the SIA concentration is usually much lower than that of vacancy when defect sinks exist. In particular, c v D v ≈ c i D i when the system is approaching the steady state. Because the evolving rate of a defect field due to diffusion depends on both the concentration and the corresponding diffusivity, a similar evolving rate for vacancies is expected to that for SIAs at the steady state. This, plus the change caused by defect production, suggest that the vacancy field is more likely dynamically evolving than being static under irradiation. Such an assumption is backed up by our AKMC simulations in situations when The recombination rate is affected by the appearance of perturbation. Consider a small perturbationc v = δ v exp(ik · r) on top of c 0 v . Due to the nature of mutual recombination, this will induce a simultaneous anti-phase perturbation in c 0 i , denoted asc i = −δ i exp(ik · r). With this, the averaged recombination rate induced by SIA diffusion is given by: Here l i is the SIA mean free life, defined as the average number of hops before an SIA is recombined. It can be estimated as in the recombination dominant situation. q 0 is the averaged probability that an SIA is not recombined after an atomic hop. λ m is the current SIA position after m hops in reference to its original position. An average over λ m is needed for the stochastic nature of atomic hops. For convenience, the symbol σ kλ is used to represent the summation term, i.e., The derivation of Eqs. 4 to 7 is given in the Appendix B of Supplementary Material. The dependence of recombination rate on wave directionk can be elucidated by numerically computing σ kλ , which represents the reduction in recombination rate caused by perturbations. The detailed calculation procedure is given at the end of Appendix B in the Supplementary Material. Unless otherwise stated, a bcc matrix with R rec being the 2 nd nearest neighbor(NN) distance (so that N new = 7) is used. A vacancy concentration of 0.001 is taken to estimate l i . For reference, the vacancy concentrations is at an order of 0.01 (converted from the void size and density) in metals where void superlattices have formed (Ghoniem et al. 2001). We note that an accurate estimate of l i is not actually needed for the purpose of predicting symmetry development. The same anisotropy in σ kλ holds as long as l i is larger than the recombination radius. In cases that SIAs perform 1D diffusion along symmetrical directions, e.g., the 111 family in bcc, or 2D diffusion within symmetrical planes, SIAs are divided into n i groups equally, with n i being the fold of symmetry of SIA diffusion. The calculation of σ kλ runs over all n i groups and is then normalized by n i . The numerical results of σ kλ are plotted as a function of wave directionsk to elucidate the anisotropy.

Atomic kinetic Monte Carlo simulations
To demonstrate the symmetry breaking predicted by Eqs. 5 and 7, rigid lattice AKMC simulations are carried out following the method in REF (Gao et al. 2018a). Here, vacancy and SIA are denoted as types of elements occupying and diffusing on a prescribed, static lattice. Three types of atoms are used to represent matrix atom, vacancy, and SIA, respectively. The interaction between atoms are described by pairwise bonds within the 2NN cutoff. The bond energies between matrix atoms and vacancies are fitted using the vacancy formation energy and the surface energy, thereby defining a free energy model analogous to that for a binary regular solution. This is consistent with the free energy model used in the rate theory description. The interactions between SIAs and between SIAs and matrix atoms are ignored for their extremely low concentration in the simulations. The rigid-lattice based AKMC method is used for two reasons. First, the current method meets the requirement on efficiency for simulating large numbers of point defects and long physical time. The more accurate on-the-fly and off-lattice KMC methods, such as the SEAKMC method (Xu et al. 2013), are more powerful in capturing individual defect evolution, but less efficient for long-term evolution of a large amount of defects. For the same reason, the interatomic interaction and the migration barriers are calculated using a pairwise bond model, instead of using an interatomic potential and the advanced barrier searching methods such as the dimer method (Gao et al. 2000). Second, the current AKMC method is more consistent with the rate theory description, thereby suitable for demonstrating the predictions by the rate theory.
For KMC simulations, the system is evolved following the residence-time algorithm (Soisson and Fu 2007). At each KMC step, a list of diffusion events, i.e., atomic jumps of vacancies and SIAs, are created. The probability of each event i is given by Here ν 0 is the attempt rate. One event is randomly selected to advance the system, and the system time is advanced by 1/ i τ i . The diffusion of vacancies and SIAs are simulated by switching with nearest neighbor (NN) atoms. For vacancy, the switching is done with a randomly selected 1 st NN, representing isotropic diffusion. Both isotropic and anisotropic SIA diffusion are simulated. For isotropic diffusion, SIAs are switched with a randomly selected 1 st NN. For anisotropic SIA diffusion, SIAs are switched with a randomly selected 1 st NN along a given crystal orientation for 1D diffusion, or within a given crystal plane for 2D diffusion, respectively. SIAs are divided into n i groups if they diffuse along n i symmetrical directions or in n i symmetrical planes, with each group diffusing along one direction or within one plane. The diffusion barrier is calculated by E a = E 0 + (E f − E i )/2, and updated once the local environment is changed (Soisson and Fu 2007). Here E 0 is the diffusing barrier at the dilute concentration regime, and E f − E i is the energetic difference between the states after and before an atomic hopping, considering the dependence on local environment. Note that this gives a constant activation barrier E 0 for SIAs because they are ignored in the thermodynamic consideration.
In addition to diffusion events, the below athermal events are also simulated: • Defect production: To simulate defect production, Frenkel pairs are inserted at a rate of one pair per time interval t fp , corresponding to a defect production rate P = (t fp * N a ) −1 , with N a being the total number of atoms used in the simulations. Different production rates can be simulated by adjusting t fp . For each insertion, two lattice atoms are randomly selected, one assigned as vacancy and the other as SIA. Such a defect production method is consistent with the rate theory description, and it corresponds to the situation of electron irradiation. • Recombination: After the production or jump of a point defect, possible recombination will be carried out. Defects of the opposite nature within the recombination radius will both be assigned as matrix atoms, to represent recombination. To be consistent with the theoretical analysis, the recombination radius is set to be the 2 nd NN distance.
• Sink absorption: A mean field sink, represented by a mean free jump (Gao et al. 2018a), is implemented for both vacancy and SIA, by using a static sink strength that can be converted from the density of pre-existing dislocations (Wiedersich 1972). Defects that have migrated for times more than their mean free jumps will be eliminated by converting them back to matrix atoms. Alternative to the mean field sink, planar sinks can also be simulated by assigning a region in the simulation cells as sink. Defects that migrate into this region will be converted to matrix atoms, representing sink absorption. The mean field sink is used in the simulations unless stated otherwise. The planar sink is used to study the effect of heterogeneously distributed sinks.
The AKMC simulations are done with either 2D square and 3D cubic cells with periodic boundary conditions (PBCs). The simulation cell sizes are described along with the results. Both the simulation cell size and shape have been varied to confirm that there is no artificial effect caused by the PBC. For each simulation, we start from defect free simulation cells. Frenkel pairs are then introduced at given dose rates for defect accumulation. When there is no defects in the system, the KMC time is advanced by the time interval t fp with a Frenkel pair introduced. We rely on the stochastic nature of KMC simulations for perturbations, meaning that no perturbation is manually introduced in addition to the above KMC events.
The material parameters for bcc Mo (Gao et al. 2018a) are used to derive the bond energies and diffusion barriers. Specifically, the vacancy formation energy is 2.9 eV; the surface energy is 2.95 J/m 2 ; the SIA (vacancy) migration barrier is 0.083 eV (1.45 eV); the attempt rate is 1.0×10 12 /s for both vacancy and SIA; the mean free jump before being absorbed by a sink is 1000 for both vacancy and SIA.
Because the AKMC simulations are used to confirm the theoretical predictions of superlattice structures, we purposely choose irradiation conditions in which superlattice are expected to form. To save computation time, a high dose rate of 9.8 dpa/s and a high temperature of 973 K are used. To confirm that the superlattice structures are not affected by the dose rate and temperature used, simulations with lower dose rates have also been performed. Note that a lower temperature may need to be used when a much lower dose rate is used, as suggested by the P-T diagram for superlattice formation (Ghoniem et al. 2001;Gao et al. 2018a). For clarity, in presenting AKMC results only vacancies are shown as dark dots to elucidate potential ordering.

Symmetry breaking induced by anisotropic SIA diffusion
Several important points are readily clear from Eq. 5. First, the occurrence of perturbations reduces the recombination rate becausec v andc i always occur simultaneously with antiphase due to the nature of mutual recombination, e.g., increasing in one will simultaneously induce a reduction in the other. The reduction is negligible whenc v << c 0 v and c i << c 0 i , with negligible changes in the critical wavelength. Second, the recombination rate depends on the wave vector k (= k 0k ) and the SIA diffusion property via λ m . When SIA diffuses isotropically, σ kλ is symmetrical aboutk and depends only on k 0 . While SIA diffuses anisotropically, the dependence onk can lead to symmetry breaking.
We start with the dependence of σ kλ on k 0 andk before applying the theory to realistic cases. For this purpose, a case study is carried out considering the [1 1 1] 1D SIA direction (not the 111 family) in a bcc crystal. Three wavelengths, 0.5l i (l i k 0 = 4π), 2l i and 10l i , are considered in calculating σ kλ . As shown in Fig. 1, σ kλ is anisotropic aboutk and increases monotonically ask moves away from the (1 1 1) plane towards the plane normal, [1 1 1]the 1D SIA diffusion direction. Following Eq. 5, k · λ m = 0 is satisfied for any k in (1 1 1) plane because λ m parallels [1 1 1], giving a minimum of -1 for σ kλ . The same anisotropy in σ kλ with respect tok holds for all wavelengths, but the degree of anisotropy decreases  Fig. 1a-c. The dependence of recombination on k has profound effects on vacancy ordering. For its negative sign in Eq. 3, recombination tends to diminish perturbation waves. With anisotropic SIA diffusion, recombination becomes anisotropic with respect tok. With the same wavelength, waves along directions that are associated with minimum recombination rates have the highest growth rates, and are thereby favored to grow over others. The development of waves along favorable directions can lead to symmetry breaking during defect self-organization. The anisotropy in recombination fades as the wavelength increases, suggesting weak ordering for long waves, e.g., when the distance between voids is large.

AKMC confirmation of void ordering
Four case studies are designed and carried out to further demonstrate the theoretical predictions of symmetry breaking and void ordering. To focus on symmetry, we further average the results over k 0 by varying k 0 from 2π/l i to π/2, to take out the dependence on wavelength. A wavenumber larger than π/2 may result in unphysically small wavelength. Accordingly, AKMC simulations are carried out to confirm the theoretical predictions. Note that superlattices form only in appropriate conditions defined by dose rate and temperature. To demonstrate superlattice structure selection, the AKMC simulations are done in conditions where superlattices are expected to form.
Case I considers 1D SIA diffusion in a 2D square lattice (N new = 3), Fig. 2a. The recombination is reduced the most along the k 2 ([0 1]), and the least along k 1 ([1 0]), corresponding to minimums and maximums in σ kλ , respectively, as shown in Fig. 2b. As a result, waves along k 2 will be favored to grow over others. When a critical vacancy concentration is reached, a wavelength will be selected at the instability point (Gao et al. 2018a). As such, a [0 1] wave will develop, and voids will precipitate at the wave peaks for the locally high vacancy concentration, with a [1 0] alignment. Eventually, stripes normal to [0 1] form, showing the favorable growth of [0 1] waves caused by [1 0] 1D SIA diffusion.
Case II considers 1D SIA diffusion along [1 1 1] in 3D. As shown in Fig. 3a, with [1 1 1] 1D SIA diffusion, waves with ak in the (1 1 1) plane are equally favored. When projected along [1 1 1], a random distribution of voids (actually a column of voids aligned along [1 1 1]) with a characteristic length should be expected after instability takes place. This is indeed observed in AKMC shown in Fig. 3d. The appearance of a ring in the fast-Fouriertransformation (FFT) image indicates the existence of a characteristic length but not a symmetry. In both case I and II, voids align along the 1D SIA diffusion directions when instability occurs, consistent with the shadowing effect hypothesis (Forreman 1972;Woo and Frank 1985).
In case III, 1D SIA diffusion along both [1 1 1] & [-1 1 1] in bcc is considered. In this case, the minimum recombination is reached atk = [0 -1 1]. This means a [0 -1 1] wave with the selected wavelength will develop at the instability point, as shown by the AKMC simulation in Fig. 3e. The FFT inset also indicates the appearance of a plane wave. It is cri- tical to note that in this case, void alignment along [1 1 1] & [-1 1 1] becomes unnecessary, in contrast to the shadowing effect. Indeed, no such alignment can be seen in Fig. 3e.
In the last case (IV) we consider realistic 1D SIA diffusion along all 111 directions, as in bcc metals suggested by atomistic calculations (Nguyen-Manh et al. 2006). In this case, the minimum recombination rate is reached along all six 110 directions, as shown in Fig. 3c (three at the edge of sphere are not visible due to projection). This is consistent with previous analysis showing that aligned voids along the closedpacked directions receive less recombination and higher vacancy flux (Semenov and Woo 2006). Therefore, six 110 waves with the same wavelength are expected, and their superposition gives a bcc superlattice, as shown in Fig. 3f and the inset FFT. While cases I-III are designed to demonstrate the theory, case IV is realistic and it explains the bcc superlattice formation widely observed in experiments (Krishan 1982;Ghoniem et al. 2001).
The above analysis is not limited to 1D SIA diffusion. Actually, any SIA diffusion mode can be considered by sampling λ m in a statistical way. By varying the SIA diffusion mode, many superlattice structures can be predicted, as summarized in Table 1, including 2D square and hexagonal, 3D simple cubic (sc), bcc and fcc, consistent with experimental observations (Ghoniem et al. 2001;Krishan 1982;Gan et al. 2010;Johnson and Mazey 1980). This confirms that the symmetry of void superlattices is not governed by the matrix structure, but by the diffusion property of SIAs (Woo and Frank 1985), as shown by previous phase field (Hu and Henager 2009;Hu et al. 2016;Gao et al. 2018b) and KMC simulations (Heinisch and Singh 2003;Evans 2006;Gao et al. 2018a). This is also consistent with Walgraef et al. (1996), which showed that a small anisotropy in SIA diffusion is needed for vacancy loop ordering. The good agreement between the theoretical predictions and the simulations and the experimental observations is very encouraging. While the validity of Eq. 4 depends on a critical assumption that the vacancy field is dynamically evolving, instead of being static as in Ref. (Amino et al. 2011), this assumption is not used in the AKMC simulations. The consistency between the predicted results and those from simulations shown in Table 1 suggest that the assumption is appropriate for the conditions considered here. fcc (Gao et al. 2018a;Hu et al. 2016) fcc (Gao et al. 2018a) sq (Gao et al. 2018a) hex (Gao et al. 2018a) Experiment bcc (Evans 1971;Sikka and Moteff 1972) fcc (Gan et al. 2010) * fcc (Johnson and Mazey 1980) * * * The SIA diffusion model in bcc UMo is yet to be confirmed ** In fcc SIA clusters rather than individual SIAs perform 1D diffusion Several important points need to be emphasized on the comparison between theoretical predictions and experimental observations and previous simulations. It is predicted that superlattice structure is dictated by SIA diffusion mode. With the same host matrix structure, e.g., bcc, both bcc and fcc superlattices can form, by 111 and 110 SIA diffusion, respectively. This is consistent with the experimental observations of bcc superlattices in bcc Mo (Evans 1971) and W (Sikka and Moteff 1972), in which SIAs are predicted to diffuse along 111 (Nguyen-Manh et al. 2006), and of fcc superlattice formation in bcc UMo (Gan et al. 2015), in which SIAs are suggested to diffuse along 110 (Hu et al. 2016). In addition, sc superlattices have been observed in the anion sublattice of CaF 2 and SrF 2 , and the cause was attributed to 100 anion SIA diffusion (Johnson and Chadderton 1983). Indeed, it is predicted here that, with 100 SIA diffusion, sc superlattices should be expected. It is also interesting to notice that the same superlattice structure can result from different SIA diffusion modes. For instance, a bcc superlattice can form with either 111 1D SIA diffusion or {1 1 0} 2D SIA diffusion, as has been shown in previous KMC simulations (Evans 2006) and confirmed by our AKMC simulations. In the latter case SIA diffusion along each {1 1 0} plane will stabilize a favorable {1 1 0} wave in c v . It should be noted that previous atomistic calculations suggest 111 SIA diffusion in most bcc metals including Mo and W (Nguyen-Manh et al. 2006).

Characteristic stages of superlattices formation
The above analysis clearly suggests that superlattices form by superposition of symmetrical concentration waves that develop at the instability point. This means that the widely accepted void alignment along the 1D SIA diffusion direction is an unnecessary condition for superlattice formation, although it may still occur in special cases where the preferred wave directions parallel the 1D SIA diffusion directions, e.g., in cases I and II. Instead, the current theory suggests that superlattice formation may possibly take a 3-stage process involving a planar void ordering. The first stage is before instability, in which unstable voids form randomly. They can be regarded as unstable perturbations that decay exponentially. The second stage starts from the instability point, at which selected waves shown as stable voids start to develop. The instability is driven thermodynamically by vacancy concentration supersaturation. At finite temperatures, these waves may not develop at the same time due to stochastic effects. One or more waves will develop first, giving a planar ordering of voids. In the third stage, all preferred waves will develop and form 3D superlattices. Once the superlattices form, 1D alignment of voids may appear as a consequence rather than the reason of superlattice formation. Therefore, its absence at the early stage should not be interpreted as an evidence against 1D SIA diffusion. At low temperatures, the second stage may not be as clear as at high temperatures due to less stochasticity.
The 3-stage superlattice formation process can be clearly demonstrated using AKMC simulations. Here, 1D SIA diffusion along all 111 in bcc Mo is considered. As seen in Fig. 4a, at 0.8 dpa small voids form with a random distribution. At about 1.9 dpa, voids start to exhibit a (1 1 0) planar ordering. Eventually, all {1 1 0} waves develop, giving a bcc superlattice, as shown in Fig. 4c. Planar ordering of voids and gas bubbles have been commonly observed experimentally (Evans 2006;Harrison et al. 2017). It has been interpreted as an evidence against the shadowing effect. However, as shown by the theoretical analysis and AKMC simulations, the planar void distribution may possibly be a characteristic feature of superlattice formation, reflecting the wave nature of superlattice nucleation, regardless of 1D and 2D SIA diffusion. The 3-stage formation process can be seen in ion irradiated Nb (Loomis et al. 1977) but was not connected to the nature of superlattice nucleation.

Discussion
Several assumptions used in the rate theory framework and the AKMC simulations warrant some discussion. The first is the mean-field production of Frenkel pairs, which is close to the condition of electron irradiation, but different from the cases of neutron and ion irradiation ). In the latter two cases a substantial fraction of defects is produced as clusters, and the fraction can be different for vacancy and SIA clusters. This has been treated in the production bias model (Singh and Foreman 1992) for better modeling of swelling induced by neutron and ion irradiation. The direct production of clusters can accelerate vacancy accumulation by reducing the recombination, as suggested by Eq. 5. Moreover, small SIA clusters can rapidly migrate to sinks without recombining with vacancies, again accelerating vacancy accumulation. As a result, at the same temperature and dose rate, vacancies accumulate with a much lower rate under electron irradiation than under neutron and ion irradiation. The slower vacancy accumulation rate effectively delays the occurrence of instability with respect to dose. This, in addition to the difficulty to reach high dose levels, explains why void ordering is rarely observed (Fisher and Williams 1977) in metals or alloys under electron irradiation. In contrast, superlattices have been widely observed under ion and neutron irradiation (Krishan 1982;Ghoniem et al. 2001). Void superlattice formation has been observed in the anion sublattice in fluoride CaF 2 under electron irradiation (Johnson and Chadderton 1983;Ding et al. 2005).
Another major assumption is the omission of SIAs from the thermodynamic description. This effectively excludes the formation of SIA clusters. It has been shown that under irradiation, SIA clusters, either mobile (e.g., loops) or immobile ones (e.g., clusters, loops, and tangled dislocations), can form. The immobile ones can act as sinks for both vacancies and SIAs. As such, the overall sink density will increase. This will affect the wavelength selection at the instability point, but not the symmetry development. The mobile SIA clusters usually take the configuration of loops and migrate along the Burger vector directions. Therefore, they either rapidly reach the sink, accelerating vacancy accumulation, or are annihilated by recombining with vacancies, thereby contributing to the anisotropic recombination due to their 1D migration. Actually, 1D SIA loop migration has been regarded as the reason for void ordering (Dubinko et al. 1989). In metals such as bcc Fe, SIAs perform 3D diffusion and SIA clusters and loops perform 1D diffusion. The void ordering could be caused by SIA clusters and loops, rather that individual SIAs. In this regard, in some material systems SIA clusters could play a significant or even a dominant role in superlattice symmetry selection. Currently, No SIA clustering is considered in the AKMC simulations. Their role in superlattice formation will be investigated in the future by extending the current theoretical analysis and AKMC simulations. This is needed to generalize the current analysis for materials under ion and neutron irradiation and materials in which only SIA clusters and loops diffuse anisotropically.
While demonstrating the theoretical prediction an unrealistically high dose rate has been used to save computation time. To minimize possible artificial effect from dose rate, the formation of bcc superlattice in a bcc matrix with 111 1D SIA diffusion (the case shown in Fig. 4) is repeated with several different dose rates. Four more simulations are performed at 773 K with the dose rate being 10 −3 , 10 −2 , 10 −1 and 1.0 dpa/s, and another simulation at 973 K with the dose rate being 10 dpa/s. The simulation cell sizes are 40×40×40 a 3 0 . Note that for lower dose rate a lower temperature is needed according to the formation window (Ghoniem et al. 2001). For all simulations, bcc superlattices are obtained. In general, at the same temperature, lowering the dose rate leads to an increase in the superlattice constant, consistent with the results in the literature (Gao et al. 2018a). To further test the effect of dose rate, the simulations shown in Fig. 3a and b are repeated using a simulation cell size of 40×40×40 a 3 0 , a dose rate of 10 −3 dpa/s and a temperature of 773 K. Again, the same types of ordering are observed in the simulations. The results from the new simulations are consistent with the theory that the ordering of voids is governed by the SIA diffusion property, as long as the instability takes place.
In the theoretical analysis we focus on the recombination in bulk. A mean-field approach is used for sink absorption. While in reality, geometrical sinks are randomly distributed. The heterogeneity of sinks may also affect void ordering. As an attempt to explore this effect, an AKMC simulation with a planar sink is carried out. The simulation cell size is 40×40×80 a 3 0 , with one atomic plane located at the middle of the Z direction assigned as a planar sink. Point defects that migrate to the sink sites are annihilated by converting them back to matrix atoms. The mean-field sink is not used in this simulation. The dose rate used is 10 −3 dpa/s and the temperature is 773 K. As shown in Fig. 5, the introduction of a planar sink induces a denuded zone in the middle along the Z direction. The width of the denuded zone decreases with increasing dose. Still, a bcc superlattice is observed at 2.0 dpa, as indicated by the projections along the [100] and 111 directions in Fig. 5c and d, receptively. This is also consistent with experiments that superlattice can form in regions very close to grain boundaries such as twins (Wang et al. 2016). We note the geometry of sinks considered here is much more simplified than that in real materials. A dedicated study may be needed to explore the role of sinks including their strength and spatial distribution.
The recombination analysis here is based on the assumption that the SIA diffusivity is much higher than that of vacancy. This is true in most bcc and fcc metals but not necessarily in all materials. When vacancy diffusivity is comparable or even higher than that of SIA, due to the isotropical diffusion of vacancies, the anisotropic effect induced by anisotropic SIA diffusion will become less significant. Different ordering phenomena may occur in such a situation (Gao et al. 2019). It should be noted that a much higher SIA diffusivity than that of vacancy has been postulated as a criteria for superlattice formation (Hu and Henager 2009).

Conclusion
To conclude, this work presents a discrete lattice theory concerning the patterning of voids caused by anisotropic SIA diffusion. The discrete lattice analysis is based on the same rate theory framework that was used to predict the critical wavelength, which corresponds to the void superlattice constant. It is now possible to use the unified rate theory framework to predict both the superlattice constant and structure. Therefore, the theory is expected to be impactful to guide experiments for tailored microstructure under radiation. With validations by atomistic simulations and experiments, several important implications can be drawn from the discovery: i) The superlattice is a superposition of stable perturbation waves at the instability point, whose symmetry is dictated by anisotropic SIA diffusion. This is consistent with the shadowing effect (Forreman 1972;Woo and Frank 1985) and the conclusion of Walgraef et al. that a small degree of SIA diffusion anisotropy is needed for vacancy loop ordering (Walgraef et al. 1996). The symmetry of these waves defines the superlattice structure, which depends on the kinetics of SIA diffusion but not necessarily the host matrices; ii) A 3-stage formation process, from random voids to planar ordering and then 3D lattice, is implied. The widely accepted theory on void alignment along 1D SIA diffusion direction is suggested to be a consequence, rather than the reason of superlattice formation; iii) The effect of anisotropic SIA diffusion is scaled by the SIA mean free life. At conditions for the instability to occur, the SIA mean free life usually varies from a few to hundreds of nanometers, rendering superlattice formation a nanoscale phenomenon; iv) The symmetry breaking is induced by anisotropic defect diffusion via a reaction term. Similar effect may exist in other systems coupling phase separation dynamics and chemical reaction when anisotropic mass transport is involved.
Additional file 1: The file Supplementary Material contains the derivation of Eqs. 4 & 5. The file stripe.mov is an animation of the formation process of Fig. 2(c). The file 3-stage.mov is an animation of the simulation in Fig. 4.