A. Birth and death process

Within the cluster dynamics (CD) approximation, defects are treated as lattice gas. These defects produced by athermal processes will diffuse in crystalline solids and agglomerate in clusters over a time scale $\tau \propto {{\rm{e}}^{{{{E_{\rm{a}}}} \over {kT}}}} \approx {10^{ - 6}}{\rm{s}}$ (where E _a the activation energy varies from few hundreds of MeV to few eV) over a time scale five orders of magnitude larger the production rate in displacement cascades (few picoseconds), and a space scale defined by the mean distance between different sinks and roughly equal to few hundreds of nanometers, i.e., two orders of magnitude larger than the displacement cascade characteristic size. When defects start to agglomerate into clusters, the size of each cluster grows and shrinks by a process of accretion and emission, controlled by the local concentration of point defects in the neighborhood of the cluster. ^{Reference Karlin and Taylor16} Figure 2 displays the coarsening kinetics of lead particles under electron irradiation.

FIG. 2.

Coagulation kinetics of Pb particles produced by 125 keV electron irradiation (at 673 K with a flux of 4.810¹⁸ cm⁻²/s) in lead glass. Reprinted with permission from Ref. Reference Biron50, Orsay university.

An appropriate framework for describing such events is the theory of stochastic processes ^{Reference Nicolis and Prigogine10,Reference Van Kampen11} for which the cluster size appears to be a random variable. The decoupling of time and space scales allows applying the master equation (ME) formalism to discuss the microstructure formation in crystalline solids induced by vacancies and interstitials. This description is extensively used to discuss a large class of phenomena like chemical reactions. ^{Reference Nicolis and Prigogine10,Reference Gillespie17} The ME refers to a Markov process in which state transitions can occur only between neighboring states (i → i + 1 or i → i − 1). Within the CD approximation, ^{Reference Ghoniem and Walgraef9} one assumes in many cases that only monomers (vacancy, interstitial, and gas fission resulting from transmutation) migrate. The probability P(n, t) to create a cluster of size n obeys to the following differential equation for (n ≥ 2):

In this expression, the cluster flux associated with the motion of clusters belonging to a class of size n to the class n + 1 is usually written as J _n→n+1 = λ _n P(n, t) − μ _n+1 P(n + 1, t). In this last expression, λ _n is the probability rate (per time unit) for one monomer to be scotched to a cluster of size n and μ _n is the probability rate for a monomer to leave a cluster of size n. As atomic clustering in many cases ^{Reference Binder and Stauffer18} is driven by the concentration of monomers (vacancies, interstitials, and gas), the complete distribution of clusters (di-mere, tri-mere…) must be computed. This distribution is given by a coupled hierarchy of these discrete equations. Equation (1) casts as

This last equation defines a birth and death (BD) process. ^{Reference Nicolis and Prigogine10,Reference Van Kampen11} The probability rates λ _i and μ _i are highly nonlinear function of the temperature and the flux as well as the cluster size.

In a pure birth process (λ _i = λ; μ _i = 0), P(i, t) is given by a Poisson distribution and the variance of the process is equal to its mean value. It results that the cluster distribution modeled by a BD process is determined not only by the average size of clusters 〈iP(i, t)〉 but also by all its higher moments (fluctuations). Fluctuations (second moments, third moments…) are obviously related to the system size and the correlation length over which two parts of the system “feel” each others. The main limitation of modeling cluster formations under irradiation by a BD process is to accurately determine different transition rates. ^{Reference Waite19,Reference Abromeit and Wollenberg20} These reaction rates do not depend on the different moments of the defect concentration but are the function of the defect concentrations.

B. Rate theory

It has been clearly shown ^{Reference Nicolis and Prigogine10} that if λ _n and μ _n in Eq. (2) follow a Poisson distribution, the cluster size can be computed assuming that λ _i and μ _i depend only the averaged value of defect concentrations. Such an assumption called the RT approximation ^{Reference Nichols21} is extensively applied to describe chemical reactions. The RT can be understood as a mean filed theory as these transition rates calculated neglecting all fluctuations of defect concentrations. Thus, the time evolution of the population of point defects [vacancy (v), interstitial (i), and sinks like free surface, grain boundaries, or dislocations loops (s)] is dictated by the following “chemical” reactions:

In this approach, the underlying lattice is neglected. Assuming a law of mass action can be applied, the number of λ _i and μ _i terms can be reduced, leading to a new set of control parameters K. This set of chemical reactions can be mathematically handled by a set of deterministic equations:

where C _i, C _v, and C _s are interstitial, vacancy, and sink volume concentrations.

Introducing the production rate K _i(v) of free migrating interstitials and vacancies escaping from displacement cascades, these terms can be written as ν_i(v)σ ^d ϕ, where ϕ is the flux of impinging particles, σ ^d is a displacement cross section ^{Reference Sigmund and Gras-Marti22} and is an estimation of the displacement cascade formation averaged over all crystalline directions, ^{Reference Simeone and Luneville23} and ν_i(v) is an efficiency factor varying from few percents to few tens of percent depending on the temperature and the displacement cascade density. ^{Reference Kiritani24} Within the effective medium theory framework, ^{Reference Brailsford and Bullough25} many studies ^{Reference Dunlop, Rullier-Albenque, Jaouen, Templier and Davenas26} were devoted to calculate vacancy–interstitial recombination ${K_{{\rm{iv}}}} = 4\pi r_{{\rm{iv}}}^2{{\left( {{D_{\rm{i}}} + {D_{\rm{v}}}} \right)} \over \Omega }$ , where r _iv is the radius capture in which a vacancy and an interstitial recombines and D _i(v) is the thermal diffusion of interstitials (vacancies) in the medium. The term Ω is the atomic volume. The term ${K_{{\rm{i}}\left( {\rm{v}} \right){\rm{s}}}} = k_{{\rm{i}}\left( {\rm{v}} \right){\rm{s}}}^2{D_{{\rm{i}}\left( {\rm{v}} \right)}}\left[ {{C_{{\rm{i}}\left( {\rm{v}} \right)}} - C_{{\rm{i}}\left( {\rm{v}} \right)}^{{\rm{eq}}}\left( T \right)} \right]$ , where $k_{{\rm{i}}\left( {\rm{v}} \right){\rm{s}}}^{ - 1}$ is a mean free path associated with the diffusion of a point defect to a sink and has the dimension of a length and $C_{{\rm{i}}\left( {\rm{v}} \right)}^{{\rm{eq}}}\left( T \right)$ is the equilibrium concentration of point defects in the medium. The temperature appears in the recombination terms assumed to be thermally activated. ^{Reference Nichols21} The term K _is (K _is respectively) controls the vacancy (interstitial) concentration in the vicinity of a given sink (dislocation line, cavity, bubble…). The effect of the sink strength mediated by the elastic/electric fields governs these terms.

This approach can be extended to study the growth rate of precipitates, dislocation, and voids by adding specific evolution laws. ^{Reference Abromeit and Wollenberg20,Reference Maydet and Russel27}

Within this approximation, the nonequilibrium concentration of point defects results from the resolution of a set of p-coupled nonlinear autonomous (the time does not directly appear in their right hand side) differential equations ^{Reference Sizmann28} :

where U is the concentration field vector (C _i, C _v, C _s) ^t . This set of equations typically describes the evolution of dissipative dynamical systems. ^{Reference Berge, Pomeau and Vidal29} The energy of the incident particles bumping the material is dissipated via the creation of point and extended defects. The flux ϕ and the temperature T can be understood as control parameters in the theory of dynamic systems. ^{Reference Berge, Pomeau and Vidal29} In the functional F(U, T, ϕ), the physics associated with three main causes of the microstructure formation under radiation damage, i.e., interaction between clusters or loops, fast interstitial migration, and interstitial loops punched out from over pressure gas bubbles must be physically sounded. Modeling these effects implies that many parameters like the bias strengths of bubbles to capture radiation-induced point defects must be known. This is the main difficulty of these approaches. ^{Reference Li, Hu and Stan2,Reference Hu, Burkes, Lavender, Senor, Setyawan and Xu5} Moreover, finding precise solutions for the equations defining the dynamical system [Eq. (5)] is often hopeless and one prefers to study the long term behavior of solutions. As this system is dissipative, it settles down to steady states defined by all possible attractors in the phase space. The number and the nature (attractive, repulsive…) of these attractors determine the long-term behavior of the system versus the initial conditions and the control parameters (the flux, the temperature, and input parameters like the activation energy). The qualitative evolution of the system can be obtained from the phase portrait of dimension p associated to a set of p autonomous differential equations. The phase portrait results from the direct calculation of the velocity field ${{\dot{\bf U}}}$ versus U via Eq. (5). Fortunately, a great deal of insight into the behavior of the system is obtained simply by knowing the sign of the time derivative of each point in the phase portrait. ^{Reference Somorjai30} This numerical procedure allows us to determine all possible attractors.

Figure 3 displays a phase portrait for radiation-induced formation of incoherent precipitates. ^{Reference Maydet and Russel27} In this example, the time evolution of coherent precipitates is dictated by two distinct variables, the number of solute atoms n, and the vacancy concentration x. Precipitates with atomic volumes greater than the matrix volume (n > n*) will increase to relieve the local strain energy, and the undersized precipitates (n < n*) associated with smaller atomic volume will vanish. Within the RT approximation, these continuum fields are given by ^{Reference Maydet and Russel27}

where λ_X, λ_V, and λ_V are accretion rates of solute and they are proportional to the flux, vacancies, and interstitials, respectively. The terms μ_X and μ_V can be written as functions of the temperature and the free enthalpy ΔG(n, x). ^{Reference Maydet and Russel27} They are the rate of loss of solutes and vacancies. All the underlying physics describing the radiation-induced precipitation are included in these terms. This set of equation can be written in a concise way by Eq. (5).

FIG. 3.

2D phase portrait of the number of solute atoms (x) and excess vacancies (n) calculated as the excess between the atomic volume of a host atom and a solute atom in an incoherent precipitate produced under irradiation. ^{Reference Maydet and Russel27} Arrows indicate particle velocities in the (x, n) space. These velocities exhibit different orientations in each quadrant (I to IV). The particle moves into the stable node (0, 0) the saddle critical point (x*, n*). They spiral into the stable focus (non-null real part of the complex conjugate eigenvalues) and circle around the center (null real part of the complex conjugate eigenvalues).

As clearly displayed in Fig. 3, the trajectories of the variables (n, u) are attracted (or repelled by) to different specific points named as attractors in different quadrants (basins of attraction). These points are obviously zeros of F(U, T, ϕ) = 0. This example highlights the interest of the phase portrait to determine the different steady states versus control parameters like the temperature and the irradiation flux.

Different questions naturally appear: how can we define the nature of different possible attractors U ₀ leading to distinct steady states and can a change of control parameters induced a bifurcation from one attractor to another one? The next section is devoted to answer such questions.

C. Bifurcation and instability

For dissipative systems, it is possible to define the nature of these different attractors calculating eigenvalues λ _j of the first derivative DF(U, T, ϕ) (a p × p matrix) of F(U, T, ϕ) calculated for different U ₀ points. As F is a real function, eigenvalues of DF(U, T, ϕ) are real or complex conjugated:

1. (i) A fixed point is called a Node when both eigenvalues are real and of the same sign. The node is stable when the eigenvalues are negative and unstable when they are positive. For the stable node, the eigenvalue(s) with minimal absolute value of the real part is called principle or leading; when the eigenvalues are different, all orbits but two tend to the node along the leading eigen-vector (the picture is reversed for the unstable node);

2. (ii) A fixed point is called a Saddle when eigenvalues are real and of opposite signs. The saddle is always unstable;

3. (iii) A fixed point is called a Focus (sometimes called spiral point) when eigenvalues are complex conjugate; the focus is stable when the eigenvalues have negative real part and unstable when they have positive real part.

It must be noticed that nodes change stability when time is reversed (i.e., when t is replaced by −t), whereas saddles are unstable regardless of the direction of time. Other kinds of fixed points (non hyperbolic ones) can exist. ^{Reference Somorjai30} Obviously, the nature of different fixed points depends on the dimension of the ordinary differential equation set. Increasing the dimension p, various types of hyperbolic attractors can appear. ^{Reference Somorjai30}

Since evolution of the system is dictated by nonlinear differential equations [see Eq. (5)], the nature of steady states may evolve for different values of the set of the control parameters (ϕ, T), labeled μ below, leading to bifurcation of different steady states for critical values of μ_c. A bifurcation leads to an exchange of the stability between two different steady states of the dynamical system. In fact, the system jumps from one attractor to another point attractor. A classification of all possible bifurcations has been dressed by physicists. ^{Reference Cross and Hohenberg31} As an illustration, we consider two types of bifurcation extensively observed in the phase transition community. ^{Reference Tolédano and Dmitriev7} In the 1D system, the time evolution of the OP U during a second order phase transition is given by

Hence, the fixed point is null for μ ≤ 0 and becomes unstable for positive values of μ while nonzero fixed points exist and are stable for positive μ values. The bifurcation is supercritical and known as a pitchfork bifurcation. Introducing a fifth order term U ⁵ in Eq. (7), it is possible to describe first order phase transition:

The null fixed points stable for $\mu \ge {{ - {\nu ^2}} \over 4}$ become unstable leading to a subcritical bifurcation showing a hysteresis as expected for a first order phase transition as displayed in the bifurcation diagram (Fig. 4).

FIG. 4.

Evolution of the fixed point U ₀ versus the control parameter μ in 1D. (a) Displays a supercritical bifurcation and (b) displays a subcritical bifurcation (ν = 2). Full curves are associated with stable solutions and dashed lines define unstable solutions. In the supercritical bifurcation, two opposite stable solutions ${U_0} = \sqrt \mu$ exist for μ > 0. The null solution is stable for μ < 0. In the subcritical bifurcation, two distinct points U ₀ = 0 and ${U_0} = \sqrt {{{\mu + \sqrt {{\mu ^2} + 4\mu } } \over 2}}$ are stable for $- {{{\nu ^2}} \over 4} \le \mu \le 0$ ensuring that two distinct steady states occur for these values of μ. The existence of these two stable points ensures that a hysteresis takes place in this domain.

The stability of a given attractor versus μ can be determined from the amplification rate ω(μ) defined as $\det \left( {{{\left. {{{\partial {F_i}\left( {{{\bf{U}}_{\bf{0}}},\mu } \right)} \over {\partial {U_j}}}} \right|}_{{U_0}}} - \omega \left( \mu \right){\delta _{ij}}} \right) = 0$ resulting from the linearization of Eq. (5) in the vicinity of μ_c. When the real part of ω(μ) vanishes, such a bifurcation is called a subcritical bifurcation diagram. Sometimes, the real part of ω(μ) vanishes but its imaginary part is non-null leading to a limit cycle. This kind of bifurcation is named an Andronov–Hopf bifurcation. ^{Reference Berge, Pomeau and Vidal29,Reference Somorjai30} The main role of nonlinear terms is to introduce a saturation of the exponential growth rate ω(μ). These nonlinearities will select different symmetry adapted points and then induce a change in the symmetry of steady states. This last remark highlights the link between the theory of dynamical systems and the theory of phase transition.

A. Instabilities: hard and soft modes

Instabilities are associated with the maximum of the real part of the amplification rate ${{{\rm{d}}{\mathop{\rm Re}\nolimits} \left[ {\omega \left( {k,\mu } \right)} \right]} \over {{\rm{d}}k}}$ . This maximum is reached for a critical value of the wave vector k _c. Depending on values of Im[ω(k _c , μ _c )], different instability modes can propagate. For k _c = 0, soft modes propagate by analogy with the theory of phase transition. For non-null k _c values, hard modes propagate and the different directions of the propagation are given by all k _c nonequivalents to the space group associated with the invariant of the initial kinetic equation. The set of such critical wave vector is named the star of k _c. ^{Reference Tolédano and Dmitriev7} Steady states can then be defined by different values of Im[ω(k _c , μ _c )]:

1. (i) The critical wave vector associated with the maximum of the real part of the amplification rate can be null and Im[ω(k _c , μ _c )] = 0. A stationary state is produced and mainly dominated by long wave instabilities. No pattern is created.

2. (ii) The critical wave vector associated with the maximum of the real part of the amplification rate can be null and Im[ω(k _c , μ _c )] is non-null. A periodic temporal state without spatial patterning is produced (limit cycle). For non-null values of Im[ω(k _c , μ _c )], the propagation front of the instability is thus defined by $v = {{{\mathop{\rm Im}\nolimits} \left[ {\omega \left( {{k_{\rm{c}}},{\mu _{\rm{c}}}} \right)} \right]} \over {{k_{\rm{c}}}}}$ .

3. (iii) For non-null k _c values and Im[ω(k _c , μ _c )] = 0, a stationary spatial pattern is created.

4. (iv) For non-null k _c and Im[ω(k _c , μ _c )] = 0 values, a double spatial and temporal pattern can be created.

The interplay between reaction diffusion processes is sufficient to explain spatial patterns and was first studied by Turing. Patterns produced do not depend on the detail of the kinetics under investigation but rather on the general mathematical expression for the kinetic equations. In this sense, these patterns exhibit a class of universality. As an example of radiation-induced Turing’s patterns, Fig. 6 displays the formation of a BCC void superlattice in Mo induced by 2-MeV N ions (at 1143 K and 100 dPa). ^{Reference Sikka and Moteff32} Similar patterns have been observed in ceramic materials irradiated by He ions. ^{Reference Beauford, Vallet, Nicolai and Bardot33}

FIG. 6.

Electron micrograph displaying an ordered array of voids in W irradiated at 823 K in a fast neutron reactor at a fluence of 1022 cm⁻². Reprinted with permission from Ref. Reference Sikka and Moteff32, Copyright 4218680185280, AIP Publishing LLC.

The linear analysis cannot by itself determine what type of microstructure is selected beyond the instability, i.e., for values of μ far from the critical control parameter value μ_c. Pattern selection is effectively determined by nonlinear terms. ^{Reference Walgraef34} Fortunately, close to the instability threshold, space-time separation takes place between stable and unstable nodes. Stable modes can thus be adiabatically suppressed and the dynamics reduces to the nonlinear unstable mode dynamics. Only these modes govern the long time evolution of the system and capture asymptotic properties of the complete kinetic model in the vicinity of its bifurcation point. Performing such an adiabatic elimination, the stability of steady states of the system is dictated by amplitude equations.

B. Amplitude equation

A k range exists for which all infinitesimal perturbations can grow in the vicinity of ε = 0 and ρ_L(x, t) can be written as $\sqrt \varepsilon {{\rm{e}}^{i{k_0}x}}A\left( {x\sqrt {{\varepsilon \over \lambda }},t\varepsilon } \right)$ , where A is an unknown function. From this analysis, new space scale and time scale emerge associated with $X = {x \over {\sqrt {{\varepsilon \over \lambda }} }}$ and $T = {t \over \varepsilon }$ .

For example, in the vicinity of ε_c = 0 (more generally in the neighborhood of the μ_c), ρ_L(x, t) exhibits a patterning of wave length $\lambda = {{2\pi } \over {{k_0}}}$ and its amplitude A(X, T) varies very slowly and satisfies simpler equations known as amplitude equations. Performing a scaling with the new variable X, T, the evolution of A(X, T) casts as ^{Reference Cross and Hohenberg31}

This amplitude equation has the structure of a Landau–Ginzburg equation (due to the |A|² term) analyzed in detail in the literature. ^{Reference Cross and Hohenberg31} As a result, Eq. (12) describing spatial inhomogeneities of loop distribution (or void densities) has the generic structure. It is known to produce stripes and hexagonal structures in 2D systems and BCC lattices or planar wall structures in 3D systems (see Fig. 6). Hexagonal and BCC structures appear via subcritical bifurcation while stripes and walls appear via a supercritical bifurcation. ^{Reference Ghoniem and Walgraef9} An expected sequence of selected patterns should develop increasing the irradiation flux from BCC as displayed in Fig. 6 to arrays of dislocations ^{Reference Ghoniem and Walgraef9} as displayed in Fig. 1.

Steady states given by a modulated wave of wave length ${{2\pi } \over {{k_0}}}$ and amplitude ${A_k} = {\left( {\varepsilon - 4k_0^2{k^2}} \right)^{{1 \over 2}}}{{\rm{e}}^{i\psi }}$ , where ψ is an arbitrary phase for ε > 0. The bifurcation from a uniform to a modulated structure generates a band structure made of one pair of opposite wave vectors (k, −k) of modulus k ₀. The final structure thus results from the growth of an arbitrary number of wave–vector pairs of arbitrary orientations leading to stripes with different orientations.

C. Eckhaus and zigzag instabilities

As the patterns under consideration spontaneously break the continuous invariant, different microstructures (Goldstone modes) characterized by their phase ψ(X, t) are expected to exist. These solutions exist in a domain of the parameter space limited by the marginal stability curve defined by $\omega \left( k \right) = \varepsilon - 4k_0^2{k^2} = 0$ . For a 2D pattern, the stability of different modulated phases is dictated by the evolution of a perturbation f(X, T)e ^{i(kx+ψ(X,T))} according to

where D _X and D _Y are the function of k ₀ and ε.

The amplitude perturbation f(X, Y, T) exhibits a relaxation and the phase ψ(X, Y, T) follows a diffusive process occurring over a larger time scale. ^{Reference Pomeau and Manneville35} Far from the bifurcation point ε = 0, the stability of patterns can thus be easily studied using this linear diffusive equation:

1. (i) D _x > 0 and D _y > 0 imply that the structure is stable. A perturbation relaxes via a diffusive process.

2. (ii) D _x < 0, stripes relax to a stable structure via compression or dilatation parallel to walls or rolls according to the exact nature of the system. This instability is known as the Eckhaus instability.

3. (iii) D _y < 0, the structure relaxes via modulations perpendicular to the structure stripes. This instability is named zigzag instability.

This stability analysis of the phase equation points out that transitions between an unstable “stripes” steady state to a stable hexagonal lattice ^{Reference Ghoniem and Walgraef9} can be induced by a set of unstable modes lying in portions of the k − ε plane (Fig. 7). As an example, Fig. 8 highlights this transition in 2D simulations.

FIG. 7.

Stability diagram of stripe structures (2D) resulting from Eq. (11). The structure is unstable with respect to infinitesimal perturbation below the marginal curve (full line). Above the marginal curve (dotted line), Eckhaus (full hatched domain associated with negative D _x values) and zigzag (dashed domain associated with negative D _y values) perturbations can develop destabilizing the pattern. A stable domain insensitive to any perturbation exists leading to roll patterns.

FIG. 8.

To illustrate the effects of instabilities on pattern formation, the transition from an unstable stripe steady state to a stable hexagonal lattice via the propagation of Eckhaus (first row) and zigzag instabilities (second row) is displayed from SH simulations. Reprinted with permission from Instabilities and self-organization in materials Volume I: Fundamentals of nanoscience by Ghoniem (2008) Figs. 9.12 by permission of Oxford University Press.

Such an analysis can be straightforwardly applied to more complex structures like honeycomb lattices for which wave vectors k _j of module k ₀ defining the structure satisfy some geometric relations (k ₁ + k ₂ + k ₃ = 0 leading to 2D honeycomb or 3D BCC structures). These structures resulting from the propagation of unstable critical modes defined by different k _j do not follow the geometrical constraints. Patterns with modulations separated by angles slightly different from ${\pi \over 3}$ can appear for some positive values of ε and k values. The final structure is no more made by hexagons but is similar to rhombuses.

D. Topological defects

This part is largely inspired from the work of Walgraef. ^{Reference Walgraef34} Although each individual spatiotemporal pattern resulting from the propagation of nonstable modes break transitional and rotational symmetries, different patterns grow in different parts of the system implying the formation of various defects connecting these different steady states. The existence of such structures generalizes the notion of topological defects defined for equilibrium structures and is generally the rule beyond pattern forming instabilities as it is the case in equilibrium phase transitions. As an illustration, Fig. 9 displays a penta–hepta pair of cells associated with a dislocation in hexagonal patterns. ^{Reference Ghoniem and Walgraef9}

FIG. 9.

Coarsening of hexagonal patterns in different directions leads to the formation of penta–hepta pairs of cells (black lines) in a hexagonal pattern. This topological defect can be understood as a dislocation. Reprinted from Instabilities and self-organization in materials Volume I: Fundamentals of nanoscience by Ghoniem (2008) Figs. 9.15 by permission of Oxford University Press.

This example points out that topological defects (antiphase boundary, disinclination, dislocation…) appear in pattern-forming systems. To clarify this point, reaction-diffusion dynamics is considered:

In the limit α → 0, the system reduces to a set of discrete equations where each point of the field U(x, t) evolves independently of its neighbors. Complex dynamics of the system may converge to different attractors such as in the case of Andronov–Hopf bifurcation for example. The basic idea ^{Reference Coullet, Emilsson and Plaza36} is that the domains of attraction of disjoint attracting set do not intersect. Considering a system with two attracting sets, one part of the system lies in the domain of attraction of one set while another part lies in the second domain. The evolution of the system leads to the formation of a singularity (for α = 0) lying on a nonattractive set and forms the core of defects. The dimension of the defect core is the dimension of the embedded space d minus the number of unstable directions of the nonattracting set labeled n _u, which can be understood as the co-dimension of the topological defect. This analysis tells us that co-dimension 1 defects correspond to points in 1D, to lines in 2D, and to planes in 3D while co-dimension 2 defect corresponds to points in 2D and filaments in 3D. When diffusion is turn on (non-null α values), the topological defects depend specifically on the system under consideration. If the domains of attraction of two different sets are identical, the diffusion effect is to smooth out the singularities. A typical example of this effect is the basic kink defect. For a stable limit cycle induced by an Andronov–Hopf bifurcation, n _u = 2 and the formation of radial equiphases are created in 2D. The effect of diffusion is to regularize the equiphase forming a regular vortex observed in the XY model of ferromagnets or superfluids. ^{Reference Walgraef34} Obvious generalizations of these vortexes in irradiated materials are the dislocations of layered or hexagonal patterns.

The key point of this analysis is that both a formation energy $\left( {F_i^d} \right)$ and a topological charge (m _i ) can be computed for each type of defects. Since this energy is nonnegligible, the probability for the nucleation of isolated defects in large systems is vanishingly small. However, a finite density of defects (dislocations pairs or dislocation clusters) of zero topological charges $\left( {\sum {{m_i} = 0} } \right)$ could nucleate in the structure leading to perfect array of defects explaining the formation of patterns as displayed in Figs. 1 and 2.

A. Gradient and nongradient dynamics

Assuming that the evolution of a system can be modeled within the PF framework, the time variation of the LRO η(x, t) is assumed to be linked to the variation of an effective real free energy functional,

Here, F was separated in two distinct terms. The first one F _R(η(x, t), μ) is associated with the relaxation part $- \Gamma {{\delta {\cal L}\left[ \eta \right]} \over {\delta \eta \left( {x,t} \right)}}$ of the dynamics (Γ being a positive operator). F _NR(η(x, t), μ) describes all nonrelaxing processes associated with the dynamics and is orthogonal to F _R(η(x, t), μ) for the given inner product. ^{Reference Ghoniem and Walgraef9}

Different cases exist:

1. (i) If F _NR = 0 and Γ is constant, the evolution of the system follows the line of steepest descent of ${\cal L}$ . The system is named a relaxation gradient flow.

2. (ii) If F _NR = 0 and Γ is not constant (system with conservation laws), transient patterns do not correspond to the lines of steepest descent of ${\cal L}$ . The system is a relaxation nongradient flow.

3. (iii) If F _NR is non-null, the status of the dynamics is still under investigation. The system is a nonrelaxation potential flow.

For systems described by relaxation flows, their behavior should be similar to those at equilibrium since a Lyapounov functional ${\cal L}\left[ \eta \right]$ can be determined. This functional plays the same role than the free energy and then can be understood as an “effective” free energy for dynamical systems. In particular, spatiotemporal patterns resulting from the minimization of ${\cal L}\left[ \eta \right]$ should appear and generate Turing-like patterns leading to the formation of topological defects. The next section is devoted to highlight this point discussing the nonequilibrium behavior of radiation-induced decomposition in nonmiscible alloys.

B. Application of the PF method: spinodal decomposition of alloys under irradiation

The spontaneous formation of patterns under irradiation was observed in many nonmiscible alloys for seventy years. The observed experimental patterns result from the existence of steady states. The main ingredient leading to steady states results from the competition between a thermally activated dynamics driving the system at temperature T to the equilibrium and an externally imposed particle exchange Γ^ext(ϕ) regardless of the chemical identity of particles. This last dynamics is by essence athermal (T = ∞) and driven by the nature, the energy, and the flux ϕ of impinging particles. Martin ^{Reference Martin38} points out that the ordering induced under irradiation is indistinguishable from that of an equilibrium Ising model assuming that this exchange has an infinite range as expected. ^{Reference Bergersen and Racz39} The connection between these dynamical systems and equilibrium system is performed via the construction of a Lyapounov functional acting as an effective free energy. Thus, the dynamics of order–disorder transformation under irradiation for binary alloys with positive heat of mixing is indistinguishable to an equilibrium Cahn–Hilliard (CH) equation, the continuum counterpart of the Ising model, with infinite-range interactions. By analogy with the CH equation including elastic energy assumed to be of infinite range, ^{Reference Khatchaturyan6,Reference Huston, Cahn and Hilliard40} a “law of corresponding states” can be introduced for a system at temperature T via T _eff = T(1 + W). T _eff can be understood as the effective temperature of the system and $W = {{{\Gamma ^{{\rm{ext}}}}\left( \phi \right)} \over {\Gamma \left( {T,\phi } \right)}}$ (Γ(T, ϕ) is the thermally activated frequency of atoms promoting equilibrium) is a kinetic enhancement factor induced by the infinite range exchanges. From this analysis, T _eff becomes important at very low irradiation temperatures and do not affect the stability of precipitates at high temperature and low doses in agreement with experimental results. ^{Reference Adda, Beyeler and Brebec41} Since athermal exchanges are assumed to take place over an infinite range, such an approach is unable to describe patterns associated with spatial modulation of the chemical composition of the alloy at the mesoscopic scale, i.e., few tens of nanometers, experimentally observed under irradiation. ^{Reference Tsaur, Lau and Mayer42} To overcome this difficulty, the exchange of atoms occurring during the irradiation process is assumed to occur over a finite range with an average relocation distance R by analogy with phase separation induced by long range Coulomb interactions describing similar patterns. ^{Reference Chen and Khatchaturyan43}

At equilibrium, the local composition (atomic fraction) given by the OP η(x, t) = c(x, t) − c ₀, null in the disordered high temperature phase of uniform concentration c ₀ in the ${A_{\bar{c}}}{B_{1 - \bar{c}}}$ alloy of average composition $\bar{c}$ is dictated by the CH equation routinely used to study the phase separation dynamics within the PF framework:

where f(η(x, t)) is the free energy density which usually results from a Landau polynomial expansion $f\left( {\eta \left( {{\bf{x}},t} \right)} \right) = {{{a_2}} \over 2}\eta {\left( {{\bf{x}},t} \right)^2} + {{{a_3}} \over 3}\eta {\left( {{\bf{x}},t} \right)^3} + {{{a_4}} \over 4}\eta {\left( {{\bf{x}},t} \right)^4}$ , the parameters of which are fitted to satisfactorily reproduce the equilibrium phases. ^{Reference Demange, Luneville, Pontikis and Simeone44} In this equation, Γ(T, ϕ) is the averaged mobility of species assumed to be concentration independent. ^{Reference Bray45}

On the other hand, the contribution of irradiation to the kinetic evolution of the alloy composition entering Eq. (16) can be expressed by the following equation ^{Reference Sigmund and Gras-Marti22,Reference Simeone and Luneville23} :

where p _R(x) ∝ exp(−|x|/R) is the density probability function describing the relocation of atoms set in motion by elastic collisions. R can be understood as a mean free path describing the relocation of atoms and is roughly equal to 1 nanometer. ^{Reference Luneville, Mallick, Pontikis and Simeone46}

By analogy with the studies of long range interactions at equilibrium, ^{Reference Chen and Khatchaturyan43,Reference Ohta and Kawasaki47} the radiation-induced dynamics can be written as

where G(x, y) = G(|x − y|) expresses the spatial dependence of the radiation-induced relocation and can be understood as an effective long range interaction. This function G(x) is related with p _R(x) via ∇² G(x) = p _R(x) − δ(x). The total dynamics resulting from the evolution in parallel to the thermal and the irradiation dynamics can thus result from the minimization of a Lyapounov functional ${\cal L}\left[ \eta \right] = F\left( {\eta \left( {x,t} \right)} \right) + {W \over 2}\left( {\int {{\rm{d}}x{\rm{d}}y\eta \left( {x,t} \right)G\left( {x,y} \right)\eta \left( {y,t} \right)} } \right)$ . This Lyapounov functional can be seen as an effective potential. Its minimization under the constraint of species conservation gives all possible steady states produced under irradiation. This effective free energy captures the physics associated with radiation damage in this particular system assessing the use of the PF method to handle the formation of the radiation-induced microstructure in nonmiscible alloys. The term W(ϕ, T) highlights the competition between relocation of atoms and their thermal diffusion.

Figure 10 displays the 2D long-term evolution of η(x, t) resulting from a PF simulation. The power spectrum of the microstructure, i.e., the structure factor of η(x, t) resulting from the average of square of the Fourier transform of η(x, t) over initial configurations, exhibits a ring shape for long times revealing the existence of wave vector modulations of modulus k ₀ assessing pattern formation.

FIG. 10.

(a) 2D numerical simulation of the long-time microstructure (1000 in reduced units) of the AgCu alloy irradiated by 1 MeV Kr ions. Black (white) domains correspond to Cu rich (poor) regions. Boundary conditions imply that the calculated microstructure exhibits a point symmetry and two mirror symmetries. (b) Power spectrum of the microstructure. This power spectrum, i.e., the structure factor of the OP exhibits a Mexican hat shape. The maximum of its intensity exists for a well-defined modulus ${k_0} = \sqrt {k_x^2 + k_y^2}$ of wave vectors as displayed in figure (c).

A careful inspection of ${\cal L}\left[ \eta \right]$ clearly shown that order–disorder transitions (k ₀ = 0) or a decomposition process occurs (non null k ₀ values) for different (Wa, R) couples as shown in Fig. 11. From this analysis, it also becomes possible to compute the evolution of k ₀ versus R for different W values as shown in Fig. 11. This evolution extracted from numerical PF simulation shows a clear linear evolution of k ₀ with R ⁻¹ in agreement with theoretical investigations resulting from the minimization of ${\cal L}\left[ \eta \right]$ . ^{Reference Simeone, Demange and Luneville48} This pattern formation is assessed by the peaked shape of the power spectrum (inset b in Fig. 10). The variation of k ₀ with the parameters R and W points out a supercritical bifurcation of the dynamical system for the control parameter $\mu = \sqrt {W{R^4}} - 1$ (inset c in Fig. 10). Below μ_c = 0, relocation of atoms are not sufficient to counterbalance the thermal diffusion and the spinodal decomposition takes place under irradiation. For μ_c > 0, a balance between the two competitive dynamics occurs leading to the formation of modulated structures of wave vectors (k) of non-null modulus k ₀ minimizing ${\cal L}\left[ \eta \right]$ .

FIG. 11.

(a) Evolution of k ₀ versus R for different W values (black squares: W = 0.04, open circles: W = 0.7) and (b) k ₀ values are extracted from the maximum of the universal structure factor obtained from simulations. (c) The evolution of k ₀ versus the control parameter $\mu = \sqrt {W{R^4}} - 1$ displays a supercritical bifurcation.

Patterns produced under irradiation result from slow composition modulations in space and time of simple basic patterns defined by the wave vector of modulus k ₀ minimizing ${\cal L}\left( t \right)$ . Although each individual stationary pattern breaks transitional and rotational symmetries (Fig. 12), different patterns grow in different parts of the system implying the formation of various defects connecting these different steady states ^{Reference Ghoniem and Walgraef9} as displayed in Fig. 12. These patterns depend not only on R and W but also on the average composition of the nonmiscible alloys.

FIG. 12.

2D phase diagram for steady states obtained from ${\cal L}\left( t \right)$ in the one-mode approximation. Depending on the averaged composition $\bar{\psi }$ (in reduced units), stripes, bubbles in hexagonal networks and constant phase can be created. Hatched areas correspond to co-existence domains. This phase diagram is validated by numerical simulations (microstructure drawings). All simulations were performed for Δ = 0.2, R = 3, t = 3000 and the domain size equal to 200 in reduced units.