Dynamic behavior of the antiferromagnetically coupled bilayer Ising model

Using the path probability and lowest approximation of cluster variation method, we study the dynamic and equilibrium properties of a bilayer magnetic system, consisting of two ferromagnetic monolayers antiferromagnetically coupled for different spins $(\sigma=1/2$ and $S=1)$. Firstly, numerical results of the monolayer and total magnetizations are presented under the effect of the diverse physical parameters, and the phase diagrams of bilayer system are discussed. Then, since it is well established that the path probability method is an effective method for the existence of metastable states, the time evolution of monolayer- and total magnetizations is investigated.


Introduction
Ferrimagnets have two or more sublattices that are magnetized in different directions.Ferrimagnetic materials remain magnetic until a critical point called Néel temperature is reached, where thermal motion randomizes the magnetic spins.Spinel ferrites exhibit the most prevalent ferrimagnetic behavior, which is characterized by adjacent antiparallel electron spins with different magnitudes.The possibility that the net moment can change the sign at a compensation temperature below the Néel temperature is intriguing, and they are important for microelectronic and magneto-optical recording [1].The Néel and compensation temperatures are different for each material and span a wide range.Since the magnetic properties of layered composite materials may differ from their corresponding bulk systems, they have received considerable attention.The mixed and pure layered Ising systems with antiferromagnetic exchange interaction provide simple but interesting models for studying layered composite materials [2][3][4][5][6][7][8][9][10][11][12][13][14][15].
The physics of first-order phase transitions is where the idea of metastability initially appeared.From one atom to statistical ensembles of molecules, such as binary fluids, binary alloys, amorphous solids, liquid crystals, minerals, magnetic systems at molecular levels or as a whole, metastability is a common state in physics and chemistry.The quantity of states increases with system size and in case the factors driving their reciprocal interaction are less uniform or more diversified in space.The time-invariance of the active or reactive patterns with respect to the external effects defines stability and metastability in dynamic systems including magnetic systems [18,[20][21][22], electrical circuits, signal trafficking, decisional, neurological and immunological systems.
Investigations of the critical phenomena for the bilayer Ising systems are needed for a basic scientific knowledge on the magnetic properties.The magnetic properties of a ferrimagnetic bilayer system consisting of ferromagnetically coupled monolayers (A and B) with mixed spins (1/2 and 1) were investigated using the exact recursion method [23].It was shown that a mixed spin ferrimagnetic bilayer system may have one tricritical point and compensation points on the recursive lattice in agreement with the previous studies.In recent times, Jiang et al. [24] investigated a bilayer nano-stanene-like structure with Ruderman-Kittel-Kasuya-Yosida (RKKY) coupling described by the Ising model.They found that the exchange coupling, external magnetic field, the number of nonmagnetic layers and the anisotropy had major effects on the magnetization process, specific heat and internal energy.In another study, the magnetic properties, magnetocaloric effect and hysteresis behaviors were investigated and predicted for a mixed-spin AFM/FM Ising BiFe 3 /Co bilayer applying the Monte-Carlo simulation [25].The magnetization, susceptibility and critical temperature were investigated under various exchange couplings and external magnetic fields.The simulation results indicated that a decrease of the exchange coupling and an increase of external magnetic fields can cause an increase of magnetic entropy change, adiabatic temperature change and relative cooling power.In addition to the equilibrium properties using several methods [26][27][28][29][30][31][32][33][34][35], it is important to examine the dynamic properties of the layered mixed system.In this way, by examining the time dependent variations of the monolayer magnetizations, the existence of metastable states as well as stable state can be examined.Thus, our main motivation in this study is to examine the dynamic properties of the bilayer system, whose magnetic properties in equilibrium state on the Bethe lattice were examined before, using the path probability method (PPM).The agreement of the results obtained from PPM will also be compared with the finite temperature results obtained from the lowest approximation of cluster variation method (LACVM).
This article is organized as follows.The bilayer Ising model and examination techniques are introduced in section 2. In section 3, the dynamics of monolayer and total magnetizations, free energy, compensation behavior and phase diagrams are discussed.The conclusions of this study are briefly presented in section 4.

Bilayer Ising model and calculation methods
We consider a bilayer Ising model consisting of two magnetic monolayers (A and B) with different spins,  and .The model is described as a two-layer system seen in figure 1 with spin variables   = ±1/2 and   ′ = ±1, 0 on the magnetic monolayers A and B, respectively.In this way, the Hamiltonian of a two-layer Ising model is given by ( The first two summations in the model Hamiltonian are taken over the nearest neighboring interacting spin pairs on the same monolayers.For the ferromagnetic interactions on each monolayer, the bilinear interaction parameters are considered as   > 0,   > 0 and if the spin interactions between adjacent layers are antiferromagnetic, they are taken as   < 0. Furthermore,  is total numbers of atoms, and  is the coordination number of the lattice.Finally, parameter  measures the strength of uniaxial single-ion anisotropy acting on only B monolayer with spin- sites.For the spin system at equilibrium, we can derive the set of self-consistent equations for the order parameters using the lowest approximation of the cluster variation method (LACVM) [16][17][18].In this method, the average value of each of the spin states can be represented by   1 ,   2 on the sites of monolayer A and   1 ,   2 and   3 on the monolayer B, which are called the state variables.  1 is the fraction of spin value 1/2 on monolayer A and   2 is the fraction of the spin that has the value −1/2 on the same monolayer.Similarly,   1 is the fraction of spin value +1,   2 is the fraction of spin value 0 and   3 is the fraction of spin value −1 on monolayer B. The above variables satisfy the normalization conditions for monolayers A and B, respectively: The problem here is the evaluation of the spin mean values.The order parameters in terms of internal variables are given by Using ( 2)-( 3), the internal variables can be calculated as linear combinations of the order parameters In order to study the equilibrium behaviors of the model, it is necessary to calculate the magnetizations and the quadrupolar moment for each monolayer A and B. In the framework of LACVM [16][17][18], which is identical to the mean-field approximation [19], a simple expression of the internal energy of the model gives the following expression for a 2D square lattice monolayer shown in figure 1: If we substitute (4)-( 5) into (7), the internal energy per lattice sites can be written as: The total entropy  and free energy  are respectively given by where Ω is the weight factor.In this approximation, the monolayer weight factors Ω  and Ω  can be expressed in terms of internal variables for monolayers A and B, respectively, as follows: Using equations ( 7)-( 9) and making use of the Stirling formula, the free energy for per site ( = −/) can be now written as where  is presented for obtaining the normalization requirement. = 1/ B ,  is the absolute temperature and  B is the Boltzmann constant.Thus, the self-consistent equations for the three long-range order parameters, namely   ,   and   are derived by using the following minimization formulae: Using equations ( 10)-( 11), the self-consistent equations or monolayer magnetizations and quadrupolar moment are found to be By solving the nonlinear algebraic equations ( 12) numerically, one can show the thermal variations of the order parameters   ,   and   for various values of   /  , /  and   /  .Examples of the solutions will be given and discussed in the next section.
For a system away from equilibrium, the general time dependence of state variables is given by the PPM [36] as, where   is the path probability rate for the system to go from state  to .Using this technique, we may examine relaxation curves of the order parameters and observe the flatness property of metastable states.The PPM has been successfully applied to describe the dynamic behavior of many homogeneous and non-homogeneous stationary systems.The coefficients    are the product of three factors:    the rate constants with    =   , a temperature-dependent factor which ensures that the time-independent state is the equilibrium state and a third factor which is the fraction of the system that is in the state .A detailed equilibration requires    =   .To determine the path probability rate, the following formula is generally employed where  is the internal energy presented in equation (8).Using equation ( 3), the time dependent rate equations are written as follows: Using equations( 3)-( 7), ( 13)-( 15), the time dependent dynamic equation for monolayer magnetization   is obtained as follows: Then, the same procedure above is applied to the monolayer B for the general Hamiltonian given in (1), the following two dynamic equations are derived for magnetization monolayer B and quadrupolar moment, respectively: where  =  2 / 1 , while the exponential expressions (   and    ′ ) and partition functions of monolayers (  and   ) are calculated from the following general expressions respectively, as: After deriving the dynamical state equations, we can determine their solutions to see the non-equilibrium properties and discuss regarding the results obtained for finite temperature behavior.To this end, ( 16)-( 17) are solved by the fourth-order Runge-Kutta method [37,38].The Runge-Kutta method is a widely used and effective tool to solve the initial value problems of differential equations.It can be used to build precise numeric methods of high-order by functions without the need for higher order derivatives of expressions.

Numerical results and discussion
In this section, we consider the antiferromagnetic interaction between monolayers in all our numerical calculations.When the interaction parameter between A and B monolayers is selected as the antiferromagnetic, the monolayer magnetizations exhibit antisymmetric orientation at ground state and finite temperatures.The variation of monolayer magnetizations with temperature can be easily obtained numerically from the equations (12).The monolayer magnetizations and quadrupolar moment derived from the dynamical equations of the model are identical to the equations obtained by minimizing the free energy with respect to the internal variables in equation (11).In the bilayer Ising model, the different spin magnitudes and interaction parameters at the monolayers significantly affect the magnetic and dynamic properties of the bilayer model.First, let us begin by examining the variation of magnetizations with temperature to analyze the time dependence of monolayer magnetizations and quadrupolar moment.By obtaining the equilibrium state relations of monolayer magnetizations in the time-independent state of the dynamical equations ( d  d = d  d = d  d = 0), their thermal variations can also be investigated.In the absence of the singleion anisotropy parameter /  , the thermal variations of the monolayer magnetizations and quadrupolar moment can be seen in figure 2(a).In this case, the symmetry breaking is evident due to the different ground state spin orientations and the spin interactions between the monolayers.Thus, the bilayer Ising model undergoes a continuous phase transition at Néel temperature   .In figure 2(b), the time-dependent variations, i.e., the dynamics of monolayer magnetizations and quadrupolar moment are given for the same parameters as in figure 2(a).For arbitrarily chosen initial values at the time  = 0, monolayer magnetizations and quadrupole moments relax after sufficient time and become time-independent at a certain fixed temperature as time progresses.If figure 2(a) and figure 2(b) are compared at the temperature  B /  = 1.8, one can see that the equilibrium value of the monolayer magnetizations and the timedependent results obtained from the relaxation curves using the dynamic equations agree exactly.The arbitrary initial values of the monolayer magnetizations were chosen as   = −0.02,  = 0.02 and   = 0.02 for monolayers A and B, respectively.In any case, if only stable states with the lowest free energy exist, the time evolution for all values of monolayer magnetizations will eventually be at the fixed equilibrium states of   = −0.123,  = 0.83 and   = 0.87.On the other hand, if there are some metastable states with local minimum free energy in the bilayer magnetic system at this fixed temperature, the order parameters should also relax at the metastable states, considering all possible initial values.The classification of these states can be done by methods such as free energy analysis, solving dynamic equations and flow diagrams.Since we have seen that the system always relaxes into just one state, the metastable state does not exist in present case.When the time dependence of the layered magnetization is calculated for various  values, it will be seen that there is an inverse relationship between the relaxation time  and parameter .When  2 >  1 , for the same initial values of monolayer magnetizations and quadrupolar moments, the bilayer system relaxes into the stable state faster than  1 =  2 .Here, the first-rate constants are  13 =  32 =  1 which is the insertion or removal of particles associated with the translation of particles through the lattice.The second-rate constant  12 =  2 is associated with the reorientation of a particle at a fixed site.More specifically, the relaxation time takes a long time for the order parameters in the case of  = 1, while the same case gets a shorter time period with increasing values of  2 .Figure 2(c) denotes the dynamic relaxation curves of the order parameters at a fixed temperature  B /  = 3.0 in the paramagnetic phase.Since there are only stable states at the selected reduced temperature, there is a relaxation at   =   = 0 and   = 0.666 for all initial values.In addition, it is seen that the relaxation time shortens with the increase of  2 .
It is known that the interaction parameters in the model Hamiltonian have a significant effect on the equilibrium and dynamic properties.For example, with inclusion of the single-ion anisotropy parameter for the monolayer B, a first-order phase transition is seen for both monolayer order parameters in figure 3(a).In the bilayer Ising model considered, the first-order phase transition is not observed for   /  = −1, but this behavior is observed at sufficiently large values such as   /  = −0.2 and −0.05.Below the first-order phase transition temperature of the monolayer magnetizations, the metastable states should be expected because the metastability is closely related to the first-order phase transitions.Using an appropriate numerical analysis method, the stable, metastable and unstable states of the monolayer order parameters at any temperature can be obtained and classified.In the figure, solid lines indicate the stable states, dashed-dotted lines correspond to the metastable states, and dashed lines also show unstable states.This classification is done by comparing the free energy values of these solutions.The relaxation curves of the monolayer order parameters obtained as a result of dynamic equation solution are shown in figure 3(b).In the solution of dynamic equations, the temperature  B /  = 0.8 is considered for arbitrary initial values.Since unstable states with the highest free energy cannot be determined from the relaxation curves, it is seen that they are stable and metastable states.Using the free energy analysis seen in the inset in figure 3(a), it is determined that the solid lines are stable and the dashed-dotted lines 0,0 0,4 0,8   are metastable states.It should be noted here that if the interaction parameters of monolayers are different from each other, that is, when   ≠   , the thermal variation of magnetizations and quadrupolar moment will be different.This situation is clearly seen for   /  = 0.2 in figure 4(a).In figure 4(b), the dynamic behavior of the order parameters for   /  = −0.2,  /  = 0.2 and  B /  = 0.7 indicates that the possible time evolution of the order parameters is compatible with the stable states observed from the finite temperature behavior seen in figure 4(a).When the temperature is less than the second-order phase transition temperature, the system always comes to a stable state.Therefore, relaxation processes are not influenced by the rate constants and initial values of order parameters.It should be mentioned that when the temperature reaches the continuous phase transition, the system requires an excessively lengthy period to relax into disorder states.On the other hand, if the temperature exceeds critical temperature, it takes short time for the system to relax into disorder states.Phase diagrams of the ferrimagnetic bilayer Ising bilayer system can be obtained from the phase transitions of monolayer magnetizations.The solid lines in the phase diagrams show the continuous phase transitions and the dashed lines show the first order phase transition.Figure 5(a) shows the phase diagram for /  = 0 and −2 in the case of   /  = 1 where the interaction parameters in monolayers are equal.The point where first-and second order phase transition lines join is the tricritical point (TCP).The phase diagram includes the ferrimagnetic region with stable-metastable states and the disordered paramagnetic phase with only stable state.On the other hand, the first-order phase transition line and TCP disappear in the phase diagram for   <   , when /  = 0, −1.5 [figure 5(b)].Next, the equilibrium and dynamic state behaviors of the bilayer Ising model are examined over total magnetization in figure 6.In figure 6(a), the thermal variation of the total magnetization   is given for the fixed   /  = 0.2 value for varying antiferromagnetic interlayer interaction parameters   /  = −1, −0.5, −0.2 and −0.05.The total magnetization exhibits standard N-type phase transition with compensation temperature according to Néel classification for   /  = −1, −0.5, −0.2, −0.05.In the bilayer Ising model, which represents a rather simple model of ferrimagnetic materials, the compensation phenomenon is more pronounced for   /  = −0.05.At the compensation temperature, the total magnetization   vanishes, as a result of cancellation between magnetic moments of different monolayers.The significance of compensation temperature is that it can lead to antiferromagnetic-like dynamics in ferrimagnets and high-speed domain walls.Magnetization switching can be easily achieved with faster domain walls, which makes these materials a promising candidate for high-density data storage.In figure 6(b), the time evolution of total magnetization is illustrated at  B /  = 0.6.As can be seen from the time evolution of total magnetizations, all relaxation curves for   /  = −1, −0.5, −0.2, −0.05 take a fixed magnetization value i.e., time independent stable state shown in equilibrium behavior.

Conclusion
In this study, the bilayer Ising model with different spin magnitudes ( = 1/2 and  = 1) for the monolayers with antiferromagnetically coupled interlayer parameter was examined using the lowest approximation of cluster variation method and path probability method.Within the framework of these methods, the dynamic equations giving the variation of monolayer magnetizations with respect to time and order parameters were obtained.Then, the numerical solution of the dynamical equations was performed by the fourth-order Runge-Kutta method and the variations of magnetizations with respect to time were examined in the case of antiferromagnetic interaction between A and B monolayers.In addition to the stable states, the existence of metastable state in the system has been revealed by solving dynamic equations.
Then, the magnetization relations obtained from the dynamical equations of the model were solved and the behavior of the monolayer and total magnetization with respect to temperature was investigated.From the obtained thermal variations, it is seen that for different values of the interaction parameters, both first-order and second-order phase transformations occur in the system.As with the layer structure with different spin magnitudes in previous studies, the phase transition temperatures of the magnetization are the same, but their value varies significantly with temperature depending on the fact that ferromagnetic interaction parameters in each monolayer with different spin values are equal and different.In addition, an increase in the antiferromagnetic interaction parameter between the monolayers reduces the secondorder phase transition temperatures and the Néel temperature.In the bilayer model, another factor that significantly affects the magnetic behavior of system is the single-ion anisotropy parameter that affects directly the spin-1 states in the monolayer The single-ion anisotropy parameter is especially effective on the phase transition in the bilayer system, that is, the first-order phase transition and the compensation temperature.In the variation of total magnetization with respect to temperature at the compensation temperature, the total magnetization first disappears below the second-order phase transition temperature, and this temperature is of great importance in experimental systems, especially in magnetic recording system.Thus, the different spin magnitudes in the monolayers and inclusion of single-ion anisotropy parameter in the system lead to richer phase transitions and the magnetic behavior such as compensation temperature and metastable state.

Figure 1 .
Figure 1.Part of the bilayer system constructed of two ferromagnetic monolayers (A and B) coupled with an interlayer interaction parameter   .

Figure 5 .Figure 6 .
Figure 5. (Colour online) The transition temperature as a function of the interlayer coupling for values of the single-ion anisotropy when   /  = 1.0 and   /  = 0.2, respectively.The black circle represents the tricritical point (TCP).