Vacancy-assisted domain-growth in asymmetric binary alloys: a Monte Carlo study

A Monte Carlo simulation study of the vacancy-assisted domain-growth in asymmetric binary alloys is presented. The system is modeled using a three-state ABV Hamiltonian which includes an asymmetry term, not considered in previous works. Our simulated system is a stoichiometric two-dimensional binary alloy with a single vacancy which evolves according to the vacancy-atom exchange mechanism. We obtain that, compared to the symmetric case, the ordering process slows down dramatically. Concerning the asymptotic behavior it is algebraic and characterized by the Allen-Cahn growth exponent x=1/2. The late stages of the evolution are preceded by a transient regime strongly affected by both the temperature and the degree of asymmetry of the alloy. The results are discussed and compared to those obtained for the symmetric case.


I. INTRODUCTION
The study of domain-growth is a prototype problem in out-of-equilibrium statistical mechanics. Besides its fundamental interest, it has many technological implications in areas like metallurgy, semiconductors, surface physics, ...etc. Despite during the last 15-20 years domain-growth has been the subject of a deep and continuous investigation [1,2,3], the models commonly used are still far from real materials. A typical experimental situation occurs when a binary alloy undergoing an order-disorder transition is rapidly cooled from the high-temperature disordered-phase to the low-temperature ordered-phase. The growth of the ordered domains, subsequent to the quench, is a complex phenomenon in which different factors may play an important role: details of the atomic diffusion, topological defects, temperature, degeneration of the ordered phase, interface structure, impurities, etc... It is important to realize that the information relative to the influence of each single factor is not always available from the experiments and therefore alternative reliable methods such as numerical simulations may be very useful.
Many Monte Carlo simulation studies based on the Ising-like lattice models assume a very simplified picture of the alloy: i.e., perfect, without defects and with an ordering mechanism based on the neighboring atom-atom exchange. The exchanges proposed are accepted or rejected according to the Metropolis algorithm. The results [4,5,6] obtained are in agreement with the phenomenological Allen-Cahn curvature-driven theory [7] for the time evolution of the mean-size of the ordered domains, i.e. R(t) ∼ t x , with x = 1/2. Moreover, evidence has accumulated that the ordering process exhibits dynamical scaling and that the late-time growth behavior is subjected to a high degree of universality [8]. Concerning the experiments, only a small number have been designed with the aim to check the validity of the growth-law and one may say that, in general, they are compatible with the growth exponent x = 1/2, at least, for quenching temperatures below but close to the equilibrium ordering temperature. In reference [9] the authors report a very complete summary of the experimental data available in the literature.
Concerning binary alloys, atomic diffusion is the underlying microscopic mechanism inherent to any ordering process. Moreover, it is accepted that the atomic diffusion occurs via vacancies. In this sense, quite recently, the possibility for the more realistic vacancy-assisted domain-growth has been investigated [10,11,12]. The model description is based on the three-state ABV lattice model [9] and the microscopic dynamical mechanism accounts for the atom-vacancy exchanges only. The most important result is the prediction for growth exponents greater than the standard Allen-Cahn value.
In this paper we make another step towards the study of real alloys. In most theoretical studies it has been assumed that the alloys are symmetric, i.e. the Hamiltonians used are invariant under the exchange of the two species, A and B, forming the alloy. In the framework of pure Ising models with pair interactions the asymmetry term appear as a field term (chemical potential) which is usually neglected. This turns out to be irrelevant for the models without vacancies since the conservation law for both A and B particles makes the field term to be constant. Nevertheless, for models including vacancies, such asymmetry terms should be considered.
The motivation for this study comes from the fact that in many binary alloys the vacancies exhibit a tendency to locate preferentially in one of the possible different equivalent sublattices (into which the ordered structure can be decomposed) [13]. Our main goal is to modelize such behavior by taking into account the corresponding asymmetry term in the ABV model, and perform Monte Carlo simulations of the vacancy assisted domain growth. Although the present simulations are restricted to a two-dimensional lattice, we expect that the general conclusions will apply to real alloys.
The results obtained may be summarized as follows: The main effect of the asymmetry parameter is to decrease dramatically the speed of the ordering process. The specific interaction between vacancies and interphases is repulsive and, therefore, the vacancy migration from the ordered bulk to the interphases, which is required for the domains to grow, needs to be thermally activated. The associated energy barrier gives rise to a temperature-dependent transient previous to the true asymptotic behavior. Moreover, the process of ordering inside the domains lasts until long times. Nevertheless, we obtain that the long-time domaingrowth behavior is algebraic, and compatible with the Allen-Cahn theory with an exponent x = 1/2 at all the temperatures and degrees of asymmetry studied. This result is different from previous studies of vacancy-assisted dynamics of domain-growth in symmetric binary alloys where the standard x = 1/2 value is found only at temperatures T → T c while, at moderate and low temperatures, the exponent is definitively much larger than 1/2 and tends to 1 for T → 0.
The remaining of the paper is organized as follows. Next, in section II, we briefly summarize the ABV model Hamiltonian and discuss the behavior of the vacancies in terms of the model parameters. In the same section we also define the model dynamics. In section III we present the details of the Monte Carlo simulations and provide the definition of the relevant quantities computed. Results are presented and discussed in section IV. Finally, in Section V we conclude.

A. The ABV Hamiltonian
The ordering configurations of an AB alloy with a constant concentration of vacancies, are described in terms of three state spin-1 variables S i = 1, −1, 0, defined on each lattice site i = 1, 2....N. S i = 1 means that site i is occupied by an A atom, S i = −1 by a B atom and S i = 0 by a vacancy. We shall restrict to 2d-square lattices with a constant number of particles (N A ,N B ) and vacancies (N V ). Considering only nearest neighbor (n.n.) pair interactions, the Hamiltonian for the ABV model can be written as [9,14]: where the sums extend over all n.n. pairs. It contains three independent parameters: J, K, and L. The parameter J determines the ordering of the system. We take J > 0 in order to ensure that the ground state will be antiferromagnetically ordered and formed by two alternating sublattices (A-sublattice and B-sublattice). In what follows we shall take J as the unit of energy. The parameter L * ≡ L/J accounts for the energy difference between AA and BB bonds, i.e.. it controls the asymmetry of the alloy. Finally, the parameter K * ≡ K/J controls the energy of the bonds involving vacancies, and therefore the vacancyvacancy specific interaction. In Table I we have summarized the different bond energies and its corresponding excess-energy with respect to the ordered AB bond. The behavior of the vacancies during the ordering process depends on both parameters K * and L * , and may be understood from the following considerations: 1. For given values of N A , N B (= N A ) and N V , the asymmetry parameter L * controls the tendency for the vacancies to locate preferentially in one of the two sublattices. If, for instance, L * > 1 (L > J) the energy of the BB bond is lower than the energy of the AB bond. The vacancies, in particular those initially located in the A-sublattice, will show a tendency to migrate to the B-sublattice. The same argument holds if L * < −1 for the A-sublattice.
2. Furthermore for |L * | > 1, the tendency for the vacancies to concentrate in a given sublattice, determines the specific interaction between vacancies and antiphase boundaries (APB). The APB's are interphases separating thermodynamically equivalent ordered domains with the same absolute value of the corresponding order parameter. These APB's are always present during the domain-growth process and are formed by a sequence of AA and BB bonds. For values of |L * | < 1, the AA and BB bonds are unfavorable with respect to the AB bond, and therefore the vacancies tend to concentrate at the APB's. But if, for instance, L * > 1 the BB bond is energetically the most favorable and the vacancies, located preferentially in the B-sublattice, will not show any natural tendency to go to the APB's (notice that the AA unfavorable bonds in the APB's can only be broken by vacancies located in the A sublattice). In short, for the asymmetric case (|L * | > 1), the vacancies will prefer to stay in the ordered bulk rather than to concentrate at the APB's, and the corresponding energy difference is proportional to (|L * | − 1).
3. Finally, we notice that the behavior described above may be modified by the specific interaction between vacancies, controlled by the parameter K * . A simple analysis shows that this interaction is attractive for K * < 1 and repulsive for K * > 1. Recently, the effect of K * , in the simplest case L * = 0, has been extensively studied [15,16,17,18].
Here we shall concentrate on the effect of L * only and leave for a further investigation the interplay between K * and L * that may be very subtle.
From the above considerations one finds six different regions in the space of model parameters, as it is illustrated in Fig.1. We notice that the different behaviors of the vacancies are separated by the lines |L * | = 1 and K * = 1.

B. Model dynamics
The non-equilibrium properties of model (1) may be studied using Monte Carlo computer simulation techniques [19]. This requires the implementation of a given microscopic ordering mechanism, compatible with the conservation laws in effect. In our case this means to preserve the number of particles of each kind as well as the number of vacancies. We have chosen the vacancy mechanism: a single vacancy is considered and only atom-vacancy exchanges are allowed. We call this a vacancy jump. Proceeding as it is usual, the vacancy jumps are proposed to n.n. and next nearest neighbors (n.n.n) sites, with equal probability. This is very convenient when performing Monte Carlo simulations in order to prevent trapping phenomena [6,10]. The attempts are accepted or rejected following the Metropolis algorithm which only takes into account the energy difference between the initial and final states in a given vacancy jump, that is: P (∆H) = min{1, exp( −∆H B T )}. The unit of time, the Monte Carlo step (MCS), is defined as N attempts of single vacancy jumps.
In the present study the word "barrier" is used to denote the extra energy needed for the vacancy to migrate from the bulk to the antiphase boundaries. This barrier is intrinsic to model (1) and should not be confused with eventual energy barriers associated to the details of the vacancy path in a given jump between two neighboring sites. There exist in the literature Monte Carlo algorithms which consider the existence of these last barriers [20]. Nevertheless, for the sake of clarity, in this first study of domain growth in asymmetric binary alloys, we will not consider them.

III. SIMULATION DETAILS
Our simulated system is a (nearly) stoichiometric AB alloy containing a single vacancy. Notice that, in this case (N V = 1, N B = N A − 1), the parameter K * is irrelevant and therefore may be neglected in model (1). The particles (N −1) are sitting on a square lattice of size ℓ × ℓ = N, subjected to periodic boundary conditions. The main results presented correspond to L * = 2 and ℓ = 600 (N = 3.6 · 10 5 ). Some additional results for other values of the asymmetry parameter (L * = 0, L * = 1 and L * = 3) are also shown. Starting from a completely disordered configuration, the ordering process is studied at different quenching temperatures (T * ≡ k B T /J = 1.0, 1.25 and 1.5) below the critical ordering temperature T * c ≃ 2.27. The simulations are extended up to ∼ 4 · 10 4 MCS and the results are averaged over ∼ 250 independent runs.
From the information given above it follows that the concentration of vacancies in our simulated system is ∼ 10 −6 . This is reasonable since typical concentrations of vacancies may range from ∼ 10 −3 to ∼ 10 −10 , strongly depending on temperature [21]. For the case with L * = 0 [9], it has been shown that the results obtained in vacancy-assisted dynamics do not change if one increases the number of vacancies (and the system size) maintaining its concentration constant. This is because of the small vacancy concentration we are dealing with which makes the vacancy-vacancy interaction term in eq.(1) to be negligible in front of the other contributions. We expect that the same will apply in the present non-symmetric study and, therefore, focus on the simplest case N V = 1.
From the simulations we have measured the following quantities: 1. Long range order parameter m. It is defined as the absolute value of the antiferromagnetic order parameter on the system.
where the function sign(i) takes alternating ±1 values on the lattice in a chessboard way.
In order to obtain information concerning the path of the vacancy during the ordering process it is convenient to define the two following local order parameters: 2. Local order parameter around the vacancy m v . It is defined as the absolute value of the antiferromagnetic order parameter in a travelling (5 × 5) cell centered at the vacancy position. In order to reduce the numerical uncertainties partial averages over consecutive MCS have been performed. In the case of L * = 2 each value is the average over 10 (consecutive) MCS while for L * = 3 the same has been done over 100 MCS. Furthermore, averages over many (30-1000) independent runs have been performed in order to minimize the dependence with the initial conditions.
3. Local order parameter m (5×5) . It is defined as the absolute value of the antiferromagnetic order parameter in a (5 × 5) cell. The computation is done by dividing the original lattice into cells of size (5 × 5) and then averaging over all of them.
Domain-growth can be monitored by measuring other quantities. We shall study the excess-energy and the width of the structure factor.

Excess internal energy ∆E(t).
It is defined as the excess internal energy per particle: where E(t) is the energy of the system (given by eq. (1)) at time t after the quench and < H > T is the equilibrium energy of the system at the temperature T . At each temperature, this equilibrium energy has been obtained on a system of ℓ = 200, starting from a completely ordered configuration and then using the standard atomatom exchange mechanism. The final value of < H > T is the average over the last 1.8 · 10 4 MCS of a single evolution of 2 · 10 4 MCS long.
5. Structure factor. The structure factor is defined as the Fourier transform of the correlation function. It is written as: where k is a reciprocal space vector, a is the lattice spacing, and r i is the position vector of site i. We have computed the profiles along the (10) and (11) (and equivalent) directions, around the superstructure peak at k = ( 1   2   1 2 ). 6. Structure factor width σ(t). It is defined as the square root of the second moment of the structure factor: where q ≡ | k − ( 1   2   1 2 )| is the distance to the superstructure peak. The sum is performed over the q values in the corresponding direction and extends up to q max ≡ 3σ(t). Note that σ(t) and q max have been obtained in a self-consistent way. By this method we avoid the problems associated to the fitting of a profile function and to the determination of the background which, in this problem, cannot be neglected and shows a time evolution.
Then, provided scaling holds, where R(t) is the mean size of the ordered domains. Finally, in order to reveal the existence of transient regimes during the growth, it is convenient to introduce a definition for the growth exponent, not affected by a priori assumptions: 7. Effective growth exponent x e (t). It is defined as the logarithmic derivative of the excess-energy according to: where, in the computation, instead of the limit, we have taken a small value of ν (ν = 2) and fitted a power-law to all intermediate data contained in the range (t, νt).

IV. MONTE CARLO RESULTS AND DISCUSSION
In this section we present the numerical results obtained from the Monte Carlo simulations. Figure 2 shows snapshots of the configurations obtained at T * = 1, for L * = 0 and L * = 2 at four different selected times. Disordered sites are shown in black. The snapshots for the two values of L * correspond to simulation times for which the system has the same excess internal energy. The overall picture of the domains is qualitatively similar. Nevertheless, for L * = 2 one observes that the amount of disorder in the bulk is larger.
Prior to the study of the growth-law itself, it is interesting to analize the behavior of the vacancy during the ordering process. More precisely, it is important to understand the characteristics of the path followed by the vacancy during the vacancy-assisted ordering process which takes place as a response to the quench, performed at T * . This is done by monitoring the behavior of m v vs. m (5×5) as it is illustrated in Fig. 3 for a system of linear size ℓ = 200, temperature T * = 1.0 and four different values of the asymmetry parameter L * . The straight line denotes a random walk eventually followed by the vacancy in case the path is uncorrelated with the state of order in the system. As it can be seen in Fig.  3, for L * = 2 and L * = 3, the order around the vacancy is clearly higher than the average local order in the system (at all times m v > m (5×5) ). This is indicative of the tendency of the vacancy to stay preferentially in the ordered regions (inside the bulk). This behavior is opposite to the one found for L * = 0 [10,11], where the vacancy clearly exhibits a natural tendency to concentrate at the interphases (APB's). Notice that, after an initial transient, the intermediate case L * = 1 corresponds to the uncorrelated random walk. Moreover, it is clear that the process rapidly slows down as L * increases (especially for |L * | > 1). This can be realized from the number of MCS (indicated by an arrow) needed to achieve the same amount of order in the system (for instance a value of m (5×5) = 0.50). Notice also, that in all cases, at short times, m v > m (5×5) . This is because the ordering mechanism involves the vacancy jumps exclusively and therefore, at short times, the system can only be locally ordered around the vacancy. The above scenario is in full agreement with the discussion given in section II i.e., for |L * | > 1 the specific interaction between the vacancy and the APB is repulsive. In what follows we shall mainly concentrate in the asymmetric case with L * = 2.
We first study the time evolution of the growing domains by following the behavior of ∆E(t). This is shown, in a log-log plot, in Fig. 4 for three different quenching temperatures T * = 1.0, 1.25, 1.50 and ℓ = 600. In order to verify that finite size effects do not mask the long time behavior, we have also measured the corresponding time evolution of the order parameter m(t), which is shown in Fig. 5, for the same three temperatures. Note that simulations have been run up to a time (∼ 4 · 10 4 MCS) for which the order parameter is m ≃ 0.1. This is enough to ensure that, in the long time stages studied, the system is still in the regime of competing domains.
Results in Figures 4 and 5 indicate that the asymptotic regime of ∆E(t) and m(t) [22] is, for the three temperatures, compatible with the algebraic Allen-Cahn growth-law, denoted in the figures by a dotted straight line. We have estimated the growth-exponents of the different curves by least-square fitting a power-law to the last ∼ 25000 MCS. They are indicated in the insets. All the values are compatible with the Allen-Cahn exponent (x = 0.50 ± 0.06). Nevertheless, notice that the values of the exponent obtained from ∆E(t) show a slight tendency to increase as one decreases the quenching temperature. This has to do with the temperature-dependent transient regime, not revealed by the behavior of m(t), that may affect the asymptotic behavior up to very long times. Moreover, from our simulations we have observed that, inside the domains, the equilibrium ordering is established asymptotically and is not complete until long times (see Fig. 2). From all these considerations it follows that, in our case, ∆E(t) has to be taken as an approximate measure of the domain size evolution and that the corresponding exponents are only orientative of the tendency of the system. We now turn back to the transient regime. It is related to the existence of thermally activated processes which hinder the growth of the ordered domains. Indeed, once the domain structure is formed, the domains grow only when the vacancy is located at the APB's. Since for |L * | > 1 the vacancy prefers to stay in the ordered bulk (as it was discussed in section II), the migration to the APB's needs to be thermally activated. The corresponding energy barrier that the vacancy has to climb can be easily evaluated and one obtains ǫ * = 4(|L * | − 1) = 4. This barrier verifies ǫ * > T * , even for the highest temperature studied.
A more quantitative analysis of the whole evolution can be extracted from the computation of the effective exponent obtained from the logarithmic derivative of ∆E(t), as it was defined in the previous section. Results are shown in Fig. 6. The general trends are, for the three temperatures studied, the same. The effective exponent first increases from very small values (at short times) up to a maximum (over x = 1/2) at a time t 0 and next decreases to reach, asymptotically the Allen-Cahn value. On the other hand, the position of the maximum (t 0 ) and the value of the maximum depend on temperature: both are larger as lower the temperature. We have evaluated the value of t 0 for the different temperatures and get t 0 = 1 · 10 3 , 0.59 · 10 3 and 0.64 · 10 3 MCS for T * = 1.0, 1.25 and 1.50 respectively. Moreover, the value of t 0 depends on L * . In fact, we have performed simulations for L * = 3 (not shown in the paper) at the same temperatures, and get t 0 = 4.1·10 4 , 1·10 4 and 0.58·10 4 MCS respectively. The value of t 0 is, for L * = 3 much larger than for L * = 2, providing an idea about the increasing times needed to reach the asymptotic regime when increasing L * . We expect that t 0 is proportional to the passing time of the barrier, We can check that ǫ * ∝ L * − 1 by representing t 0 vs. the factor (L * − 1)/T * in a linearlog plot. This is shown in Fig. 7 for different values of L * . The square corresponds to L * = 1, the diamonds to L * = 2 and the dots to L * = 3. The dashed line shows an exponential fit, which renders ǫ * = 3.7(L * − 1), in agreement with the theoretical estimation given previously, ǫ * = 4(|L * | − 1). The deviation of the points corresponding to the highest temperature (T * = 1.5) may be attributed to the existence of higher order corrections in Eqn. (8) due to entropic effects. We now present the data obtained from the study of the structure factor. In Fig. 8 we show the time evolution of the width of the superstructure peak, σ(t), in the two relevant directions and for the same three temperatures and model parameters as in Figures 4 and 5. The long-time behavior is, for the three temperatures, algebraic (σ ∼ t −x ) and characterized by an exponent x = 0.50 ± 0.05. The inset contains the values of the growth exponent obtained by fitting the algebraic law σ(t) ∼ t −x to the data at each temperature.
Moreover we have studied the existence of dynamical scaling. As an example, in Figure  9 we show the structure factor profiles at T * = 1.0 scaled with the corresponding σ(t) along the two directions. The overlap of the data is very satisfactory. Moreover the shape of the scaling function can be compared with the universal scaling function corresponding to a non-conserved order parameter system [23], which is also plotted with a continuous line. The agreement is excellent for large and intermediate values of q/σ. For small values q/σ < 1 the behavior is slightly different in both directions. This may be associated to the existence of anisotropic correlations in the domain shape. This will be analyzed in detail in a future work.

V. SUMMARY AND CONCLUSIONS
The results presented in the previous section show that during the asymptotic regime of vacancy-assisted domain-growth ∆E(t) −1 ∼ m(t) ∼ σ(t) −1 . In spite of some details observed for small values of q, one may conclude that dynamical scaling holds and that the long-time behavior can be characterized by an unique length which evolves with time according to R(t) ∼ t 1/2 . This asymptotic regime is preceded by a long transient that depends on both the temperature and the asymmetry character of the alloy. In the ABV model the asymmetry is controlled by the parameter L * which in turn determines the energy barrier ǫ * = 4(|L * | − 1) that the vacancy has to climb in its migration from the bulk to the interphases. The characteristics of the transient regime are determined by the ratio ǫ * /T * . This behavior is definitively different from the one found in symmetric binary alloys with L * = 0. The origin of this difference lies in the fact that for symmetric binary alloys the interaction between the vacancies and the interphases is attractive. The first consequence of this is to speed up the global process. The growth exponent was found to be temperature dependent [9] and clearly larger that 1/2 at low and moderate temperatures. Only at T * → T * c the Allen-Cahn exponent is recovered. The fast growth predicted in the symmetric case has never been found experimentally. The authors [10,11] argued that the temperatures at which experiments are usually performed are too high. In the light of the present results it seems more reliable that it has to do with the intrinsic asymmetric character of real alloys. Actually the asymmetry parameter L * can be estimated from ab-initio calculations. We have used the WIEN97 code [24] to evaluate L * for the β-CuZn alloy and have obtained L * ∼ = 3.       FIG. 8. Evolution of the structure factor width σ(t), for a system with L * = 2 and ℓ = 600, along the two directions (10) (continuous lines) and (11) (discontinuous lines). Data corresponds to three different temperatures, as indicated. The dotted straight line indicates the Allen-Cahn behavior σ(t) ∼ t −1/2 . The inset shows the fitted power-law exponents to the last 25000 MCS.