Properties of spatial arrangement of V-type defects in irradiated materials: 3D-modelling

We consider the dynamics of pattern formation in a system of point defects under sustained irradiation within the framework of the rate theory. In our study we generalize the standard approach taking into account a production of defects by elastic fields and a stochastic production representing internal multiplicative noise. Using 3D-modelling we have shown that with the damage rate growth, a morphology of clusters composed of vacancies changes. The same effect is observed with variation in the multiplicative noise intensity. Stationary patterns are studied by means of correlation analysis.


Introduction
It is well known that metals and alloys under irradiation are typical examples of far-from-equilibrium systems with nonlinear coupling between their elements. The simplest ones are vacancies (V-defects) and interstitials (I-defects). Depending on irradiation conditions related to both defects replacement rate and temperature, the point defects can arrange into defects of higher dimension such as clusters, defect walls with vacancy loops [1,2], voids [3], precipitates [4], bubble lattices [5][6][7]. At the same time, irradiation can modify a surface of targets due to sputtering resulting in the formation of ordered nano-size surface structures widely used in modern electronics [8][9][10]. Such patterns are considered to be dissipative structures. Their theoretical and experimental investigations are widely discussed in literature during last four decades (see, for example, [7,[11][12][13][14][15]). The production of defects and their microstructure can principally modify the physical and mechanical properties of an irradiated material leading to swelling, hardening, etc. As far as defect clusters are responsible for material embrittlement [16,17], a study of microstructure of materials, defects and their arrangement into clusters is an urgent problem in modern material science and in statistical physics. In the present paper we consider the processes of the formation of V-type defect structures using both analytical treatment and numerical modelling.
In theoretical investigations of microstructure formation, a formalism known as multiscale modelling serves as an effective tool. In the framework of this approach, properties of the system on lower (microscopic) hierarchical level are self-consistently taken into account in the modelling procedures at higher levels of description (mesoscopic or macroscopic). On the other hand, one can use the so-called "hybrid" methods such as phase field crystals method where the dynamics of microscopic crystal system is described in terms of atomic densities covering time scales for elastic interactions between atoms and their diffusion [18,19]. It can be used to study the diffusion of defects [20], dislocation dynamics [21,22] and structural transitions [23,24]. At a mesoscopic hierarchical level, one can consider the dynamics of defects using the rate theory where concentration of both interstitials and vacancies are considered (see, for example, references [11,13,15,25]). Usually such models take into account irradiation production rate and interactions between defects except defect production caused by elastic deformation of a medium. An idea of stochastic production of defects during irradiation was exploited to study the ordering, chemical patterning and phase separation of irradiated materials [26][27][28][29][30][31][32].
The main goal of our paper is to consider the dynamics and spatial arrangement of an ensemble of point defects using the rate theory describing the evolution of point defects. We take into account the production of defects by irradiation, relaxation (the effect of sinks) and the production of defects by elastic deformation of a medium caused by point defects. We generalize the standard approach [33] considering a system with two spatial scales related to the diffusion and microscopic interactions. Within the framework of the previously derived approach used to study nanosize patterns in stochastic systems (see references [34][35][36][37]) we consider the most probable stationary defect structures. In numerical investigations we solve a three-dimensional problem and discuss the effect of the system parameters onto the dynamics and morphology change of defect structures. We describe stationary patterns by means of correlation analysis.
The paper is organized as follows. A generalized model described by the main mechanisms of the formation of defects with two spatial scales related to diffusion scale and defect interaction scale is presented in section 2. The analysis of stationary states and their stability is provided in section 3. In section 4 we numerically study the pattern formation where we discuss the dynamics of statistical averages and analyze stationary patterns. We conclude in section 5.

Model and basic equations
Following the standard approach [39] the dynamics of point defects is described by a two-component model Here, c v corresponds to populations of vacancies (c 0 v is the equilibrium vacancy concentration) and c i relates to interstitials. The first terms in equation (1) relate to the displacement damage rate and take into account a production of defects due to irradiation. The second terms describe the effect of sinks (S i and S v ) related to bias factors Z i,v , network dislocation density ρ N , vacancy loops ρ v , and interstitial loop densities ρ i as follows: Diffusivities of vacancies and interstitials are defined in the standard manner: D {v,i} = D 0 {v,i} e −E m{v,i} /T , where E m{v,i} is the migration energy for vacancies and interstitials, respectively, T is the temperature. The last terms govern the nonlinear contribution caused by point defects annihilation with the recombination coefficient α = 4πr 0 (D i +D v )/Ω given in terms of recombination radius r 0 and atomic volume Ω.
In metallic systems due to the difference between migration energies of point I-and V-defects (for example, for pure nickel E mv = 1.04 eV, E mi = 0.3 eV with r 0 = 1.5 · 10 −9 m) there is a difference between their diffusivities. It allows one to introduce a small parameter ν ≡ D v /D i 1. In such a case, a renormalization of time scales allows us to put ν∂ t c i 0 and eliminate the adiabatically fast variable c i from the reduced second equation of the system (1). Next, rewriting , where x relates to the vacancy concentration. Here, the last term is related to nonlinearity caused by interaction of interstitials and vacancies governing the corresponding annihilation processes. In our previous study (see reference [25]) we neglected the effect of interstitials assuming their fast dynamics and motion to sinks. An allowance for interstitials in the following derivation of the model does not change the qualitative picture of the microstructure evolution, but this allowance encumbers the physical model. This statement was discussed previously in reference [12].
We have to note that the obtained one-component model does not incorporate the production of defects by an elastic field caused by the presence of defects. Following reference [38] this effect can be considered by introducing an additional source term into dynamical equation for x as e −E /T , where

33001-2
Properties of spatial arrangement of V-type defects Here, φ(r ) is the potential of the deformation field caused by the presence of defects. It facilitates the processes of defect generation by decreasing the activation energy E ; E f is the defect formation energy in the absence of deformation field; r = x −1/3 is the averaged distance between the defects. At r → 0, the field φ has the asymptotics [38]: φ E e 0 r 3 /(1 + r 6 ). Therefore, the source term of defects can be rewritten as follows: G exp[εx/(1 + x 2 )], where ε ≡ 2Z E e 0 /T is defined through the defect formation energy E e 0 0.01 eV and coordination number Z ; the renormalized constant G is proportional to the probability of defect generation by an elastic field caused by the presence of defects. It is determined by means of Debye frequency and atomic volume; it exponentially depends on the relation between the defect migration energy and temperature. The introduction of an additional term means that the deformation field created by defects facilitates the defect generation due to a decrease in the activation energy. This is essential in laser radiation due to temperature instabilities, whereas at particle irradiation its efficiency is small compared to the defect production in cascades. However, in our case, without loss of generality, we retain this term assuming G 1.
As far as point defects are mobile species of a microstructure the corresponding rate equation would include spatial operators defined by introducing a flux of defects. This diffusion flux has an ordinary Fick It can be rewritten in the canonical form J has the density f (x) = x(ln x − 1). Here, the second part relates to pair interactions in a self-consistent manner: [25,36,40]. We assume the attraction potential −ũ(r ) in a symmetric form, i.e., ũ(r )r 2n+1 dr = 0. If the field x(r ) does not vary essentially on the interaction range of defects r 0 Ω 1/3 , then one can use the approximation The first term in equation (4) leads to the standard relation between U and x in the framework of the elasticity theory [33]. Indeed, the effective flux takes the form J = −D(1− 2 κx/T )∇ ∇ ∇x, where κ is the bulk modulus, is the dilatation parameter. The second term in the above expansion is responsible for microscopic properties of defect interactions described by interaction radius r 0 . Under normal conditions, this term is negligible compared to the ordinary diffusion one. However, as far as one has a density de- , it can be negative in some interval for x. It means that a homogeneous distribution of defects, starting with some critical speed of its formation related to the temperature, sinks the density, and the dilatation volume becomes unstable. The emergence of a directional flux of defects results in supersaturation of vacancies and in the formation of voids. From mathematical viewpoint such a divergence appearing at short time scales cannot be compensated by a nonlinear part of the dynamical equation for x. The second term in the expansion (4) can prevent divergencies of the derived model due to microscopic properties of defect interactions. Therefore, the term with the second derivative should be retained. This term governs the typical sizes of defect clusters. Inserting the above expressions into a rate equation for a vacancy concentration we arrive at a deterministic reaction-diffusion equation of the form

33001-3
Here, we use renormalization of a spatial coordinate as r = r/L d and introduce a dimensionless length = r 0 /L d .
Considering real systems we have to take into account the effect of fluctuations or stochastic sources representing the microscopic action onto the system dynamics. In our further study we introduce a random source treated as an internal noise resulting in dissipative dynamics of the whole system. Such a noise should satisfy the fluctuation dissipation relation. Acting in the standard manner, we rewrite the original deterministic model in the form where for the free energy functional F [x] we have Formally, equation (7) can be represented in the canonical form where for the functional U [x] we know only its first derivative, i.e., Following references [34][35][36] we can introduce a fluctuation source obeying the fluctuation-dissipation relation in an ad hoc form: proportional to the bath temperature. Hence, equation (11) treated in the Stratonovich interpretation represents the generalized model considered below.

Stationary states analysis
In our study, the main attention is paid to stationary patterns. Therefore, to describe the stochastic system behavior in the stationary limit we have to obtain a stationary distribution P s [x]. To this end, we need to find a stationary solution of the corresponding Fokker-Planck equation satisfying the Langevin equation (11) [41]. As it was shown previously (see references [35,36,[42][43][44][45][46]) the functional of the stationary distribution of the vacancy concentration field has the explicit form where the effective potential is Here, Σ is a renormalized parameter proportional to σ 2 . It is seen that the second term in the effective potential (13) serves as entropy contribution. It may lead to the so-called entropy-driven phase transitions [42,[45][46][47], phase decomposition [32,44] and patterning [35,36,43].
Firstly, let us consider the stationary states x s of the homogeneous system determined as solutions of the equation The corresponding solutions are shown in figure 1 (a). It is seen that at large defect damage rate the system is always in monostable state, whereas at small K the bistable regime is observed. The last effect  means that at small K the formation of defects by elastic field plays an essential role and leads to bistability of a system. At large K , the main mechanism of defect production relates to defect generation in cascades. At small K , the bistability domain for ε at fixed K is [ε b1 , ε b2 ]. In the plane (K , ε) binodals ε b1 (K ) and ε b2 (K ) form a cusp binding the bistability of the system states [see figure 1 (b)]. Here, in the cusp, the system is bistable, while outside the cusp the system is monostable. Below the cusp, the system is considered to be in a depleted state of defects, while above the cusp, an enriched state of defects is realized. For a homogeneous system, the depleted state can be related to a "crystalline" phase with a small amount of defects, whereas enriched state corresponds to some kind of "amorphous" phase with a large number of defects. Therefore, in the cusp, two possible phases ("crystalline" and "amorphous") can be observed 1 .
As far as the functional of stationary distribution P s [x] is known, the solution of the variational problem δP s [x]/δx = 0 allows one to find stationary structures x s (r) as its solution. Following references [34][35][36][37] the most probable solutions x(r) corresponding to the minima of U ef [x] can be found as solutions of the equation Substituting the necessary expressions we get According to this equation and stationary states behavior we can study the stability of stationary states in linear analysis by introducing small fluctuations δx = x − x s . In such a case, the linearized equation (16) has an exponential solution δx ∝ e [Λ+ω(k)]t , where Λ is the Lyapunov exponent responsible for homogenous perturbations, whereas ω(k) takes care of the stability of stationary states due to inhomogeneous perturbations. It gives dispersion law allowing one to define critical wave-numbers k c as solutions of equation ω(k) = 0. Expressions for Λ and ω(k) are: It follows that the noise action stabilizes the stationary state with respect to homogeneous perturbations.
In figure 1 thin lines bind a domain of inhomogeneous distribution of vacancies: at ω(k) < 0 no patterns 1 To distinguish real crystalline and amorphous states the long range order parameter should be used. In our approach we do not have such a criteria. Therefore, two macroscopic phases can be distinguished only by the value of the population x s of stationary defects.

33001-5
can be realized (below thin curves), in the opposite case dissipative structures of point defects emerge (above thin curves). The analysis of the dispersion relation with respect to three different solutions of equation (14) realized in bistable domain allows one to set that for stable solutions x s , the noise extends the domain of unstable modes whereas in the vicinity of unstable branch x s , the noise shrinks the domain for wave-numbers related to unstable modes. Moreover, here the noise is capable of decreasing the value for the most unstable mode k 0 = k c / 2.

Numerical results
Let us study the spatial arrangement of defect structure by means of computer simulations. To this end, we make discretization of the system in 3-dimensional space with N × N × N sites, where N = 128. Let us consider the behavior of a deterministic system, firstly, setting Σ = 0. Snapshots of the typical evolution of the system at K = 0.075 are shown in figure 2. Here, in figure 2 (a) the total distribution of the vacancy field is shown [blue domains (online) relate to a small vacancy concentration limit, red ones (online) correspond to high concentration limit]. In figure 2 (b) spatial defect structures are shown as a cross-section at a fixed concentration value x = 0.9. Here, the concentration of initially Gaussianly distributed defects is small. During the system evolution, the concentration increases essentially. It is stipulated by the production of defects through irradiation and elastic field in the system having an annealed concentration of defects at previous stages. When supersaturation is reached, the defects start to condense into phases (structural type of grain boundaries) with high concentration of defects. At next  stages of the system evolution one gets grains without defects. From thermodynamical viewpoint such "ideal" grains are unstable and at the next time steps one can observe the formation of additional defects inside the grains. Grains evolve according to Ostwald ripening mechanism. At late stages corresponding to stationary limit one has a net of vacancy loops and vacancy walls with voids. The quantitative picture of the system evolution is shown in figure 3. Here, the averaged concentration 〈x〉 [see figure 3 (a)] increases at small times and after supersaturation it decreases toward stationary value 〈x〉 s due to the formation of defect clusters and motion of defects to sinks. It is seen that with the growth of the defect production rate, an average concentration 〈x〉 takes up elevated stationary values. The principal information about the ordering process is provided by the dispersion behavior 〈(δx) 2 〉 [see figure 3 (b)]. This quantity plays the role of an order parameter in phase transition, phase separation and pattern formation processes [29,31,36,48]. An increasing dynamics of this quantity means the ordering of the system with the formation of different phases (in our case, phases with low and high vacancy concentrations). If it attains a non-zero stationary value, then an ordered phase organizes. In our case it is seen that 〈(δx) 2 〉 increases toward maximal value corresponding to the formation of well-structured grains: if grains have a small number of defects inside, then the order parameter takes a large maximal value. At late stages, due to reconstruction of defect clusters, it decreases toward its stationary non-zero values. It is seen that with the growth of defect production rate, the reconstruction of clusters is faster.
Moreover, at elevated K , the order parameter takes up larger values meaning the formation of well organized structure of point defect clusters. The corresponding stationary snapshots [see figure 3 (c)] illustrate different structures of defects at different values for K . The dynamics of spatial arrangement of defects can be seen from the structure function S(k, t ) time behavior as the Fourier transformation of the twopoint correlation function C (r, t ) ≡ 〈δx(r, t )δx(0, t )〉. The spherical analogue of S(k) is shown in figure 4. It is seen that during the system evolution, the peak of S(k) that denotes the period of spatial structures increases and shifts toward large wave-numbers k meaning the formation of a structure with a smaller period compared to structures related to the early stages when supersaturation is reached. An increase in the peak height corresponds to the formation of well-organized phases enriched by V-type defects. More information on spatial arrangement of V-type clusters in the stationary limit is provided by a two-point correlation function C (r, t → ∞). It is a spherical analogue calculated in 3D space at different values for defect production rate shown in figure 5 (a). It follows that C (r ) exponentially decays from its maximal value at the point r = 0. The decaying rate depends on the damage rate K . It is worth noting that with an increase in K , spatial correlation function manifests an oscillatory behavior meaning the formation of ordered structures with a fixed period. To study stationary patterns in detail we use an . The predicted behavior of the period of patterns versus damage rate is in good correspondence with well-known experimental data [49].
Next, let us study the most probable stationary patterns in stochastic case by varying the noise intensity Σ. For a homogeneous system one can find most probable values for the vacancy concentration as solutions to the stationary homogeneous equation (14). It follows that the noise action leads to an increase in the most probable vacancy concentration due to a stochastic mechanism of defect production compared to the deterministic case. Studying the stationary patterns by a numerical solution of the original equation (16) one can find dependencies of 〈x〉 s and 〈(δx) 2 〉 s versus the noise intensity Σ. The corresponding results are shown in figure 6. It is seen that with the growth of noise intensities the averaged concentration increases. The same effect is observed in the dependence 〈(δx) 2 〉 s (Σ). It means that large fluctuations in the production of defects promote the organization of well-ordered structures of V-type defects due to an entropy-driven mechanism of pattern formation [25].
From the naive consideration one can predict a decrease of the period of structures and an increase in their correlation radius. Indeed, analyzing the corresponding structure function shown in figure 7 (a) one can find that with the noise intensity growth the position of its maximum is shifted toward large wavenumbers, i.e., the period of defect structures decreases. A growth of the peak at elevated Σ corresponds to an increase in the correlation radius. The corresponding dependencies of the correlation radius and the period of patterns are shown in figure 7 (b). Snapshots of stationary patterns at different noise intensities and fixed main system parameters are shown in figure 7 (c) in order to show the morphology changes of patterns. It is seen that in the noiseless case one has a large amount of voids and a small number of linear defects. If we introduce the noise, then the effective potential (13)

Discussion
In our study we have considered the pattern formation in systems of point defects under the action of irradiation. We have generalized a standard approach of point defects dynamics by taking into account the production of defects by an elastic field caused by the presence of defects, their interactions and stochastic production of defects as a microscopic effect on the system described at mesoscopic level.
The main result of our paper relates to identifying the effect of irradiation conditions and microscopic processes of defect production and their interactions onto morphology changes of the V-type defect structure. Considering the deterministic picture of pattern formation in the framework of correlation analysis we have shown that the period of patterns decreases with the displacement damage rate growth. It is compared with the diffusion length and according to the data for pure nickel with the net dislocation density ρ N ∼ 10 14 m −2 it takes up values around 10 −7 m. For a cold worked nickel characterized by ρ N ∼ 10 15 m −2 , the period of patterns is around 100 Å. This result was discussed previously by experimental investigations [49]. We have studied the morphology changes of stationary patterns by measuring the correlation radius of stationary structures. It was shown that when the morphology of patterns changes crucially the correlation radius manifests an anomalous behavior. In a deterministic system, this effect is well seen when transition from linear structures toward the formation of planar structures (defect walls or grain boundaries) is realized. Here, the correlation radius takes up large values for pure planar defect structures.
Considering stationary patterns in a stochastic system it was shown that the noise action plays a role similar to the displacement damage rate which increases the number of defects. Due to supersaturation of point defects their organization into clusters of higher dimension can be observed. We have compared stationary patterns for a deterministic system and the corresponding stationary patterns for a stochastic one. Here, the period of patterns decreases with the noise intensity growth. We have shown that due to reconstruction of stationary distribution of the vacancy concentration field caused by the noise effect, the morphology of patterns changes essentially: there is a transition from patterns with voids to patterns with linear defects. Such a morphology change corresponds to a decrease in the correlation radius of spatial structures. With further noise intensity growth, spatial arrangement described by an increasing correlation radius is observed.
We have considered the dynamics of vacancies only assuming fast relaxation of interstitials and neglecting the dynamics of other elements such as di-, tri-, tetra-vacancies and dislocations. The proposed approach can be generalized by taking into account the dynamics of all the above elements. In our study we have used material constants for pure nickel. However, the obtained results, which are quite general, can be applied to the study of an ensemble of point defects in certain material under sustained irradiation.