Competing spin dynamics in ising systems

Competing spin dynamics in Ising systems

W. Figueiredo and B.C.S. Grandi

Departamento de Física, Universidade Federal de Santa Catarina

88040-900, Florianópolis, SC, Brazil

[e-mail]

Received 15 October 1999

We study ferromagnetic and antiferromagnetic Ising models in contact with a heat reservoir, and subject to an external source of energy. The contact with the heat bath is simulated by a single-spin flip Glauber dynamics, while the flux of energy is simulated by the two-spin exchange Kawasaki process. Pair approximation and Monte Carlo calculations are employed to determine the phase diagram for the stationary states of the model. We report the results we have obtained in one, two and three dimensions. For instance, in one dimension, while the pair approximation predicts a phase transition for the ferromagnetic case, this is not corroborated by the Monte Carlo simulations. We also use Monte Carlo simulations to evaluate the static and dynamic critical exponents at the transition lines between nonequilibrium steady states. We show that the critical exponents agree with those of the corresponding equilibrium Ising model, for which detailed balance is obeyed.

I Introduction

The study of thermodynamic properties of nonequilibrium systems is a challenging subject because there is still no complete theory to account for these phenomena. That is, formalisms like that of equilibrium thermodynamics, or based on the ensembles of equilibrium statistical mechanics, are still lacking. The behavior of systems out of equilibrium may appear quite complex; many interesting models have been studied. We refer to the books of Marro and Dickman [1] and Privman [2] for a survey in several interacting particle systems on a lattice. For a deeper insight into mathematical questions, and development of the formalism related to phase transitions of interacting particle systems, the books of Ligget [3] and Konno [4] are of fundamental importance. Kinetic Ising models on the lattice are very often used to describe the time evolution and the corresponding steady states of a great variety of interacting particle systems, as for example catalysis, contact process, domain growth, phase separation, and transport phenomena. Aside from a few approximate analytical methods used to investigate these systems, computational methods have been the main tool to study nonequilibrium steady states. In this work we report some studies on kinetic Ising spin systems where the dynamics is dictated by competing stochastic process. Each of the dynamical process we consider satifies detailed balance, which drives the system towards equilibrium states. An interesting situation appears when the system is under theaction of time varying external fields, or when two or more dynamical processes act simultaneously. In these cases, detailed balance is no longersatisfied. The physical properties are very similar to those observed in equilibrium, as for instance the appearance of ordering and phase transitions, but now the system isforced into a nonequilibrium stationary state. The competitive spin dynamics we study here is a combination of single-spin flip Glauber stochastic process [5] and a Kawasaki type two-spin exchange dynamics [6]. Gonzalez-Miranda et al. [7] studied a similar competing system where the Glauber dynamics mimics contact with a heat bath at a fixed temperature, while the Kawasaki dynamics simulates contact with a heat bath at infinite temperature. In this latter case the transition rates are assumed independent of the spin configurations, correlations between sites are absent, and the stationary states are of the Bernoulli type. In the Kawasaki process, the magnetization is conserved by the spin exchanges; the number of microscopic states for a given value of the magnetization is very large. In their work, Gonzalez-Miranda et al. obtained the phase diagram in the T-q plane, where T denotes the temperature and q, the probability of the Kawasaki process, represents a "competition parameter." The phase diagram, calculated through Monte Carlo simulations in two dimensions, shows a line of continuous transitions between the ferro- and paramagnetic phases up to the value q » 0.83 of the competition parameter, where a dynamical tricritical point appears. This result confirmed a previous pair-approximation calculation by Dickman [8] on the same model. In his mean field pair-approximation treatment, he found the value q = 0.72 for the tricritical point. Models of combined Glauber and Kawasaki dynamics were investigated previously by DeMasi et al. [9, 10] for the particular case of an infinitely large value of the Kawasaki transition rate. In this limit, the fast exchange of spins leads to a reaction-diffusion equation for the local magnetization at infinite temperature.

If the two-spin exchange depends on the spin configuration, a new behavior emerges from the competition between Glauber and Kawasaki dynamics: a ferromagnetic Ising system, depending on the value of the competition parameter, can exhibit three different types of magnetic ordering: a ferromagnetic, paramagnetic or antiferromagnetic, as was shown by Tomé and de Oliveira [11]. In their model, the Glauber process simulates the contact with a heat bath at a fixed temperature, and the two-spin Kawasaki exchange simulates a continuous flux of energy into the system. They show that, as the flux of energy is increased, the system goes continuously from the ferromagnetic to the paramagnetic state, and, for a further increase in the energy flux, it goes continuously to an antiferromagnetic stable state. They employed a dynamical pair approximation in their calculations in two dimensions; the results show that the model can be seen as a kinetic Ising model where ferromagnetic Glauber dynamics at a given temperature competes with an antiferromagnetic Kawasaki dynamics at zero temperature. This model is not symmetric to the antiferromagnetic Ising case: by employing the same competitive dynamics for the two-dimensional antiferromagnetic Ising model, the pair approximation gives only a transition line between the antiferromagnetic and paramagnetic phases at low flux of energy [12]. If we increase the energy flux, the only stable states we find are of the paramagnetic type. The ferro- and antiferromagnetic two-dimensional Ising models with competing Glauber and Kawasaki dynamics were also studied by Monte Carlo simulations [13, 14]. The phase diagrams and the corresponding stationary critical exponents were determined. The values found for the critical exponents support the idea that equilibrium and nonequilibrium Ising models, which both exhibit up-down symmetry, belong to the same universality class [15]. For the particular case where the Glauber and Kawasaki dynamics compete at T = 0, in the antiferromagnetic Ising model, the pair approximation predicts ferro, para and antiferromagnetic phases as the strength of the exchange between two nearest neighbor spins in the Kawasaki dynamics changes [16]. Using a quantum formulation of the master equation, Artz and Trimper [17] studied the kinetic Ising model with long-range interactions, subject to competing Glauber and Kawasaki dynamics. They showed that in the exchange-dominated case the system is strongly correlated for each value of temperature.

In Section II we establish the equations of motion for the kinetic Ising model with Glauber and Kawasaki competitive dynamics, and show some results obtained within the dynamical pair approximation. In Section III we present the corresponding results from Monte Carlo simulations, and in Section IV we present our conclusions.

II Equations of motion and the dynamical pair approximation

We consider an Ising model on a hypercubic lattice with N sites. The energy of the system in the state s = (s₁, s₂, ¼, s_N), where the spin variable assumes the values s_i = ±1, is given by

where the summation is only over spins that are nearest neighbors on the lattice, and the cases J > 0, ferromagnetic, and J < 0, antiferromagnetic, are considered. Let P(s, t) be the probability of finding the system in state s at time t. The evolution of P(s, t) is given by the following master equation :

where W(s¢, s) gives the probability, per unit time, for the transition from the state s¢ to state s. In order to take into account the two competing processes, we assume that

In this equation,

is the single-spin-flip Glauber process, which simulates the relaxation of the system towards the equilibrium state at the temperature T, and

is the two-spin exchange Kawasaki process, which mimics the flux of energy into the system. In the above summation, only nearest-neighbor pairs of spins are considered.

In these equations, w_i(s) is the transition probability, per unit time,of flipping spin i, while w_ij(s) is the transition probability, per unit time, of exchanging nearest-neighbor spins i and j. We use the following prescriptions for w_i(s) and w_ij(s):

which is the Metropolis transition rate [18], and

where DE_i is the change in energy when spin i is flipped and DE_ij is the change in energy after exchanging the nearest neighbor spins i and j.

From Eq. (2) it is easy to derive expressions for the evolution of the magnetization, ás_iñ, and for the correlation function between nearest-neighbor spins, ás_js_kñ. They are given by

where

where by (NN of i) we denote a summation over the nearest neighbors of site i.

Although Eqs. (8) - (13) are exact, the mean values of the right-hand sides of these equations cannot be calculated, because we do not know the exact expression for the probability P(s, t). Thus we need to consider an approximate expression for P(s, t). We employ the pair approximation [19] to evaluate the mean values on the right hand sides of Eqs. (10) - (13). We first divide our lattice into two interpenetrating sublattices, and look for solutions such that ás_iñ = m₁ for any spin i belonging to sublattice 1, and ás_jñ = m₂ for any spin j belonging to sublattice 2. The correlation function for any pair of nearest-neighbor spins i and j is written as ás_is_jñ = r. If we then perform these calculations employing the pair approximation, we easily obtain expressions for the time evolution of m₁, m₂ and r. These expressions are very lengthy, and can be found, for the ferromagnetic case, in the work of Tomé and de Oliveira [11].

We look for stationary solutions of Eqs. (8) and (9) characterized by constant values of m₁, m₂ and r. We expect three different solutions: m₁ = -m₂¹ 0 (antiferromagnetic stable states), m₁ = m₂¹ 0 (ferromagnetic stable states), and m₁ = m₂ = 0 (paramagnetic stable states). The equation = 0 , with m₁ = m₂ = 0, gives the paramagnetic stationary states. If we write the order parameter for the ferromagnetic phase in the form , and for the antiferromagnetic phase as , we can expand the right-hand side of Eq. (8) for i = 1 and 2. Retaining only terms linear in m_F and m_A, we can find the linearized equations of motion for these order parameters: their stationary solutions give the transitions between ordered and disordered states. We exhibit, in Fig. 1, the phase diagram obtained by Tomé and de Oliveira [11] for this open ferromagnetic Ising model. The ferromagnetic state is stable for large values of p and at low temperatures; as p decreases, which is equivalent to increase the flux of energy into the system, the ferromagnetic state disappears, giving place to the paramagnetic state. But if the flux of energy is further increased we finally encounter a new ordered state, of the antiferromagnetic type. All the transitions are continuous. A beam of electromagnetic radiation incident on the magnetic system, where the radiation fields interact with the magnetic system through a Hamiltonian involving a nearest-neighbor pair of spins and their neighborhood, can account for the two-spin process we are considering. In this way, the absorption of energy by the system will depend on the local environment [20-23]. Let us consider the system in the stationary state represented by A in Fig. 1. We want to follow its evolution toward the new steady state, B, in the same figure. The magnetic system, in contact with the heat bath at temperature T, is bombarded by a flux q of radiation (the parameter q can give us a measure of the energy flux into the system). Then the system initiates its evolution towards the stationary state B, absorbing energy from the electromagnetic beam, in such a way that, locally, this favors an antiferromagnetic alignment of the spins. Despite absorbing energy from the beam, increasing the antiferromagnetic correlation, the stationary state B is still a globally paramagnetic state. When, finally, the system arrives in the stationary state B, the system becomes transparent to the electromagnetic radiation of the beam (indeed, there occur small energy fluctuations, the energy of the system sometimes increasing due the absorption, at other times decreasing due to the presence of the heat bath). If, for instance, the incident beam is switched off when the system is in steady state B, it will return to the initial state A, by giving up energy to the heat bath; therefore the paramagnetic states A and B are different, because at the steady state B, the antiferromagnetic correlations are higher than in the equilibrium state A. This higher antiferromagnetic correlation is maintained by the constant flux of energy into the system. In order for the system to arrive at steady state C from state B, we must increase the flux of energy, which in our case is equivalent to increasing the parameter q: the system will absorb energy from the beam until it reaches the state C, when will again become transparent to the radiation beam, although the antiferromagnetic correlation and energy are now higher than in state B. For an extremely high flux of energy, that is, q almost equal to 1, the system will reach a fully antiferromagnetic ordered state for whatever value of the heat bath temperature. That is, the correlations induced by the single spin-flip Glauber transitions are immediately destroyed by the interaction of the system with the radiation beam, which returns it to its highest possible energy state per spin, E = 2J. Thus, for a given stationary state, characterized by the parameters T and q, the fluctuation of the energy is due interactions with the heat bath and with the electromagnetic radiation. Therefore, in its steady states, the system becomes transparent to the radiation: antiferromagnetic correlations are created (due the flux of radiation) and destroyed (due the one-spin-flip Glauber transitions). If p = 1, the stationary states are the equilibrium ones, and the critical temperature is given by k_BT_c = 2.885J, which is the equilibrium critical temperature in the Bethe-Peierls approximation.

Figure 1. Phase diagram of the two-dimensional ferromagnetic Ising model with competing dynamics in the pair approximation. T is the temperature, in units of and P = measures the competition between the Kawasaki and Glauber dynamics. F, P and AF denote the ferromagnetic, paramagnetic and antiferromagnetic phases, respectively. A, B and C represent selected steady states.

We present in Fig. 2 the phase diagram of the competing model for antiferromagnetic coupling. This phase diagram was obtained using the dynamical pair approximation to solve the coupled system of equations on the square lattice. It is interesting to note that the stable antiferromagnetic region is very small compared with the ferromagnetic case shown in Fig. 1. The antiferromagnetic phase is destroyed by a small input of energy. The paramagnetic-ferromagnetic transition line is absent in this antiferromagnetic model. As expected, the flux of energy into the system breaks the symmetry between the ferromagnetic and antiferromagnetic Ising models observed in equilibrium.

Figure 2. Phase diagram of the kinetic antiferromagnetic Ising model in two dimensions. T is the heat-bath temperature and P = is related to the flux of energy. The system exhibits only the antiferromagnetic (AF) and paramagnetic (P) phases, separated by a line of continuous non-equilibrium transitions.

In Fig. 3, we exhibit the phase diagram of the one-dimensional version of the ferromagnetic Ising model with competing Glauber and Kawasaki dynamics [24]. In this pair approximation calculation, we see that the phase diagram exhibits a continuous phase transition between the paramagnetic and the antiferromagnetic phases for the nonequilibrium case (p ¹ 1). As is to be expected, for the p = 1, the pair-approximation gives no phase transition at finite temperature for the Ising model. The ferromagnetic ordered state appears only at zero temperature. Similar calculations performed for the one-dimensional antiferromagnetic Ising model shows that for any value of p, and T ¹ 0, the only phase that remains is the paramagnetic one. Although this one-dimensional nonequilibrium Ising model has no exact solution, our Monte Carlo simulations for this model clearly shows that the transition line of Fig. 3 is absent. The phase transition observed in this figure is indeed an artifact of the pair approximation. This approximation underestimates fluctuations and can give wrong results when applied to nonequilibrium one-dimensional situations.

Figure 3: Nonequilibrium phase diagram of the one-dimensional ferromagnetic Ising model, as obtained in the pair approximation. Here h = exp and P = is related to the energy flux. The stationary states of the system are represented by the paramagnetic (P) and antiferromagnetic (AF) phases, separated by a line of continuous nonequilibrium transitions. The singular point P = 0 and h = 0 corresponds to the fully ordered ferromagnetic state.

III Monte Carlo simulations

As we have just seen, Monte Carlo simulations of the nonequilibrium one-dimensional Ising model gave results completely different from the predictions of the pair approximation. In this section, therefore, we present phase diagrams of the ferro- and antiferromagnetic Ising model in two dimensions, obtained by Monte Carlo simulations, in order to compare them with those of Figs. 1 and 2. We also calculated the static and dynamic critical exponents at the steady state phase transitions, and extended the simulations to the three-dimensional version of the competing model we are considering.

Let us first consider the square lattice. We take lattices with L × L = N sites, with L ranging from L = 6 up to L = 80. We used periodic boundary conditions in all of our simulations. We started the simulations with different initial states, in order to guarantee that the final stationary states we use in our calculations are the correct ones. For a given temperature T, and a chosen value of the probability p, we choose at random a spin i, from a given initial configuration. We then generate a random number x₁ between zero and unity. If x₁£ p we perform the Glauber process; in this process, we calculate the value of w_i(s). We generate another random number x₂: if x₂£ w_i(s), we flip spin i, otherwise we do not. If x₁ > p we perform the Kawasaki process. We generate another random number x₃ in order to select one of the four nearest neighbors of spin i, say j. Then we find w_ij, and exchange the selected spins only if w_ij = 1 . We note that the stationary regime is established after 10⁴×N Monte Carlo steps, for all lattice sizes considered. One Monte Carlo step equals N spin-flip or exchange trials. In order to estimate the quantities of interest, we used 5 × 10⁴ Monte Carlo steps to calculate averages.

In order to locate the critical temperature for each value of p, we plotted the reduced fourth-order cumulant [25] U_L(T) (see Eq. (16) below) as a function of temperature T, for several values of L. Once this value is independent of lattice size at the critical temperature T_c, the crossing point of these curves gives T_c [26]. In Fig. 4 we exhibit the reduced fourth-order cumulant for p = 0.5. From this figure we estimate T_c = 2.42 ± 0.01 in units of . We have also considered other values of p in our analysis, in order to determine the complete phase diagram of the model. The resulting phase diagram can be seen in Fig. 5, where we plot h = exp (-) as a function of (1 - p). Clearly, this phase diagram is rather different from that obtained through the dynamical pair approximation, Fig. 1. Here we find a very small region corresponding to the antiferromagnetic phase. The transition between the disordered paramagnetic phase and the ordered antiferromagnetic phase, occurs only for high values of the energy flux. This ordered phase occupies a narrow region of the phase diagram, with p between 0 and p = 0.075. Unlike the work of Gonzalez-Miranda et al. [7], where the transition rate associated with the Kawasaki process was independent of the spin configurations, here we do not observe any dynamic tricritical behavior. On the other hand, the critical temperature exhibits a slight maximum near p = 0.3 . For p < 0.3, where the Kawasaki process is the dominant one, the critical temperature approaches zero as p ® 0. For the pure Kawasaki case, p = 0, the evolution of the system is the same for whatever temperature, that is, we always go to the state of maximum energy compatible with a given initial magnetization. For instance, if the initial state is a paramagnetic one, the final state will be the one where the staggered magnetization per spin reaches its maximum value, that is, 1.

Figure 4. Reduced fourth-order cumulant U_L(T), for p = 0.5, as a function of temperature T for several lattice sizes L. Circles correspond to L = 6, triangles to L = 10, squares to L = 20, crosses to L = 40, downward triangles to L = 60, and diamonds L = 80. We join the data points for each lattice size by a broken line to guide the eye. The critical temperature is T_c = 2.42 ± 0.01, in units of .

Figure 5. Phase diagram of the two-dimensional competing ferromagnetic Ising model. h = exp, and (1-p) is related to the energy flux. F, P and AF refer to the ferromagnetic, paramagnetic and antiferromagnetic phases, respectively.

By employing finite-size scaling relations [25, 26], we can evaluate the stationary critical exponents associated with these transitions. For a system with L × L = N spins, with periodic boundary conditions, we define, at the stationary states, the "magnetization" per spin M_L and the"susceptibility" per spin c_L as

where . We also define the reduced fourth-order cumulant U_L as

The above-defined quantities obey the following finite-size scaling relations in the neighborhood of the stationary critical point:

where , T_c being the critical temperature for each value of p.

Taking the derivative of Eq. (19) with respect to temperature T, we obtain the following scaling relation :

so that . We can then find the critical exponent n from the slope of the straight line which is the best fit to the data points of U_L¢(T_c) for each value of L. Fig. 6 is a log-log plot of U_L¢(T_c) versus L, for p = 0.5: from the best fit to the data we find n = 1.13 ± 0.04. In Fig. 7, we exhibit a log-log plot of the magnetization M_L(T_c) at the critical temperature, T_c, versus L, for p = 0.5. From the slope of the straight line, which is the best-fit to the Monte Carlo data points, and using Eq. (17), we obtain the stationary critical ratio : our estimate in this case is = 0.13 ± 0.01. Another stationary critical exponent of interest is that associated with the susceptibility. We can find the ratio by employing two different approaches based on the scaling relation given in Eq. (18). First, we can construct a log-log plot of c_L(T) versus L, at the critical temperature T_c; then from the slope of the best-fit straight line to the data, we obtain, for p = 0.5, = 1.70 ± 0.05, as can be seen in Fig. 8. We can also estimate the same ratio by a log-log plot of the maximum value of the susceptibility versus L. It is easy to see that if is the value of T for which c_L(T) is maximum, then = T_c + , where u^max is a constant, independent of L, which maximizes c₀(u). Based on these arguments, we immediately see that the maximum of the susceptibility also scales as L. In this way, from Fig. 8, we obtain this = 1.65 ± 0.05 for p = 0.5.

Figure 6:
Log-log plot of U_L¢(T_c) versus L for p = 0.5. The straight line is the best fit to the data, which gives n = 1.13 ± 0.04.

Figure 7: Log-log plot of the magnetization M_L(T_c) versus L, for p = 0.5. From the slope of the straight line, which is the best fit to the data points, we find = 0.13 ± 0.01.

Figure 8: Log-log plot of the susceptibility c_L(T) versus L. Circles: c_L(T) at T = T_c; triangles: c_L(T) at its maximum. The straight lines are best fits to the data points. From the slopes we obtain = 1.70 ± 0.05 (circles), and = 1.65 ± 0.05 (triangles).

In Fig. 9 we exhibit the behavior of n as a function of (1 - p). For this Ising model with competing Kawasaki and Glauber dynamics, the stationary critical exponent n is almost equal to 1. This interesting behavior is in agreement with the arguments given by Grinstein et al. [15] that equilibrium and nonequilibrium stochastic spin systems which present up-down symmetry fall in the same universality class. In order to corroborate these arguments, we also calculate the stationary critical exponent ratios and from the log-log plots of the magnetization M_L and of the susceptibility c_L, respectively, as functions of L. For each log-log plot we calculated M_L and c_L at the critical temperature T_c(p), as exhibited in Fig. 5. In Figs. 10 and 11 we exhibit our results for and , respectively, for several values of p. The exact values for the equilibrium Ising model exponents are well known: n = 1, b = 1/8 and g = 7/4. As we can see, our estimated values for and , reported on Figs. 10 and 11, are in accord with the corresponding values at equilibrium. We have also calculated the critical exponents at the continuous transition line between the paramagnetic and antiferromagnetic phases. For p = 0.03, for instance, we found the following values:T_c = 6.82 ± 0.02, n = 0.97 ± 0.05, = 0.13 ± 0.01 , and = 1.83 ± 0.06 .

Figure 9.
Critical exponent n as a function of the parameter (1-p) at the ferromagnetic- paramagnetic transition line of Fig. 5.

Figure 10. Ratio as a function of the parameter (1-p) at the ferromagnetic- paramagnetic transition line of Fig. 5.

Figure 11. Ratio as a function of the parameter (1-p) at the ferromagnetic- paramagnetic transition line of Fig. 5.

In Fig. 5, at higher temperatures, we only see a transition between the paramagnetic and the antiferromagnetic phases. When the Glauber dynamics dominates, only paramagnetic stationary states are possible. However, as the internal energy of the system increases, which is equivalent in our description to the predominance of the Kawasaki dynamics, long-range antiferrromagnetic order appears. A detailed Monte Carlo analysis along with finite-size scaling arguments, at infinite temperature, shows that the transition occurs for p_c = 0.073 ± 0.003, that is, the system is driven almost exclusively by Kawasaki dynamics [27]. Within the accuracy of our Monte Carlo statistics, the value we found is the same as the one obtained for the isotropic majority-vote model on a square lattice [28]. In this model a spin variable s_i = ±1 is chosen at random on a square lattice at each discrete time. The chosen spin adopts, with probability p, the sign of the majority of its four nearest neighboring spins, and the sign of the minority with probability (1 - p). As the number of nearest neighbors is even, the chosen spin flips with probability 1/2 when there are equal numbers of positive and negative spins in its neighborhood. The only transitions in the isotropic majority-vote model are single-spin flips. As was shown by de Oliveira [28], his model can be regarded as composed of two processes: a noiseless zero-temperature process, and an infinite-temperature one, where the spins are distributed at random. In this sense, our model shows some resemblance to the isotropic majority-vote model, where only single-spin flips are considered. The noiseless zero-temperature process in that model is replaced here by a spin-exchange Kawasaki dynamics at an infinitesimally negative temperature; on the other hand, the effect of the infinite-temperature process in the majority-vote model is equivalent in our case to a Glauber dynamics with a transition rate of 1/2. From the simulation data, we obtained the critical exponents at the transition line between the paramagnetic and the antiferromagnetic phases, at infinite temperature, for this nonequilibrium stochastic Ising spin system. The log-log plot of U_L^¢ (p_c) versus L gives n = 1.0 ± 0.1. In addition, the log-log plot of the staggered magnetization, M_L(p_c) versus L, at the critical probability p_c, gives us the ratio ; the estimated value obtained from the slope of the best-fit straight line to our data is = 0.13 ± 0.01. The stationary critical exponent associated with the susceptibility can be found by employing two different approaches, as was explained before: we can construct a log-log plot of c_L(p) versus L, at p_c; then, from the slope of the straight line, which is the best fit to the data points, we obtain the value = 1.71 ± 0.01. We can also estimate this ratio from a log-log plot of the maximum value of the susceptibility versus L. It is easy to see that if p_L^max is the value of p for which c_L(p) is maximum, then p_L^max = p_c + , where u^max is a constant, independent of L, which maximizes c₀(u). Based on these arguments, we immediately see that the maximum of the susceptibility also scales as L. In this way, we obtain = 1.70 ± 0.01. The values we have found for the stationary critical probability p_c and for the set of critical exponents are in agreement with those obtained by de Oliveira [28] for the isotropic majority-vote model on the square lattice. As we have pointed out, this fact must be related with the similarities between the competing processes in the two models.

We have also determined the dynamical critical exponent z for this competing Glauber and Kawasaki dynamics [30]. Following Suzuki [29], dynamic finite-size scaling theory asserts that the magnetization of a system of linear size L, at its critical point, evolves in time according to the scaling relation

If we consider very large lattices, it is expected that the magnetization does not depend on the lattice size. Then, it is easy to see that M(t,L) can be written as

where A is a constant that does not depend on L, and the last equation is valid only for large values of L. Then, taking into account the last equation, we can evaluate the exponent z froma log-log plot of M(t, L) versus t, for a fixed lattice size L, once we know , given in Fig. 10.

The Monte Carlo method was used again to follow the evolution of the magnetization in time for the competing model we are studying. First of all, we select a given value of competition parameter p. For this value of p we read in Fig. 5 its critical temperature, corresponding the transition between the ferromagnetic and paramagnetic phases. After initializing the system in its ground state, it evolves in time (measured in Monte Carlo steps); we record the magnetization at intervals of 10 Monte Carlo steps.

In Fig. 12 we exhibit a log-log plot of M(t) versus t for two lattice sizes, L = 160 and L = 320, for p = 0.5. We see that the decay of M(t) is almost independent of L, which allows us to use Eq. (22) to evaluate the critical exponent z. In these calculations we used 100 and 50 samples for the small and large lattices, respectively. We followed the decay of magnetization up to 320 Monte Carlo steps. If we discard the first fifty Monte Carlo steps, we can fit our data points to a straight line, and we obtain from its slope the value = 0.065 ± 0.001. We have discarded the initial data points because we want to study the system in its second regime, where a power-law decay of the magnetization is expected [31]. Taking (from Fig. 10) our result = 0.13 ± 0.01 for p = 0.5, we estimate z = 2.0 ± 0.2.

Figure 12.
Magnetization as a function of time, measured in Monte Carlo steps (MCS) for p = 0.5. Measurements were made every 10 MCS, between 10 and 320 MCS. Lattice sizes: (160 ×160) (small circles) and (320 × 320)(small triangles).

In Fig. 13 we show a plot of the exponent z as a function of (1 - p). For all values of p we used lattices of size L = 320, and runs of up to 320 Monte Carlo steps. As we can see, the values of the dynamical critical exponent fluctuate around z = 2. This indicates that the underlying symmetries of this model are not affected by the energy flux. We would like to point out that since in our simulations the magnitude of our estimated errors can be as large as 0.2, the above assumption can not be taken as a rigorous statement. But we expect that for whatever value of the energy flux, the critical exponent z must remain the same, because the intensity of the flux of energy cannot change the universality class of this competing model. For the special case of p = 1, that is, when the system satisfies detailed balance, we found z = 2.0 ± 0.1. This value must be compared with the best estimates for this exponent obtained through intensive use of large-scale simulations: for instance, Stauffer [32] found z = 2.18 for a square lattice with L = 496640 and Linke et al. [33] found z = 2.16 for a square lattice with L = 10⁶. For the continuous transition between the paramagnetic and antiferromagnetic phases, the dynamical critical exponent can be obtained in a similar manner as we have done for the ferromagnetic to paramagnetic transition. For instance, for p = 0.03, we found z = 2.0 ± 0.2. In particular, for the transition at infinite temperature, where the critical probability is p_c = 0.073, the value of the dynamical critical exponent z is 1.9 ± 0.1.

Figure 13.
Dynamical critical exponent z as a function of (1 - p) at the transition line between ferro and paramagnetic phases. The estimated values of z fluctuate around the value 2.0.

For J < 0 (antiferromagnetic coupling), the phase diagram for competing Glauber and Kawasaki dynamics was also obtained through Monte Carlo simulations [14]. The phase diagram is very similar to the one found in the dynamical pair approximation (Fig. 2). At T = 0, for instance, the critical value of p is 0.965 while in the pair approximation its value is 0.968. We also evaluated the critical exponents along the antiferromagnetic-paramagnetic transition line [14]: as for the ferromagnetic case, the critical exponents compare very well with the analogous static exponents of the corresponding two-dimensional equilibrium Ising model.

Finally, we studied the three-dimensional version of the competing Glauber and Kawasaki dynamics for the ferromagnetic Ising model [34]. We used Monte Carlo simulations to find the phase diagram for the stationary states and the corresponding critical exponents. The phase diagram is similar to that observed in two-dimensions (Fig. 5): although the stationary states are mainly ferromagnetic at low temperatures, an antiferromagnetic phase appears for extremely high values of the flux of energy. The region of the phase diagram occupied by the antiferromagnetic phase is larger than in the two-dimensional case. For instance, at very high values of temperature the critical value for the paramagnetic to antiferromagnetic transition is p = 0.3, while the value in two dimensions model is p = 0.073. At zero temperature, we observe only a ferromagnetic steady state, for any value of the competition parameter. We have also calculated the static critical exponents as a function of the competition parameter p. For p = 0.5, for instance, we find n = 0.67 ± 0.04, = 0.53 ± 0.01 and = 1.93 ± 0.03. These exponents are almost independent of p. We would like to stress that the values we have obtained for these critical exponents compare very well with the analogous static exponents of the equilibrium three-dimensional Ising model. For instance, numerical investigations of the equilibrium model by Ferrenberg and Landau [35] yield n = 0.6289 and = 0.518. As this nonequilibrium model preserves up-down symmetry, it is expected that it belongs to the same universality class as the equilibrium Ising model [15].

IV Conclusions

We have studied the behavior of ferromagnetic and antiferromagnetic Ising models subject to two competing stochastic dynamics. The system is contact with a heat bath at fixed temperature, and under the action of an external energy flux. The exchange of energy with the heat reservoir is assumed to be represented by the single-spin-flip Glauber process, while the energy flux into the system is simulated by a kind of two-spin-flip Kawasaki diffusive process. We considered the model in one, two and three dimensions, and have employed the pair approximation and Monte Carlo simulations in our analyses. In one dimension, while the pair approximation predicts a phase transition between the paramagnetic and antiferromagnetic phases for the ferromagnetic Ising model, the phase transition does not appear in the simulations. For the square lattice, in the ferromagnetic case, the predictions the of pair approximation for the phase diagram differ from simulation results. For example, at T = 0, the simulations give only a stable ferromagnetic phase for any value of the competition parameter, while the pair approximation predicts a paramagnetic phase, between the ferro- and antiferromagnetic phases. For T ¹ 0, three different phases appear in the Monte Carlo simulations; however, the antiferromagnetic phase occupies only a small area of the phase diagram. In the limit of very large temperatures, the model becomes similar to the isotropic majority-vote model on a square lattice. For this ferromagnetic model the static and dynamic critical exponents were calculated at the steady-state transition lines; the values we found are in accord with those of the equilibrium Ising model, where detailed balance is obeyed. On the other hand, for the competing two-dimensional antiferromagnetic Ising model, both the pair approximation and Monte Carlo calculations give only a continuous transition line separating the antiferromagnetic and paramagnetic phases. At zero temperature, the critical value of the competing parameter is essentialy the same in both calculations. Also, the critical exponents along the antiferro-paramagnetic transition line are those expected for the equilibrium Ising model in two dimensions. We have also studied the three-dimensional version of the competing Glauber and Kawasaki dynamics for the ferromagnetic Ising model. Our simulations gave a phase diagram which is similar to that observed for the two-dimensional ferromagnetic Ising model. At zero temperature the stable states are of the ferromagnetic type, regardless of the value of the competition parameter; the antiferromagnetic phase occupies a larger region of the phase diagram than in the two-dimensional case. The static critical exponents were calculated along the transition lines. The values we found are independent of the competition parameter, and agree with the static exponents of the equilibrium three-dimensional Ising model. As the Glauber and Kawasaki dynamics preserve the up-down symmetry of the Ising model, even when both operate simultaneously, it is expected that the resulting nonequilibrium model belongs to the universality class of the equilibrium Ising model.

Acknowledgements

We are very grateful to Tânia Tomé and Ron Dickman for the careful revision of the manuscript. This work was partially supported by the Brazilian agencies CNPq and FINEP.

• [1] J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999).
• [2] V. Privman, ed., Nonequilibrium Statistical Mechanics in One Dimension (Cambridge University Press, Cambridge, 1996).
• [3] T. M. Liggett, Interacting Particle Systems (Springer-Verlag, New York, 1985). 
• [4] N. Konno, Phase Transitions of Interacting Particle Systems (World Scientific, Singapore 1994). 
• [5] R. J. Glauber, J. Math. Phys. 4, 294 (1963). 
• [6] K. Kawasaki, in Phase Transitions and Critical Phenomena, ed. C. Domb and M. S. Green (Academic, London, 1972), vol. 4.
• [7] J. M. Gonzalez-Miranda, P. L. Garrido, J. Marro and J. L. Lebowitz, Phys. Rev. Lett. 59, 17 (1987). 
• [8] R. Dickman, Phys. Lett. A 122, 463 (1987). 
• [9] A. DeMasi, P. A. Ferrari and J. L. Lebowitz, Phys. Rev. Lett. 55, 1947 (1985). 
• [10] A. DeMasi, P. A. Ferrari and J. L. Lebowitz, J. Stat. Phys. 44, 589 (1986). 
• [11] T. Tomé and M. J. de Oliveira, Phys. Rev. A 40, 6643 (1989). 
• [12] B. C. S. Grandi and W. Figueiredo, Phys. Rev. B 50, 12595 (1994). 
• [13] B. C. S. Grandi and W. Figueiredo, Phys. Rev. E 53, 5484 (1996). 
• [14] B. C. S. Grandi and W. Figueiredo, Phys. Rev. E 56, 5240 (1997).
• [15] G. Grinstein, C. Jayaparash and Yu He, Phys. Rev. Lett. 55, 2527 (1985).
• [16] Yu-qiang Ma, Ji-wen Liu and W. Figueiredo, Phys. Rev. E 57, 3625 (1998). 
• [17] S. Artz and S. Trimper, Int. J. of Mod. Phys B 12, 2385 (1998). 
• [18] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953).
• [19] H. Mamada and F. Takano, J. Phys. Soc. Jpn. 25, 675 (1968). 
• [20] S. M. Rezende, Phys. Rev. B 27, 3032 (1983). 
• [21] K. C. Johnson and A. J. Sievers, Phys. Rev. B 10, 1027 (1974). 
• [22] P. Moch, G. Parisot, R. E. Dietz and H. J. Guggenheim, Phys. Rev. Lett. 21, 1596 (1968). 
• [23] P. A. Fleury, S. P. S. Porto and R. London, Phys. Rev. Lett. 18, 658 (1967). 
• [24] B. C. S. Grandi and W. Figueiredo, Mod. Phys. Lett. B 10, 945 (1996). 
• [25] V. Privman, Finite-Size Scaling and Numerical Simulation of Statistical Systems (World Scientific, Singapore, 1990). 
• [26] K. Binder, Z. Phys. B 43, 119 (1981). 
• [27] B. C. S. Grandi and W. Figueiredo, Physica A 234, 764 (1997). 
• [28] M. J. de Oliveira, J. Stat. Phys. 66, 273 (1992). 
• [29] M. Suzuki, Prog. Theor. Phys. 58, 1142 (1977). 
• [30] B. C. S. Grandi and W. Figueiredo, Phys. Rev. E 54, 4722 (1996). 
• [31] Z. Li, L. Schülke and B. Zheng, Phys. Rev. E 53, 2940 (1996). 
• [32] D. Stauffer, Physica A 244, 344 (1997). 
• [33] A. Linke, D.W. Heermann, P. Altevogt and M. Siegert, Physica A 222, 205 (1995).
• [34]  J. R. S. Leăo, B. C. S. Grandi and W. Figueiredo, Phys. Rev. E 60, (1999). 
• [35] A. M. Ferrenberg and D. P. Landau, Phys. Rev. B 44, 5081 (1991).

Publication Dates

• Publication in this collection
  17 Oct 2001
• Date of issue
  Mar 2000

History

• Received
  15 Oct 1999