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. 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. ͓S0163-1829͑99͒09329-7͔


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 the fact that during the last 15-20 years domain growth has been the subject of a deep and continuous investigation, 1-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 Isinglike 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 curvaturedriven 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. 8Concerning 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 Ref. 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.1][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͒. 13Our 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 domain growth 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 Sec.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 Sec.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 Sec.IV.Finally, in Sec.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 ͑NN͒ pair interactions, the Hamiltonian for the ABV model can be written as 9,14 HϭJ ͚ ͗i j͘ where the sums extend over all NN 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 vacancy-vacancy specific interaction.In Table I we have summarized the different bond energies and 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.6][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 nonequilibrium properties of model ͑1͒ may be studied using Monte Carlo computer-simulation techniques. 19his 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 NN and next-nearestneighbors ͑NNN͒ sites, with equal probability.This is very convenient when performing Monte Carlo simulations in order to prevent trapping phenomena. 6,10The attempts are ac-  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. 20Nevertheless, 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 lϫlϭN, subjected to periodic boundary conditions.The main results presented correspond to L*ϭ2 and lϭ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. 21For 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͒ negligible in front of the other contributions.We expect that the same will apply in the present nonsymmetric 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 sgn(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 traveling (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 ͑con-secutive͒ 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 l ϭ200, starting from a completely ordered configuration and then using the standard atom-atom 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 analyze 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 lϭ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 Sec.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 lϭ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 Figs. 4 and 5 indicate that the asymptotic regime of ⌬E(t) and m(t) ͑Ref.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 ϳ25 000 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 Sec.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, t 0 ϳe ⑀*/T* .͑8͒ We can check that ⑀*ϰL*Ϫ1 by representing t 0 vs the factor (L*Ϫ1)/T* in a linear-log 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 Eq. ͑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 Figs. 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 Fig. 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 correspond- ing to a nonconserved 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 than 1/2 at low and moderate temperatures.Only at T*˜T c * the Allen-Cahn ex- ponent 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. 25 cepted 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/k B T)͖.The unit of time, the Monte Carlo step ͑MCS͒, is defined as N attempts of single vacancy jumps.

FIG. 1 .
FIG. 1. Space of the parameters K* and L* divided into the regions where different dynamics are expected.The solid line K* ϭ1 separates the region of vacancy attraction (K*Ͻ1) from the region of vacancy repulsion (K*Ͼ1).The solid lines L*ϭϮ1 separate the region where the vacancies tend to precipitate at the antiphase boundaries (͉L*͉Ͻ1) from the region of vacancyantiphase boundary repulsion (͉L*͉Ͼ1).

FIG. 2 .
FIG. 2. Sequence of snapshots of the evolving domain structure for a system with lϭ600, T*ϭ1.0, and two values of L*.Black regions represent disordered sites.

FIG. 3 .
FIG. 3. Evolution of the local order parameter defined around the vacancy m v vs the local order parameter m 5ϫ5 for a system with lϭ200, T*ϭ1.0, and different values of L*, as indicated.The dotted straight line shows the behavior corresponding to a random walk of the vacancy.The arrows indicate, for each case, the number of MCS necessary to reach a value m 5ϫ5 ϭ0.5.

FIG. 6 .
FIG. 6. Effective growth exponent as a function of time obtained from the logarithmic derivative of the behavior of ⌬E(t), corresponding to the same cases as in Fig. 4. Estimations of t 0 and typical error bars are shown.

TABLE I .
Absolute bond energies and excess energies with respect to the AϪB bond for the ABV model as a function of the parameters J, K, and L.