KMC Simulation of 2D Immiscible Binary Alloys Driven by
Irradiation
Yaofeng Chen and Mohsen Dadfarnia
Abstract
KMC simulation is performed to investigate the ion irradiation effect on a simple 2D AB alloy, which undergoes phase separation at low temperature in its thermal equilibrium phase diagram. At high ballistic jump frequency, ion irradiation can produce a solid solution-liked dynamic steady state. At intermediate ballistic jump frequency, a nanoscale compositional patterning steady state is found to exist in the 2D AB alloy under ion irradiation, which consists with the finding of the 3D KMC simulation in the literature.
Introduction
The phase stability under irradiation has received much attention motivated by technological problems such as the long term integrity of nuclear waste confinement materials and alloy preparation by ion implantation or ion beam mixing ý[1]ý[2]. Irradiation may derive a system away from thermodynamic equilibrium by local mixing of the components. For a solid under irradiation, colliding particles can produce the atomic rearrangements. The magnitude of the atomic displacements vary according to the energy transferred by the impinging particle. Above a threshold energy a sequence of collisions or ballistic displacement is originated.
An interesting physical situation arises when irradiation induced mixing acts in opposition to the thermodynamically driven kinetics (e.g. ý[3] and ý[4]). A mixture of immiscible alloys without external force undergoes a phase separation. This competition between mixing by irradiation and separation by thermodynamic kinetics can be led to different steady-state configurations under the continuous driving forcing.
There are two general approaches to the simulations: deterministic and non-deterministic. Molecular dynamics (MD) simulations is a type of deterministic simulation. In Molecular Dynamics (MD) simulations, one integrates the equation of motion for individual atoms in a many body system and therefore, those simulations follow the real dynamics of the system. However, time steps should be taken very small to capture even high frequency dynamics of individual atoms. Thus, MD simulations are restricted to small systems simulation in small time span.
Direct Monte Carlo (MC) simulation is a quasi-random or non-deterministic approach. According to a stochastic algorithm, MC method evolves to explore the phase space of a system. The MC method follows a Markov process to evolve a system to an equilibrium state. This method can be used to obtain static or equilibrium properties of model systems.
Utilizing a continuous model Enrique and Bellon ý[4] showed analytically that under certain conditions (see Fig. 1) it is possible to get a steady state patterning for a immiscible binary alloy under irradiation. They ý[7] also performed a three dimensional KMC simulation to confirm the prediction of their analytical model (Fig. 2).
Figure 1. Steady state regimes as a function of relocation range for ballistic
mixing R and frequency of exchange per mobility γ. A and C are material
constants ý[4] . 
Figure 2. Temporal evolution of the system starting from different initial
conditions: (a) random solid solution and (b) pure A/B bilayer,
at an intermediate frequency of atomic relocations Г=200 s-1. with Relocation distance R=3ann. The picture shows (111)
cuts of the 3D KMC simulation cell ý[7].
In terms of experiment, Wei and Averback ý[8] investigated ion-beam mixing of the immiscible Ag-Cu alloy system using 1.0 MeV Kr ion irradiation at different temperatures. They observed that below room temperature the Ag-Cu thin film under radiation was mixed completely with the formation of single phase homogeneous alloy and at moderate temperature they were able to obtain a solubility enhancement. At intermediate temperature two or three phases are formed.
Enrique et al. ý[9]ý[10] performed a systematic study of ion beam mixing in Ag-Cu multilayered thin films irradiated from room temperature to 225°C, based on the examination of cross-sectional samples by electron microscopy. They found that independent of irradiation temperature extensive grain growth and breakage of the layered structure takes place. At intermediate temperatures, they observed that competition between thermal decomposition and irradiation mixing resulted in a nanometer scale phase separation (see Fig. 3).
Figure 3. The top row shows bright-field (BF), diffraction pattern (DF) and
dark-field (DP) images obtained at the TEM. The bottom row shows a
compositional analysis along a line shown in a micrograph, performed at STEM.
(a) The as-deposited sample consists of ten period of (6.7 nm Cu/ 11.2 nm Ag)
layer on a (111) Si substrate. (b) Sample irradiated at room temperature to 
ý[9]ý[10].
In this project, we study a two dimensional immiscible binary alloy of Ag-Cu which continuously undergoes ballistic atomic exchange due to irradiation. KMC method is utilized to simulate the system and study the compositional patterning.
Kinetic Monte Carlo
The main idea behind Kinetic
Monte Carlo is to use transition rates based on a dynamic model, with time
increments formulated so that they relate to the microscopic kinetics of the
system. The
where 
and 
are successive states
of the system, 
is the probability
that the system is in state at time t, and 
is the transition
probability per unit time that the system will undergo a transition from state to state .
The Monte Carlo method may
be utilized to simulate a Poisson process and therefore an exact correspondence
between 
is
given as
where 
is a random variable
counting the number of events 
which have occurred
within a time 
. From Eq. the expected number of events occurring within a
time is 
. The probability density of time between successive events in Poisson process is
.From this probability density, the mean time
between successive events is 
. Then if a system under study can be cast into a set of
independent events with known rates, then the time between events has an
exponential distribution given in Eq. that can be sampled to
determine the time that it takes for an event to occur. This allows us to get
the real time 
between successive
events as
.
where 
is the uniform random
number distributed between 
. This random sampling of the Poisson time distribution for
each chosen event ensures that a direct relation between
Fichthorn
and Weinberg ý[5] showed that if three criteria are met, the
- Transition probabilities should reflect a so-called dynamical hierarchy in addition to obeying detail balance. This means they should be ordered and divided by the maximum of the transition probabilities such that the maximum be unity,
- Time increments upon successful events are formulated correctly in terms of the microscopic kinetics of the system,
Monte Carlo Algorithm
Our approach follows the same methodology provided by Enrique and Bellon ý[3] for irradiation and Abinandanan et al. ý[11] for diffusional phase transformation.
1.
Alloy model with ion irradiation
(a) Alloy model
In our simulation,
we use an A50B50 alloy model with one vacancy. The
details of the model are described in the paper by Abinandanan
et al. ý[11]. A rigid 2D square lattice NxN (in our simulation N =64 or 128) is
considered whose sites are occupied by atoms of Ag or Cu or by vacancy
designated as V. A single vacancy is introduced in the lattice for the
simplicity. The atoms and vacancies interact with each other through bonds of energy 
, 
, 
, 
, and 
. Only the nearest neighbor interactions are considered. The
energy of the system at any configuration, i,
is given by:
where 
is the label of each site, 
is the label of
nearest neighbor (n.n.) site, and X and Y can be Ag,
Cu or V. 
is an occupation
operator: 
equals unity if site is occupied by species
X in configuration i , and equals zero
otherwise.
For the thermal diffusion the transition probabilities are given by:
where 
is represents the
attempt frequency for vacancy jump, 
is the saddle point
energy, 
is the energy of the configuration. The saddle point energy
is independent of the direction of the transition. (i.e. 
). The exponent in Eq. is simply a sum of a constant saddle-point energy plus the energy required to
break the bonds linking the vacancy and the exchanging atom to their
environment. The constant saddle-point term can be fixed by setting a value for
the vacancy migration energy at a given concentration.
For simplicity, our 2D simulation use the
same materials parameters in ý[3] and ý[12] that actually perform the 3D KMC simulation, i.e., 
, 
, 
, 
and 
.
(b) Modeling the ion irradiation
To model the events of ion
irradiation on the AB alloy, we use the similar method that are described by Enrique
and Bellon ý[3], especially in the reference ý[12]. Here we only consider the effect of the heavy ions,
like Ar and Kr. The light ions like He have a
different behavior ý[12].
To model the effect of the irradiation, they introduce forced atomic
relocations at a controlled frequency per atom 
. For simplicity, instead of modeling the cascade created by
every ion irradiation event, they perform only one forced exchange at one time
when one ion hit on the alloy. In fact, for every irradiation event, one atom
and a jump distance are randomly selected, in compliance with the distribution
of atomic relocations obtained by their MD simulation. Based on their MD simulation,
the distribution can be fitted well into an exponential decay function, i.e.,
0.07Å
and R=3.08 Å. R is approximately the nearest neighbor
distance of the alloy ý[12] (see Fig. 4).
Figure 4. Normalized sistribution of
atomic relocations. Ne, Ar
and Kr give essentially the same normalized distribution. He, on the other
hand, prodices a relocation distribution that weighs
short-range jumps more heavily. Also shown are fittings to the data based on
experimental decays:
0.07 Å
Å) for the case of Ne, Ar, and Kr, and 0.08 Å
Å) for the case of He ý[12].
In our 2D model, we will set R=ann, the first nearest
neighbor distance of our 2D square lattice. Noticing KMC simulation is
performed on a rigid lattice, one need to lump the data from MD into different
neighboring shells.
When a ballistic jump will occur in the simulation, first we need to randomly choose a site in the lattice. And then a jump distance with the distribution of the exponential decay is randomly chosen and then a random direction (a random angle) is assigned to this distance. The closest lattice point to that location is set as the final destination of the relocation atom ý[12]. The transition probabilities for the ballistic mixing caused by irradiation are
where is the frequency for
ballistic exchange per atom.
2. Kinetic
In order to have an efficient simulation, we
used the resident time algorithm ý[3] ý[11] or called n-fold way method. This algorithm can
ensure each MC step have a successful transition from one configuration to
another configuration. The total probability 
is the sum of the
probability of thermal jump of the only vacancy to its four nearest neighbors, i.e., 
, and the probability of the ballistic jump caused by the
irradiation, i.e., 
Number of atoms ý[11]. To perform the simulation, for every MC step, we
first choose a random number 
(between 0 and 1). If
the value of 
is larger than , then a irradiation event, i.e., a ballistic jump will be
chosen. Otherwise, the system will have a thermal jump and the actual
transition 
which occurs can be
found as
.
Figure 5 shows schematically how the transition is chosen. The flowchart of the algorithm is also depicted in
Fig. 6.
Figure 5. Schematic of the probability array composed of segments whose lengths
are proportional to the transition probabilities.
Figure 6. Flowchart of the Kinetic Monte Carlo Algorithm for
simulation 2D immiscible binary alloy under irradiation.
Results and Discussion
1. Simulation
results
We perform the simulation of the A50B50
alloy on the 2D 64×64 square lattice. The melting temperature for the bulk AgCu alloy is about 1053K. Here we use the temperature 
(about 330 K) that is far below
the melting temperature, at which the system goes phase separation.
At the beginning, we don’t introduced the
ion irradiation, i.e., 
. Starting from a random configuration shown in the Fig. 7,
we can see the system evolve into separated phase. The typical configuration of
the system is shown at different time. The configuration at t = 1.52s indicated
a macroscopic phase separation state of the AB alloy.
Figure 7, the evolution of the
phase separation of the AB alloy model with time increasing.
After the first step, we turn on the ion
irradiation on the system. Staring with a macroscopic phase separation, the
system can have different dynamic phase with different irradiation rate . In Fig. 8, we show that start form an arbitrary
configuration of a macroscopic phase separation, like Fig. 8 (a), the system
can come to three kinds of dynamic steady state under irradiation. As
shown in Figure 8, the system under the ion irradiation with a high irradiation
rate, can become miscible, like the Fig. 8 (f),
although the system should have phase separation in its thermal equilibrium
phase diagram. And as we can expect, low irradiation rates are not able to
change the tendency of macroscopic phase separation as shown in Fig. 8 (b) and
(c).
(a) 
, t = 0.0s (b) 
, t = 4.29s (c) 
, t = 4.50s
(d) 
, t = 3.13s (e)
, t = 4.12s (f)
, t = 3.10s
Figure 8. (a) The starting configuration with a macroscopic phase separation. (b) and (c) show the feature of a macroscopic phase separation at a small irradiation rate. (d) and (e) show the feature of a nanoscale compositional patterning at intermediate values of irradiation rate. (f) shows the phase separated alloy becoming miscible at a high rate of irradiation rate.
Beside the above two kinds of states that are similar to those in the thermal equilibrium phase diagram, we found that an additional dynamic state exists in this 2D alloy model under irradiation, similar to what they found in the 3D AB alloy ý[3] ý[11]. We can observe the configurations shown in Fig. 8 (d) and (e) being very similar to that of compositional patterning state shown in Fig 1.
Figure 9 shows the evolution
of the nanoscale patterning dynamic steady states. In
this case, we start from the same configuration Fig. 8 (a), then turn on the
irradiation with and monitor the
changes of configuration in a short time interval.
(a) t = 0.0s (b) t = 0.027s (c) t = 0.096s
(d) t = 0.33s
Figure 9. The evolution of a patterning state with 
, which starts for the configuration(a) the
same as Fig. 8(a).
If we start from a configuration with the
features of solid solutions that is obtained from irradiating the alloy at a
high irradiation rate , we can also get the patterning state by lowering the
irradiation rate. Figure 10 (b) and (c) shows that the configurations obtained
from decreasing rate from to and 
, respectively.
(a) t=0.0s (b) , t = 3.57s (c)
, t = 2.77s
Figure10. Starting form a configuration (a) having
features of solid solution at 
, the patterning
states can also be reached by simply decreasing the irradiation rate to the
intermediate values, like 
(b) and 
(c), which
have given the patterning states starting from a macroscopic phase separation
state.
2.
Discussion
From our simulations of 2D alloy under
irradiation, we find an additional steady state, which have the similar microstructure
of the 3D alloy model under irradiation that have been found and named as nanoscale compositional patterning. Starting either at a
macroscopic phase separation state or a fairly mixed state, the new
compositional dynamic steady state can always be reached as it is shown in Figs.
8 and 10.
In our 2D
simulation, we model the irradiation as the ballistic jumps with the relocation
distance r having an exponential decay distribution which is proportional to
exp(-r/R). As we already pointed out, R usually has a value about ann or 2ann,
where ann is the first nearest neighbor
distance [12]. So we do not need to try many different values of R and explore
the whole dynamic phase diagram theoretically, which are shown in the Figure 1
obtained from 3D KMC simulation. In our investigation, we only set R=ann, the simulation results already shown the
existence of compositional patterning in the 2D AB alloy systems under ion irradiation
with intermediate ballistic jump rate .
To
confirm our findings, we also performed the simulations on a 128×128 2D square
lattice. As shown in Figure 11 (a) and (c), the larger system also has the
patterning states. Because the number of
vacancy in our code is always one, the used in the 64 square
lattice need to be divided by (128/64)2=4 before its being used in
the 128 square lattice for comparison [3, 11]. The small images (b) and (d)
obtained in 64*64 square lattice show similar features of patterning to (a) and
(c) respectively.
From the visualized
configurations of the systems under irradiations with different irradiation rate (we actually run the
simulations with =0 
, 10, 100, 500, 1000, 1500, 2000, 2500, 3000, 5000 and 10000) like those images shown in Fig. 8, the system
seems like undergoing a continually changed length scale of the separated phase
in the 2D systems. However, future work is needed to find a quantitatively
method to describe the patterning state and to locate the boundaries between
the three dynamic steady states.
Figure 11. (a) and (c) the patterning of the 128 square lattices. The
small images(b) and (d) are put on the top of (a) and
(c) for comparison and show similar features of patterning to them respectively.
(a) 128×128 lattice, 
; (b) 64×64 lattice, 
; (c) 128×128 lattice, 
;
(d) 64×64 lattice, 
.
Conclusion
We
performed KMC simulation to investigate the ion irradiation effect on a simple
2D AB alloy, which undergoes phase separation at low temperature in its thermal
equilibrium phase diagram. At high ballistic jump frequency, ion irradiation
can produce a solid solution-liked dynamic steady state. Moreover, a nanoscale compositional patterning steady state is found to
exist in the 2D AB alloy, which consists with the finding of the 3D KMC
simulation in the literature. For the future work, we need to find a
quantitatively method to describe the patterning state and to locate the boundaries
between the three dynamic steady states.
Reference
[1] Russel, K.C., “Phase Stability Under Irradiation,” Progress in Material Science, 28(3-4), pp.229-434, 1984.
[2] Martin, G., “Phase Stability Under Irradiation: Ballistic Effects,” Physical Review B, 30(3), pp. 1424-1436, 1984.
[3] Enrique, R.A. and Bellon, P., “Phase Stability Under Irradiation in Alloys with a Positive Heat of Mixing: Effective Thermodynamics Description,” Physical Review B, 60(21), pp. 14649-14659, 1999.
[4] Enrique, R.A. and Bellon, P., “Compositional Patterning in System Driven by Competing Dynamics of Different Length Scale,” Physical Review Letter, 84(13), pp. 2885-2888, 2000.
[5] Fichthorn K.A and Weinberg W.H., "Theoretical Foundations of Dynamical Monte Carlo Simulation," J. Chem. Phys., 95(2), pp.1090-1096, 1991.
[7] Enrique, R.A. and Bellon, P., “Compositional Patterning in Immiscible Alloys driven by Irradiation,” Physical Review B, 63, pp. 134111 1, 2001.
[8] Wei, L.C. and Averback, R.S., “Phase Evolution During Ion-Beam Mixing of Ag-Cu,” J. Applied Physics, 81(2), pp.613-623, 1997.
[9] Enrique, R.A. and Bellon, P., “Self-organized Cu-Ag nanocomposites synthesized by intermediate temperature ion-beam mixing,” Applied Physics Letter, 78 (26), pp. 4178- 4180, 2001.
[10] Enrique, R.A., Wu, F. and Bellon, P., “A New Approach for the Direct Synthesis of Nanocomposite Thin Films by Ion-Beam Mixing,” Surface and Coatings Technology, 150, pp. 1-7, 2002.
[11]
Abinandanan, T.A., Haider, F. and Martin, G., “Computer Simulations of
Diffusion Phase Transformations:
[12] Enrique RA, Nordlund K, Averback RS and Bellon P, “Simulations of dynamical stabilization of Ag-Cu nanocomposites by ion-beam processing”, Journal of Applied Physics, 93, pp. 2917-2923, 2003.