From the Atomistic to the Mesoscopic Scale Modeling

Diffusion Foundations ISSN: 2296-3642, Vol. 12, pp 111-126 doi:10.4028/www.scientific.net/DF.12.111 © 2017 Trans Tech Publications, Switzerland

Online: 2017-09-08

From the Atomistic to the Mesoscopic Scale Modeling of Phase Transition in Solids H. Zapolsky1,a , G. Demange1,b and R. Abdank-Kozubski2,c 1 GPM,

UMR CNRS 6643, University of Rouen, 76575 Saint Étienne du Rouvray, France

2 Institute a

of Physics, Jagiellonian University in Krakow, Lojasiewicza 11 30-348 Krakow, Poland [email protected], b [email protected], c [email protected]

Keywords: phase transition, self-assembly kinetics, nanostructuration, irradiation effects, atomic density functional, atomic fraton theory, phase field.

Abstract. The phase-field method is a very powerful tool to model the phase transformation and microstructural evolution of solids at mesoscopic scale. However, several important phenomena, like defect formation, grain boundary motion, or reconstructive phase transitions require an atomic scale study. Recently an approach called the quasi-particle approach, based on the Atomic Density Function theory was developed to incorporate the atomic-level crystalline structures into standard continuum theory for pure and multicomponent systems. This review focuses on the description of different computational methods used to model microstructural evolution and self-assembly phenomena at mesoscopic and atomistic scales. Various application examples of these methods are also presented. Introduction A relaxation phenomenon in materials that induces a kinetic process and the temporal evolution of microstructure is a fascinating theme from the viewpoint of non-equilibrium thermodynamics. At the same time, the knowledge of the microstructure is important to determine many properties of materials (mechanical, thermal, electrical, magnetic etc.). The term "microstructure" is very general., and refers to the spatial distribution of micro-sized structural features. The size, shape and spatial arrangement of such features in the microstructure determines the physical and mechanical properties of materials. Computational modelling techniques are now widely employed in material science to study a microstructural evolution in different types of materials. Due to recent advances in computing power and simulation methodologies, modelling can enable the rapid testing of theoretical predictions, or the understanding of complex experimental data at a relatively low cost. In most cases, kinetic phenomena in solids, related to changes in microstructure, occur at various time scales, from a fraction of a second to a year, and evolve along several orders of magnitude in length, starting from the atomic scale to the mesoscale. In the past, several simulation methods were introduced to address such a multiscale phenomenon. In the seventies, the atomic density function (ADF) approach was proposed by A.G. Khachaturyan [1]. A theoretical foundation of the ADF theory is the non-equilibrium Helmholtz free energy of a system that is a functional F (ρ) of an arbitrary distribution of atomic density function, ρ(r, t). This function is an occupation probability to find an atom at the site r of the underlying Ising lattice. The ADF kinetic equation is based on the assumption that the relaxation rate, ∂ρ(r, t)/∂t, where t is the time, is proportional to the transformation driving force δF (ρ)/δρ. Using this theory, the kinetics of order-disorder phase transformations and the microstructure evolution at atomic scale was modelled [2, 3, 4, 5]. Recently, the continuum version of the atomic ADF theory, called the quasi-particle (QA) approach, was proposed to model the atomic-scale evolution of crystalline patterns [6, 7]. This model opens a way to answer numerous outstanding questions concerning the atomistic mechanisms of the formation of defects (dislocations, grain boundaries, etc.), the nucleation in solid-solid transformations, the formation of polymers due to the aggregation of monomers in their solution, the folding and crystallization of polymers, and their responses to external stimuli. All rights reserved. No part of contents of this paper may be reproduced or transmitted in any form or by any means without the written permission of Trans Tech Publications, www.ttp.net. (#100690166-08/09/17,10:12:22)

112

Diffusion Foundations Vol. 12

To describe the microstructural evolution at a mesoscale, the phase field model was extensively used for microstructure simulation studies because of its ability to reproduce complex microstructures, without any a priori assumptions [1, 8]. The phase field method allows arbitrary precipitate morphology, and a manifold of thermodynamic driving forces, including the bulk chemical free energy, the interfacial, magnetic and elastic strain energies, as well as various external applied fields. In this method, each phase or domain in a microstructure is characterized by a set of field variables (compositions and/or order parameters). This allows to avoid the explicit tracking of interfaces, and to alleviate the computation load of simulations. The purpose of the paper is to give a short description of the ADF and phase field approaches, and to illustrate these methods through concrete examples. The Atomic Density Function Theory Discrete ADF theory. The atomic processes during ordering and precipitation can be modelled within the atomic density function theory [1], in which phase fields are represented by the atomic occupation density functions PA (r, t). PA (r, t) is the probability of finding an A atom at a given lattice site r, at time t. In this method, each site of the Ising lattice is occupied by an atom. For example, in the case of a binary substitutional solid solution AB, the density functions (i.e. the occupation probabilities) PA (r, t) and PB (r, t) are related by the identity PA (r) + PB (r) = 1, and are thus not independent. Using this identity, one can describe the atomic configuration of the system using one density function P (r, t) only. A temporal and spatial dependence of P (r, t) describes the evolution of the microstructure on the atomic scale. The temporal evolution of the density function P (r, t) can be evaluated using the microscopic diffusion equation [1]: [ ] δF dP (r, t) ∑ ′ L(r − r ) = , (1) dt δP (r ′ , t) r′ where L(r − r ′ ) is the effective rate of elementary diffusion jump between the nearest lattice sites r and r ′ . The kinetic matrix, L(r − r ′ ), depends on the density field, P (r, t), but during the growth stage, this dependence is reduced to a renormalization of the effective values of L(r − r ′ ). F is a total free energy of the system which is a function of the atomic density function, P (r, t), and can include the different contributions such as chemical, magnetic and elastic. The summation is carried out over all N crystal lattice sites, and the different contributions such as chemical, magnetic and elastic. The equation 1 approximates the evolution rate by the first non-vanishing term of its expansion with respect to the thermodynamic driving force, δF /δP (r, t) (small driving force). It is significantly non-linear with respect to the density field P (r), although it is linear with respect to the driving force. In the framework of the mean field theory, the configurational chemical free energy of a binary system can be expressed by: Fchem =

∑ 1∑ W (r − r ′ )P (r)P (r ′ ) + kB T [P (r) ln(P (r)) + (1 − P (r)) ln(1 − P (r))] , (2) 2 r,r′ r

where W (r−r ′ ) = WAA (r−r ′ )+WBB (r−r ′ )−2WAB (r−r ′ ) is the mixing energy. The summation is carried out over all sites of the Ising lattice. The first term in this expression is a configuration part of the free energy and the second one is the entropy contribution. In the disordered state, the density function P (r, t) is a constant and is equal to c¯A , that is an average density of atoms of sort A in system. Using the static concentration wave approach [1] the density function P (r, t) of ordered phase generated by the wave vectors ki can be written: ∑ P (r) = c¯A + γi ηi e−iki ·r , (3) i


113

where ηi are the order parameters, which vary from 0 (disordered state) to 1 (fully ordered state). As for γi , it is the coefficient of symmetry for a given ordered phase i. Using this representation for P (r, t), the free energy functional in expression 2 can be expressed as: Fchem =

∑ 1 ∑ V (k)|P (k)|2 + kB T [P (r) ln(P (r)) + (1 − P (r)) ln(1 − P (r))] , 2N k r

(4)

where V (k) and P (k) are the Fourier transform of the mixing energy, W (r −r ′ ) and the density function, P (r)), respectively. The summation on k is performed over all points within the first Brillouin zone of a considered lattice. To illustrate this method, let us consider the kinetics of precipitation and microstructure development of the ordered L12 (γ ′ -Ni3 Al phase) particles in Ni-Al alloys. This ordered structure is a superstructure based on the fcc lattice. There are four translation variants of the L12 structure. This superstructure is generated by the three wave vectors k1 = 2π/aγ (100), k2 = 2π/aγ (010), and k3 = 2π/aγ (001). Here, k = (kx , ky , kz ), and kx , ky and kz are components of the vector k along the x , y and z axes parallel to the [100], [010] and [001] directions in the reciprocal space of fcc lattice. Besides, aγ is the lattice parameter of the fcc lattice. We assume that the function P (r, t)Al = P (r, t). Then, applying the static concentration wave formalism: ) η ( i2πx e + ei2πy + ei2πz . (5) 4 Substituting equation 5 into equation 4, the expression of the chemical free energy of the L12 becomes: P (r) = c¯Al +

1 3 ( η )2 kB T [ ( η) FL12 η) ( = c2 V (k = 0) + V (k) + ln c − 3 c− N 2( 2 4 4 4 4 η) ( η) +3 1−c+ ln 1 − c + (6) ( )4 ( ) 4 ( ) ( ) 3η 3η 3η 3η ] +3 c+ ln c + +3 1−c− ln 1 − c − , 4 4 4 4 where c is the solute composition (c = cAl ). Assuming the interactions within the 3 coordination shells, V (k) and V (k = 0) are defined as: { V (k = 0) = 12W1 + 6W2 + 24W3 (7) V (k) = −4W1 + 6W2 − 8W3 . Here W1 , W2 and W3 are the effective interchange interaction energies (mixing energies), for the first, second and third-nearest neighbours, respectively. In our simulations, the following set of parameters was used: W1 = −16.38 meV, W2 = −73.08 meV, and W3 = 20 meV. According to previous works [9, 10], a lattice mismatch between γ ′ ordered and γ disordered phases, which induces the elastic strains in Ni-Al alloys, engenders the morphological changes of ordered precipitates from a spherical to a cuboidal shape as the particles grow to a larger size. Then, to describe the morphological evolution of the microstructure in the Ni-Al compounds, the elastic energy contribution to free energy should also be taken into account. The elastic energy term Eelas can be calculated using the elastic theory of multiphase coherent solids proposed in [1], in the homogeneous modulus case approximation: ν0 ∑ Eelas = B(n)|c(k)|2 . (8) 2N k∈BZ 1

Here, the function B(n) describes the elastic properties and the transformation crystallography. n = k/k is a unit vector in the reciprocal space and c(k) is the Fourier transformation of the local concentration c(r). The summation is carried out on the first Brillouin zone (BZ1 ). Following [1], the function B(n) can be written:

114


Fig. 1: Microstructural evolution of the γ ′ precipitates (are shown in white) in disordered matrix (shown in black) in Ni-13 % at Al at T=1023 K. (a)-(d) simulation results , (e)-(h) TEM images obtained in [10]. (a) & (e): t = 15 min, (b) & (f): t = 2 h, (c) & (g): t = 4 h, and (d) & (h): t = 8 h.

B(n) = 3(C11 + 2C12 )ϵ20  2 C11 + 2C12 2  ( − ϵ0 C11 1+ξ 1+

 2ξ(n21 n22

n21 n23

n23 n22 )

3ξ 2 n21 n22 n23

+ + + 1+ ( ) C12 2 2 2 2 2 2 2 1 + 2 C12 + ) + ξ n + n n + n n (n 3 2 1 3 1 2 C11 C11

C44 C11

)

 , (9) n21 n22 n23

where ϵ0 = (aL12 − aγ )/(aγ [cγ − cγ ′ ]) is the concentration coefficient of crystal lattice parameter extension, cγ and cγ ′ are the solute concentrations in disordered and ordered phases, C11 , C12 and C44 are elastic constants of cubic fcc phase, n1 , n2 and n3 are the components of a unit vector n, and ξ = (C11 − C12 − 2C44 )/C44 is an elastic anisotropy factor. In our simulations, we assumed that the equilibrium compositions of γ and γ ′ phases at T = 1023 K were cγ = 0.115 and cγ ′ = 0.235.


115

This gives ϵ0 = 0.046. The values of the elastic constants were extracted from [11]: C11 = 209 GPa, C12 = 149 GPa and C44 = 96 GPa at 1000◦ K. Simulations were performed by numerically solving the equation 1, using the semi-implicit FourierSpectral method [12] in two dimensions, on a 1024x1024 simulation grid with periodic boundary conditions along x and y axes. We assume that L(r − r ′ ) does not vanish only for the nearest neighbor sites of the bcc lattice, r and r ′ , and that it is equal to the effective rate of the elementary diffusional jumps L1 between the nearest neighbor sites of the underlying fcc lattice. In the calculations, we used a reduced form of equation 1. The reduced time is t∗ = t/τ , where τ = (L1 |∆F |)−1 is the typical time for an elementary diffusion event and ∆F is the driving force of phase transformation. The energy is measured in units of ∆F . The constant time step ∆t∗ was chosen to be 0.1, and the mesh size was ∆x = 1.0. The initial state was a homogeneous solution with small composition fluctuations around the average composition. Spherical nuclei with an average radius greater than a nucleation radius were randomly introduced to the computation domain. Since our focus is only on the coarsening kinetics and not the sequence of transformations, the artificially chosen initial configuration does not affect our results. Approximately a thousand particles were formed during the nucleation stage. The temporal evolution of microstructure in Ni-13 % at Al alloy at 1023 K obtained in simulations and observed in experience is shown in figure 1, respectively. The coherent ordered γ ′ precipitates are shown in white in the γ ordered matrix is presented in black. In figure 1, the ordered particles are identified using the local long range order parameter. During the microstructure evolution, it could be observed that the anisotropic elastic energy induced cuboidal shape for γ ′ precipitates, and the anisotropic elastic interaction induced alignment along elastically soft [10] and [01] directions, even for the precipitates with the small size. Wang and Khachaturyan [3] were the first to suggest that antiphase relationships among nearest neighbours γ ′ precipitates were responsible for their resistance to coalescence. Figure 1 shows that our simulation results are in good agreement with the experimental data reported [10]. Using the ADF theory, we were able to reproduce not only the morphology of microstructure, but also to give the quantitative description of this kinetics. Although the discrete ADF method was successfully applied to describe the kinetics of orderdisorder phase transitions in Ni-Fe [13], Ni-V [14], Fe-Ga [15], Al-Li [16] and many other compounds, it cannot be applied to model phase transitions which are accompanied by the displacements of atoms. For example, the discrete model is unable to describe the atomic configuration near grain boundaries, solid-liquid transitions, self-assembly phenomena nor austenite-martensite transformations. Recently, the ADF continuous model was developed to model all these phenomena [7]. Continuous version of the ADF theory. The ADF theory on continuum space is a particular case of the ADF theory on Ising lattice. In this new approach, the distance between two atoms is bigger than the computation grid, and plays the role of Ising lattice in the previous ADF approach. Then, each atom occupies a certain number of simulation grid sites. We assume that in this case, each atom is a sphere made of finite elements, called atomic fragments or fratons. Such fragments are treated as pseudo-particles. In the multi-component systems, the number of types of fratons is equal to the number of different atomic components of the system. A small parameter determining the limit transition between the ADF theory on Ising lattice and its continuous version is the ratio of the underlying Ising lattice parameter to the atomic radius. So, if the atomic radius (interatomic interaction radius) is significantly greater than the Ising lattice parameter, the ADF theory correctly describes the atomic evolution, as the kinetics is decoupled from the underlying Ising lattice. Thereafter, we will call the continuous version of the ADF theory the QA (quasi-particle approach). Unlike the conventional approach describing the configuration of a multi-atomic system by coordinates of atomic centers, the proposed QA describes atomic configurations by occupation numbers of fratons. The proper choice of a model Hamiltonian describing the interaction of fratons should account

116


for both their "condensation" into atomic spheres and the movement of the spheres into the desirable equilibrium atomic configuration driven by the spontaneous minimization of the free energy. In the proposed fraton model, the configurational degrees of freedom are occupation numbers, s(r), where the vector r describes a site of the computational grid lattice. The function s(r) is equal to either 1 if the site r is occupied by a fraton, or 0 otherwise. At a finite temperature T , the system is described by averaging over the time-dependent ensemble. The averaging gives the occupation probability, ρ(r, t) =< s(r) >t , 0 ≤ s(r) ≤ 1, where the symbol < · > designates the average over the time-dependent ensemble. The occupation probability, ρ(r, t), in fact, is the probability that a point r is anywhere inside the volume of any sphere describing an atom at time t. In the QA, the Helmholtz energy F is determined by the equation 2, where now F is a functional of the occupation probabilities {ρ(r, t)}. The function w(r − r ′ ) can be interpreted as a fraton-fraton interaction. The chosen fraton-fraton interaction model potential w(r − r ′ ) should guarantee the selfassembly of initially randomly distributed fratons into a desirable structure. The Fourier transform (FT) of such a potential is: w eαβ (k) =

1 ∑ wαβ (r)e−ik·r , N0 r

(10)

where the summation is carried out over all sites of the Ising lattice (grid). The wave vector, k, is defined at all quasi-continuum points, k, of the first Brillouin zone of the computational grid, that is, at all the N0 points in the k-space permitted by the periodical boundary conditions. To model the spontaneous self-assembly of a disordered atomic distribution into the desired 3D atomic structures of arbitrary geometric complexity, we introduced a new approach in a formulation of model potential, based on the concepts of structural clusters and cluster amplitudes. A structural cluster is defined as a minimum size group of geometric points (cluster), whose size is just sufficient to fully reproduce the topological features of the desirable final configuration of the α-atoms. An amplitude of structural cluster of the kind α is defined as: Ψclst α (k) =

∑

ωj (α, k)e−ik·rjα ,

(11)

jα

where summation is carried out over all points of the α-kind cluster. The index ja numbers these points, and ωj (α, k) is the weight of the contribution of each point of the cluster to its amplitude. The constants ωj (α, k) are chosen to reproduce the desirable thermodynamic and mechanical properties of the simulated atomic aggregate. With these definitions, we present the Fourier Transform of the model potential as the sum of what we call the short-range and long-range interactions: clst ∗ Wαβ (k) = λ1 θα (k)δαβ + λ2 (k)Ψclst α (k)Ψβ (k) .

(12)

Here, δαβ is the Kronecker delta function. The first term in equation 12 describes the short-range fratonfraton pair interaction. It is responsible for the spontaneous "condensation" of fratons into atomic spheres. The simplest choice of the short range part of the interaction potential is a step function. The function θ(r) and its Fourier Transform θˆα (k) are the spherically symmetric functions schematically presented in figure 2, where the constant r1 is the radius of attractive part of short range interaction potential. r1 sets the atomic radius. ∆r is the width of the repulsion part, and ξ = |θmax (r)|/|θmin (r)| is the ratio between the modules of minimum and maximum value of the shape function presented in figure 2(a). In equation 12, λ1 is a constant which determines the strength of the short-range fratonfraton interaction. The second term of equation 12 describes the long-range part of the fraton-fraton pair interaction, which is responsible for the mutual arrangement of atoms in the final desired configuration. It is presented as a bilinear expansion in cluster amplitudes, Ψclust α (k). λ2 in equation 12 is a fitting


117

parameter determining the strength of the long-range interaction. The index α, can be dropped in a single-component system. Then equation 12 is simplified to: 2 W (k) = λ1 θ(k) + λ2 (k) Ψclst (k) .

(13)

Since the FT of fraton-fraton interaction is formulated for a specifically oriented cluster, this orientation automatically lifts the angular isotropy of the system.

Fig. 2: Left: Schematic representation of the short-range potential θ(r). Right: Fourier transform of θ(r), using the following input parameters: ξ = 4, ∆ˆ r = 0.25. The temporal evolution of the density function {ρ(r, t)} is described by the microscopic kinetic equation in equation 1. Solving this equation, self-assembly phenomena in different systems can be modelled. Here, we used reduced parameters, and in particular, the average density, defined as ˆ ρ¯α = ρatα 4πRα3 /3, where ρˆ¯α = ρatα 4πRα3 /3. ρatα = Nα /V is the density of α atoms in the ground state, Nα is the number of the atoms of sort α, V is the total volume of the system, and Rα is the atomic radius of this atom. According to this definition, the reduced density ρˆ¯α is also a fraction of all computational grid sites occupied by fratons of the kind α. The input parameter ξ of the energy α(r) is measured in units of kB T0 , where T0 is the solidification temperature. The lengths are measured in units of r1 , which is very close to the atomic radius. The grid lattice increment ˆl, (the spacing of the underlying Ising lattice), is defined as a fraction of the atomic radius. The temperature Tˆ is also measured in units of T0 . The reduced time tˆ is measured in units of typical atomic migration time τ0 . The reduced kinetic ˆ coefficients L(r) are measured in units of τ0−1 . To illustrate the versatility and efficiency of the QA, we tested its application to the modelling of the self-assembly of three groups of 3D structures of increasing complexity, namely the non-ideal gas/liquid, single-component crystals, two-component crystals, and polymers with a helix and a double helix structure mimicking biological macromolecules [7]. For example, we will consider a selfassembly of atoms to diamond structure. The diamond lattice has the fcc Bravais lattice with a two-atom basis, the atoms of the basis being displaced by the vector, a(1/4, 1/4, 1/4), where a is the crystal lattice parameter. The coordinates are given as fractions of the crystal lattice parameter along the cube sides. A chosen structural cluster in this case is a cubic unit cell of the diamond lattice. It consists in eight points in the positions of the atoms in this unit cell. The four points forming a fcc lattice cubic cell (0, 0, 0), a(1/2, 1/2, 0), a(0, 1/2, 1/2), a(1/2, 0, 1/2), and additional four points obtained from them forgoing by the basis shift a(1/4, 1/4, 1/4). The function Ψclust (k) was constructed by using its definition in equation 11 and the coordinates of the chosen structural cluster points:

118


( )( ) a a a a Ψclst (k) = 1 + e−i 4 (kx +ky +kz ) 1 + e−i 4 (kx +ky ) + e−i 4 (ky +kz ) + e−i 4 (kx +kz ) .

(14)

Here, ki = 2πmi /(aN ), where i = x, y or z, mi = 1, N , N is the number of simulation grid in a given direction. Using this function and the short range interaction θ(r), the Fourier transform w(k) e ˆ ˆ was defined. In this simulation, the next set of parameters was used: λ1 = 14.085, λ2 = −7.042, ˆ = 1, ρˆ¯ = 0.07, ˆl = 0.286, ξ = 2, ∆ˆ D r = 0.17, a = 4.57, and Tˆ = 0.732. The spontaneous self-organization of fratons into the diamond structure is shown in figure 3. The initial configuration was a single unit cell embryo of the diamond embedded into the one-component gas consisting of disordered fratons. A very interesting aspect of this self-assembly of the diamond crystal is the occurrence of the transient cubic body centered (bcc) structure at the early stages of evolution (tˆ = 1000, step (c) in figure 3). The lattice parameter of this bcc structure is half that of the diamond structure. This transient state gradually transformed to the diamond structure by a is followed by the gradual disappearance of some atoms in the bcc structure and the formation of the diamond structure.

Fig. 3: Example of a self-assembly of fratons into diamond structure at reduced times (a) tˆ = 0, (b) tˆ = 600, (c) tˆ = 1000, (d) tˆ = 2800, and (e) tˆ = 3000. The initial configuration is the atomic cluster of a diamond structure placed in the centre of the simulation box. The size of the simulation box is 64 × 64 × 64. The QA opens a way to answer numerous outstanding questions concerning the atomistic mechanisms of the formation of defects (dislocations, grain boundaries, etc.), nucleation in solid-solid transformations, the formation of polymers due to aggregation of monomers in their solution, folding and crystallisation of polymers, and their responses to external stimuli. This list can be signicantly extended. Especially interesting are the modelling results describing the spontaneous self-assembly of monomers into a single-stranded polymeric helix and the formation of a double-helix structure obtained by aggregation of complementary monomers on the single-stranded helix playing the role of a template. This result may be also considered as an attempt to formulate and execute the simplest prototyping of the spontaneous formation of homopolymeric DNA from a liquid solution of monomers playing the role of nucleotides.


119

Phase Field Method The ADF theory can be seen as an atomic scale phase field. However in many cases, especially for industry application, it is important to predict and to understand the microstructural evolution at mesoscale. To achieve this goal, the standard phase field model was extensively used for microstructure simulation studies. In this approach, each phase or domain in a microstructure is characterized by a set of field variables (compositions and/or order parameters). Within the phase field model [8], the microstructure can be characterized by the coarse grained concentrations cα (r, t) of components α, and by the long range order (LRO) parameters ηi characterising the crystal symmetry of coexisting phases (i = 1, n, where n is equal to the number of coexisting phases in considered system). In general, in the phase field method, the microstructural evolution can be obtained by solving the time dependent Ginzburg-Landau equation for the LRO parameters: ∂ηi (r, t) ∂F = −L + ζi , ∂t ∂ηi (r, t)

(15)

and the non-linear Cahn-Hillard diffusion equation for the concentration field: [ ( )] ∂F ∂cα = ∇ · M (r)∇ + ξα . ∂t ∂cα (r, t)

(16)

Here L and M are kinetic coefficients, ζi and ξα are the Langevin's noise terms describing the compositional and structural thermal fluctuations, respectively. F is the total free energy. In general, the total free energy contains different contributions such as elastic, magnetic or electric. An application of this method is provided by the disrupted phase separation in AgCu system under ion irradiation. To illustrate this approach, the separation kinetics in the AgCu system under irradiation will be considered. The phase diagram of the AgCu compound predicts that at low temperature, this system undergoes a phase separation, which leads to the formation of Ag-rich and Ag-poor phases (see figure 4). 1357.8 1200

T (K)

1052.25

0.601

1234.9

1000

800

liquid/solid binodal spinodal

600 0

0.1 0.2 0.3 0.39 0.5 0.6 0.7 0.8 0.9

1

Ag average concentration c¯

Fig. 4: Equilibrium phase diagram for the AgCu alloy. However, it was shown that under irradiation, a patterned structure may be observed [17], resulting from a ballistic disordering induced by the slowing down of impinging ions [18]. The manufacturing of such steady state patterned microstructures, using the ion beam mixing process [19], has sparkled interest from the microelectronic industry notably [20]. Indeed, this fabrication process has the advantage to ensure the control of the microstructure features by such parameters as irradiation temperature and ion flux.

120


To understand the reason of this pattern formation, and to predict the spatial periodicity of this patterning, the phase field model can be used. In this case, both coexisting phases have the same crystal structure, and the kinetics of phase separation can thus be simulated by the sole Cahn-Hillard equation (equation 16). In binary alloys, only one concentration field variable c is needed to describe the morphology of microstructure. Here c(r, t) designates the silver concentration cAg (r, t), and its evolution under irradiation is governed by equation 16. The free energy of the heterogeneous system F (c) can be written as: ∫ [ F (c) =

] κ fh (c) + |∇c|2 dΩ. 2

(17)

Ω

fh (c) is the local free energy density, which can be written under the Landau polynomial form: fh (c) =

a2 (T ) a3 (T ) a4 (T ) (c − c¯)2 + (c − c¯)3 + (c − c¯)4 , 2 3 4

(18)

where c − c¯ is the deviation of concentration from the average concentration of the system. In equation 18, a cubic term is added to fh , in order to account for first order phase transitions [21]. fh displays ¯ two minima corresponding to two stable equilibrium states. In equation 16, the atomic mobility M plays a dominant part in the definition of the time scale of the dynamics. An average expression of ¯ , derived by Martin in [22] is usually used: M ¯ = c¯(1 − c¯) D. ˜ M kB T

(19)

˜ theoretically ˜ is the chemical interdiffusion coefficient [23]. D Here, c¯ is the mean concentration, and D depends on the concentration. Yet, in the case of similar isotopic diffusion coefficients of both species in the alloy [24], it can be approximated by an average isotopic diffusion coefficient, thus independent of c [25]. Coefficients a2 , a3 and a4 , defining the free energy density fh (c) in equation 18, can be estimated using an experimental AgCu phase diagram [26]. The equilibrium interdiffusion concentration profile of silver, and the stiffness parameter κ can be determined using a mixed molecular dynamics (MD)-Monte Carlo (MC) approach. The phase field interface profile intrinsically matches the physical interface given by MC simulations. This correspondence is illustrated in figure 5. From an average over every MC simulation, κ can be considered constant: κ = 0.15 (±0.005) eVÅ−1 . 1 Silver concentration c

Monte Carlo Fitted profile

0.5

0 70

80

90 ˚ Depth z (A)

100

Fig. 5: Silver concentration profile at the interface, determined by our Monte Carlo simulation (circle) for T = 700 K, and corresponding phase field profile (line) based on a fit by a Jacobi elliptic sine.


121

The versatility of the phase approach allows to add irradiations effects quite easily. In a first approximation, the impact of irradiation on the microstructure of an alloy results from two main mechanisms: ballistic effects induced by atomic relocations, and diffusion enhancement due to point defects recreation. Martin suggested that ballistic effects could be described by an athermal term [27]. The local fluctuations of concentration induced by irradiation can then be expressed by a local matter balance [28]: ∂c = Γ × (pR ∗ c − c). (20) ∂t Here, Γ is an atomic relocation frequency, and pR is a probability distribution of displacements in subcascades. It involves Sigmund and Gras Marti ion mixing formalism [29]. Ion mixing asserts that ballistic effects can be treated as an homogeneous source term. The evolution of the microstructure under irradiation is thus the result of two mechanisms, acting separately on two different time scales. This makes the microstructure only sensible to the average effect of numerous cascades covering one another. The strong heterogeneity due to subcascades [30] can hence be passed over. Theoretical an numerical studies [31] suggest that pR followed an exponential decay, mimicking small displacements in subcascades: ( ) 1 −r pR (r) = exp . (21) 2 2πR R In that case, R is the mean relocation distance. As for Γ, it is the product of the irradiation flux Φ, and a relocation macroscopic cross section σ d : Nd Vcell × , (22) Ncell Lc where Vcell and Ncell are the cell volume and number of atoms per cell. Nd is the number of relocated atoms per cascade. The second phenomenon related to irradiation is the enhancement of atomic diffusion, due to point defects recreation. Introducing additional vacancies and interstitials concentration irr irr δcirr to the mobility can be written [32]: V and δcI , the irradiation contribution M Γ = σ d Φ,

σd =

] c¯(1 − c¯) [ irr δcV DV + δcirr (23) I DI , kB T where DV and DI are the vacancy and interstitial diffusion coefficients in the alloy. This expression irr requires that DV and DI are very close for each species of the alloy [33]. The values of δcirr V and δcI can be estimated by the chemical model of Sizmann [32], taking thermal recombination, and sink absorption at interfaces into account [34]. As point defects migration takes place on a very small time scale (∼ 10−9 s) in comparison with precipitates growth (∼ 10−6 s), their concentration can be set at their equilibrium value. The total mobility M then reads: M irr =

¯ + M irr . M =M

(24)

Irradiation parameters can be computed using the Binary Collision Approximation [35]. Using for instance the BCA code MARLOWE [36] on polycrystalline AgCu under 1 MeV Krypton ions, the simulations provides with the relocation cross section σ d = 2.01 × 10−13 cm2 , the mean relocation range R = 3.04 Å, and the number of Frenkel pairs per cascade before thermal annihilation and sink absorption, NV = 19500 [37]. For an experimentally reachable incident flux Φ = 1.21×1013 cm−2 s−1 [38], the ballistic relocation frequency Γ = σ d Φ = 2.53 s−1 , and the total mobility M can be computed [25]. Finally, the phase field model for irradiated binary alloys relies upon the modified Cahn-Hilliard equation: [ ] ∂c (25) = M ∇2 fh′ (c) − κ∇2 c + Γ(pR ∗ c − c), ∂t

122


φ

φ

φ

where the ballistic term in equation 20 was simply added to the thermodynamic contribution in equation 16. Based on phase field simulations displayed in figure 6, three evolutions of the microstructure can be distinguished, depending on the irradiation flux Φ, and the temperature T . This is due to Φ and T stepping in the mobility M (Φ, T ), the free energy potential coefficients a2,3,4 (T ), and the relocation frequency Γ(Φ).

y

y x

(a) Phase separation Γ ≪ 1/t0

y x

(b) Patterning Γ ∼ 1/t0

x

(c) Solid solution Γ ≫ 1/t0

Fig. 6: Reduced concentration ϕ from our phase field simulations: aleatory initial condition with ϕ¯ = 0 (c(t = 0) = c¯) and λ = 0.25 (T < 300 K), on a 1000 × 1000 mesh grid, with ∆t = 0.1 and tfinal = 106 . Top left: demixtion (W, R′ ) = (0.005, 2.4). Top right: patterning (W, R′ ) = (0.12, 2.4). Bottom: solid solution (W, R′ ) = (1, 2.4). The microstructures displayed in figures 6(a), 6(b) and 6(c) define three microstructural "pseudophases" in the temperature-flux plane shown in figure 7. For high fluxed and low temperatures (figure 6(c), domain (1) in figure 7), we observe a homogeneous system corresponding to the formation of a solid solution. In this case, the disordering due to long range displacements induced by irradiation overcomes the expected phase separation driven by thermal diffusion: Γ ≫ 1/t0 . This corresponds to precipitates radius close to 0 at all time. On the contrary, for high temperature and moderate fluxes (figure 6(a), domain (3) in figure 7), diffusion ordering effects are predominant as Γ ≪ 1/t0 . The system undergoes a macroscopic phase separation, similar to Ostwald ripening [39]. This corresponds to an unbounded growth of precipitates radius. Between these two domains lies the patterning domain (2) (figure 6(b)). In this domain, precipitates reach a finite size and remain stationary. This is due to ballistic and diffusional effects coming to equilibrate one another. It is associated to a mean precipitates radius, quickly reaching its finite asymptotic value. It should be noted that contrary to chemical disordering, ballistic disordering is of finite range R, and patterning can only be achieved for a sufficient value of R. Indeed, it allows ballistic relocations to display longer range than diffusion migrations. This numerical diagram is consistent with diffraction experiments on the silver-copper system [38]. Summary In summary, we have presented the atomistic and mesoscopic approaches to model phase transformations in solids. It was shown that the ADF theory can be used to model the isostructural phase transformations governed by the different driving forces: chemical, elastic, magnetic etc. This method is based on the resolution of microscopic diffusion equation on the Ising lattice, where each site is occupied by atoms. For instance, the evolution of the morphology of the microstructure in Ni-Al alloys was considered. It was shown that the elastic contribution in the free energy determines the shape and alignment of the ordered Ni3 Al particles. Then, the continuous version of the ADF theory was introduced. In this approach, the simulation grid spacing is smaller that the interatomic distance. Then, each atom occupies a certain number of simulation grid points. The grid points which are occupied by one fraction of atom are called fraton. The short and long range interactions between fratons first induce the coagulation of fratons in the atomic spheres, and second the periodic alignment of these spheres. Using this approach, the self-assembly of a random distribution of atoms into various complex atomic configurations of increasing complexity were successfully modelled [7]. The structures

Irradiation flux Φ × 10−10 (cm−2 .s−1 )


104

123

Solid solution (1)

103 Patterning (2)

Macroscopic phase separation (3)

102

Solid solution (exp.) Patterning (exp.) Demixtion (exp.) 101 250

300 320

350

400

450

500

Temperature T (K)

Fig. 7: Phase diagram in the temperature-flux plane for the silver-copper system subjected to Krypton ions (1 MeV) irradiation, displaying three microstructure behaviours: solid solution (top left), patterning (top center) and phase separation (right). For low temperatures and moderate flux, the system displays stationary precipitates (disrupted coarsening). Circle, triangle and cross: diffraction results from Wei and al. [38]. ranged from single-component and multi-component crystals, to single- and double-helix polymers. In this paper, we selected one example of self assembly, namely the growth of diamond structure from completely disordered distribution of atoms. In principle, this approach can even be used to describe the 3D pattern formation by any macroscopic objects and optimization of their properties. In this case, the "fratons" being fragments of these objects are also macroscopic. We believe this approach can be used to describe a wide manifold of phenomena related to self-assembly. In the last part of this paper, the standard phase field method was introduced. This phenomenological approach can be used to model the microstructural evolution at mesoscale. The advantage of this method is that many different phenomena related to diffusion of atoms and phase changing can be easily included in the phase field equations. For example, the kinetics of phase separation in Ag-Cu alloys under irradiation was considered. It was shown that the late stage microstructure in these compounds depends on the temperature and ion flux. At low temperature and high fluxes, the disordered phase is observed. Then, at higher temperature and moderated fluxes, a patterned microstructure appears. Finally, for high temperatures and low fluxes the phase separation is recovered. All the methods we presented in this paper, have been the main approaches for the last twenty years to model the microstructure evolution and pattering kinetics in different types of materials. However, to construct a true multi-scale model, the link between the phenomenological parameters used in these approaches and physical quantities should still be established. Acknowledgment The research was funded from the European Community's Seventh Framework Programme (FP7PEOPLE-2013-IRSES) under EC-GA no. 612552.

124


References [1] Long-Qing Chen and AG Khachaturyan. Computer simulation of structural transformations during precipitation of an ordered intermetallic phase. Acta metallurgica et materialia, 39(11):25332551, 1991. [2] Armen G Khachaturyan. Theory of structural transformations in solids. Courier Corporation, 2013. [3] Y Wang, L-Q Chen, and AG Khachaturyan. Kinetics of strain-induced morphological transformation in cubic alloys with a miscibility gap. Acta Metallurgica et Materialia, 41(1):279-296, 1993. [4] R Poduri and L-Q Chen. Computer simulation of atomic ordering and compositional clustering in the pseudobinary ni 3 al-ni 3 v system. Acta materialia, 46(5):1719-1729, 1998. [5] LQ Chen and AG Khachaturyan. Dynamics of simultaneous ordering and phase separation and effect of long-range coulomb interactions. Physical review letters, 70(10):1477, 1993. [6] Yongmei M Jin and Armen G Khachaturyan. Atomic density function theory and modeling of microstructure evolution at the atomic scale. Journal of applied physics, 100(1):013519, 2006. [7] Mykola Lavrskyi, Helena Zapolsky, and Armen G Khachaturyan. Quasiparticle approach to diffusional atomic scale self-assembly of complex structures: from disorder to complex crystals and double-helix polymers. npj Computational Materials, 2:15013, 2016. [8] Long-Qing Chen. Phase-field models for microstructure evolution. Annual review of materials research, 32(1):113-140, 2002. [9] Y Wang, D Banerjee, CC Su, and AG Khachaturyan. Field kinetic model and computer simulation of precipitation of l1 2 ordered intermetallics from fcc solid solution. Acta materialia, 46(9):2983-3001, 1998. [10] Yong Ma and Alan J Ardell. Coarsening of γ (ni-al solid solution) precipitates in a γ ′ (ni 3 al) matrix. Acta materialia, 55(13):4419-4427, 2007. [11] Heinrich Pottebohm, Günter Neite, and Eckhard Nembach. Elastic properties (the stiffness constants, the shear modulus and the dislocation line energy and tension) of ni al solid solutions and of the nimonic alloy pe16. Materials Science and Engineering, 60(3):189-194, 1983. [12] LQ Chen and Jie Shen. Applications of semi-implicit fourier-spectral method to phase field equations. Computer Physics Communications, 108(2):147-158, 1998. [13] IV Vernyhora, HM Zapolsky, R Patte, and D Ledue. Atomic density function modeling of microstructure evolution in ni 3- x fe x alloys. Journal of Magnetism and Magnetic Materials, 351:52-59, 2014. [14] Helena Zapolsky, Sebastien Ferry, Xavier Sauvage, Didier Blavette, and Long-Qing Chen. Kinetics of cubic-to-tetragonal transformation in ni-v-x alloys. Philosophical Magazine, 90(14):337-355, 2010. [15] J Boisse, H Zapolsky, and AG Khachaturyan. Atomic-scale modeling of nanostructure formation in fe-ga alloys with giant magnetostriction: Cascade ordering and decomposition. Acta Materialia, 59(7):2656-2668, 2011.


125

[16] R Poduri and L-Q Chen. Computer simulation of morphological evolution and coarsening kinetics of δ ′ (al 3 li) precipitates in al-li alloys. Acta materialia, 46(11):3915-3928, 1998. [17] R. Enrique and P. Bellon. Compositional patterning in systems driven by competing dynamics of different length scale. Phys. Rev. Lett., 84(13):2885-8, March 2000. [18] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium. Rev. Mod. Phys., 65:851-1112, Jul 1993. [19] H. Bernas, J.-Ph. Attané, K.-H. Heinig, D. Halley, D. Ravelosona, A. Marty, P. Auric, C. Chappert, and Y. Samson. Ordering intermetallic alloys by ion irradiation: A way to tailor magnetic media. Phys. Rev. Lett., 91:077203, Aug 2003. [20] Y. Cheng. Mater. Sci. Rep., 45, 1990. [21] P. Tolédano and V. Dmitriev. Reconstructive phase transitions: in crystals and quasicrystals. World Scientific, 1996. [22] G. Martin. Atomic mobility in Cahn's diffusion model. Phys. Rev. B, 41(4):2279-2283, 1990. [23] A. R. Allnatt and A. B. Lidiard. Atomic transport in solids. Cambridge University Press, 2003. [24] D. B. Butrymowicz, J. R. Manning, and M. E. Read. Diffusion in copper and copper alloys part iv. diffusion in systems involving elements of group viii. J. Phys. Chem. Ref. Data, (1), 1976. [25] G Demange, L Lunéville, V Pontikis, and D Simeone. Prediction of irradiation induced microstructures using a multiscale method coupling atomistic and phase field modeling: Application to the agcu model alloy. Journal of Applied Physics, 121(12):125108, 2017. [26] J. L. Murray. Calculations of Stable and Metastable Equilibrium Diagrams of the Ag-Cu and Cd-Zn Systems. l(February):261-268, 1984. [27] G. Martin. Phase stability under irradiation: Ballistic effects. Phys. Rev. B, 30(3):53, August 1984. [28] D. Simeone, L. Lunéville, and G. Baldinozzi. Cascade fragmentation under ion beam irradiation: A fractal approach. Phys. Rev. E, 82, July 2010. [29] P. Sigmund and A. Gras-Marti. Theoretical aspects of atomic mixing by ion beams. Nucl. Instrum. Methods, 182–183, Part 1(0), 1981. [30] E. Antoshchenkova, L. Luneville, D. Simeone, R.E. Stoller, and M. Hayoun. Fragmentation of displacement cascades into subcascades: A molecular dynamics study. J. Nucl. Mater., 458:168 - 175, 2015. [31] R. A. Enrique, K. Nordlund, R. S. Averback, and P. Bellon. Simulations of dynamical stabilization of Ag–Cu nanocomposites by ion-beam processing. J. Appl. Phys., 93(5):2917, 2003. [32] R. Sizmann. The effect of radiation upon diffusion in metals. J. Nucl. Mater., 69:386-412, 1978. [33] G. Demange. Mise en œuvre d'une approche multi-échelles fondée sur le champ de phase pour caractériser la microstructure des matériaux irradiés : application à l'alliage AgCu. PhD thesis, École Centrale de Paris, 2015. [34] A. Barbu and G. Martin. Materials under irradiation.

126


[35] M. T. Robinson and I. M. Torrens. Computer simulation of atomic-displacement cascades in solids in the binary collision approximation. Phys. Rev. B, 9(12), 1974. [36] M. T. Robinson. The binary collision approximation: background and introduction. Proceedings of the International Conference on Computer Simulations of Radiation Effects in Solids, 141, August 1992. [37] G. Demange, E. Antoshchenkova, M. Hayoun, L. Lunéville, and D. Simeone. Simulating the ballistic effects of ion irradiation in the binary collision approximation: A first step toward the ion mixing framework. J. Nucl. Mat., 486:26 - 33, 2017. [38] L. C. Wei and R. S. Averback. Phase evolution during ion-beam mixing of Ag–Cu. J. Appl. Phys., 81(2):613, 1997. [39] I. M. Lifshitz. The kinetics of precipitation from supersaturated solid solutions. J. Phys. Chem. Solids, 19(1-2):35-50, April 1961.