Monte Carlo study of an anisotropic Ising multilayer with antiferromagnetic interlayer couplings

TL;DR

This study uses Monte Carlo simulations to analyze the magnetic properties and compensation phenomena of an anisotropic Ising multilayer ferrimagnet with antiferromagnetic interlayer couplings, highlighting parameter effects.

Contribution

It provides detailed Monte Carlo analysis of a complex multilayer Ising model, comparing results with pair approximation and previous studies.

Findings

Identification of compensation and critical temperatures as functions of model parameters

Analysis of parameter influence on the presence of compensation phenomena

Comparison of Monte Carlo results with pair approximation and literature

Abstract

We present a Monte Carlo study of the magnetic properties of an Ising multilayer ferrimagnet. The system consists of two kinds of non-equivalent planes, one of which is site-diluted. All intralayer couplings are ferromagnetic. The different kinds of planes are stacked alternately and the interlayer couplings are antiferromagnetic. We perform the simulations using the Wolff algorithm and employ multiple histogram reweighting and finite-size scaling methods to analyze the data with special emphasis on the study of compensation phenomena. Compensation and critical temperatures of the system are obtained as functions of the Hamiltonian parameters and we present a detailed discussion about the contribution of each parameter to the presence or absence of the compensation effect. A comparison is presented between our results and those reported in the literature for the same model using the pair approximation. We also compare our results with those obtained through both the pair approximation and Monte Carlo simulations for the bilayer system.

Tables

Table 1. Table 1: Results of the fits of the size dependence of T c ​ o ​ m ​ p subscript 𝑇 𝑐 𝑜 𝑚 𝑝 T_{comp} for p = 0.80 𝑝 0.80 p=0.80 , J A ​ A / J B ​ B = 0.65 subscript 𝐽 𝐴 𝐴 subscript 𝐽 𝐵 𝐵 0.65 J_{AA}/J_{BB}=0.65 , and J A ​ B / J B ​ B = − 0.80 subscript 𝐽 𝐴 𝐵 subscript 𝐽 𝐵 𝐵 0.80 J_{AB}/J_{BB}=-0.80 , using Eqs. ( 10 ) and ( 11 ) for different values of L min subscript 𝐿 min L_{{\mbox{\scriptsize{min}}}} . The numbers in parentheses are the statistical errors given by the least-square fit. These fits were made using the data presented in Fig. 7 .

	Fits to $T_{comp}\left(L\right)=a+b{L}^{-\epsilon }$ (Eq. (10))				Fits to $T_{comp}\left(L\right)=a$ (Eq. (11))
$L_{\text{min}}$	$a$	$b$	$\epsilon$	${\chi }^2/n_{DOF}$	$a$	${\chi }^2/n_{DOF}$
20	$3.0669825087\left(7\right)$	$2\left(3\right)\times {10}^3$	$4.35258009946$	$1.8\times {10}^{-2}$	$3.06702996\left(1\right)$	$3.4\times {10}^{-1}$
30	$3.0669845581\left(6\right)$	$4\left(6\right)\times {10}^4$	$5.11999990708$	$1.5\times {10}^{-2}$	$3.067011268\left(3\right)$	$6.6\times {10}^{-2}$
40	$3.0669505972\left(5\right)$	$1.60381821004\left(5\right)\times {10}^{-15}$	$-5.32766666961$	$3.1\times {10}^{-3}$	$3.0669965247\left(2\right)$	$5.4\times {10}^{-3}$
50	$3.0669457156\left(7\right)$	$2.057928392\left(1\right)\times {10}^{-14}$	$-4.78753257781$	$3.8\times {10}^{-3}$	$3.0669970562\left(3\right)$	$6.4\times {10}^{-3}$
60	$3.066935560\left(1\right)$	$1.10365686\left(4\right)\times {10}^{-12}$	$-3.94933526946$	$4.9\times {10}^{-3}$	$3.0669983378\left(3\right)$	$7.7\times {10}^{-3}$
70	$3.066944621\left(3\right)$	$5.11588460\left(2\right)\times {10}^{-14}$	$-4.59247329163$	$7.4\times {10}^{-3}$	$3.0670022149\left(4\right)$	$8.6\times {10}^{-3}$

Equations

$$-\beta \mathcal{H}=\sum _{⟨i\in A,j\in A⟩}{K}_{AA}{s}_i{s}_j+\sum _{⟨i\in A,j\in B⟩}{K}_{AB}{s}_i{s}_j{ϵ}_j+\sum _{⟨i\in B,j\in B⟩}{K}_{BB}{s}_i{s}_j{ϵ}_i{ϵ}_j,$$

$$m_A=\frac{1}{N_A}\sum _{i\in A}{s}_i,$$

$$m_B=\frac{1}{N_B}\sum _{j\in B}{s}_j{ϵ}_j,$$

$$m_{\text{mbox}tot}=\frac{1}{2}\left(m_A+p{m}_B\right).$$

$${\chi }_{\gamma }=N_{\gamma }K\overline{\left(⟨m_{\gamma }^2⟩-⟨|m_{\gamma }|{⟩}^2\right)},$$

$${\bar{f}}_{\text{mbox}}(t,h,L)\sim L^{-(2-\alpha )/\nu }{f}^0(t{L}^{1/\nu },h{L}^{(\gamma +\beta )/\nu })$$

$${\chi }_{\text{mbox}tot}=L^{\gamma /\nu }\mathcal{X}(x_t),$$

$${\frac{d\mathcal{X}(x_t)}{d{x}_t}|}_{T=T_c(L)}=0,$$

$$T_c(L)=T_c+A{L}^{-1/\nu },$$

$$T_{comp}(L)=a+b{L}^{-\epsilon },$$

$$T_{comp}(L)=a=\text{mbox}constant,$$

Full text

I Introduction

The study of ferrimagnetic materials has attracted considerable attention in the last few decades, especially since a number of phenomena associated with these materials present a great potential for technological applications Connell et al. (1982); Grünberg et al. (1986); Camley and Barnaś (1989); Phan and Yu (2007). In such systems the different temperature behavior of the sublattice magnetizations may cause the appearance of compensation points, i. e., temperatures below the critical point for which the total magnetization is zero while the individual sublattices remain magnetically ordered Cullity and Graham (2008). Although unrelated to critical phenomena, at the compensation point there are physical properties, such as the magnetic coercivity of the system, that exhibit a singular behavior Connell et al. (1982); Shieh and Kryder (1986); Ostorero et al. (1994). The fact that the compensation point of some ferrimagnets occurs near room temperature makes them particularly important for applications in magneto-optical drives Connell et al. (1982).

Initially, compensation effects were theoretically studied in bipartite lattices with different spin magnitudes in each sublattice Boechat et al. (2002); Godoy et al. (2004). However, this is not the only possible geometry which may lead to compensation phenomena. In particular, layered ferrimagnets have been extensively studied in the recent past Balcerzak and Szałowski (2014); Szałowski and Balcerzak (2014); Diaz and Branco (2017). In the latter case, systems are composed of stacked planes with different magnetic properties and the realization of antiferromagnetic couplings between adjacent layers has important technological applications such as in magneto-optical recordings Connell et al. (1982), spintronics Grünberg et al. (1986), the giant magnetorresistance Camley and Barnaś (1989), and the magnetocaloric effect Phan and Yu (2007). In addition, the study of the magnetic properties of these systems is of great theoretical interest since it can provide insight into the crossover between the characteristic behavior of two- and three-dimensional magnets.

In recent experimental works, we find examples of realization and study of such bilayer Stier and Nolting (2011), trilayer Smits et al. (2004); Leiner et al. (2010), and multilayer Kepa et al. (2001); Chern et al. (2001); Sankowski and Kacman (2005); Chung et al. (2011); Samburskaya et al. (2013) systems. From the theoretical stance, a bilayer system with Ising spins and no dilution has been studied via transfer matrix (TM) Lipowski and Suzuki (1993); Lipowski (1998), renormalization group (RG) Hansen et al. (1993); Li et al. (2001); Mirza and Mardani (2003), mean-field approximation (MFA) Hansen et al. (1993), and Monte Carlo (MC) simulations Ferrenberg and Landau (1991a); Hansen et al. (1993). A similar system with both Ising and Heisenberg spins has been studied in the pair approximation (PA) both without dilution Szałowski and Balcerzak (2013) and with dilution Balcerzak and Szałowski (2014). There is also a recent work considering an Ising bilayer with site dilution in an MC approach Diaz and Branco (2017). For the multilayer system, the PA has also been applied to the model with Ising and Heisenberg spins, both with no dilution Szałowski and Balcerzak (2012) and with dilution Szałowski and Balcerzak (2014).

If all layers have the same spin (e. g., spin-$1/2$) and we have an even number of layers, it is necessary that different layers have different number of spins in them for the existence of a compensation effect, as discussed in Refs. Balcerzak and Szałowski, 2014; Szałowski and Balcerzak, 2014; Diaz and Branco, 2017. It is also necessary that the layers have asymmetric intralayer exchange integrals and that the layers with stronger exchange integrals have less atoms than their weak-exchange counterparts. That is easy to achieve with site dilution. However, to the best of our knowledge, no numerical simulation methods have yet been applied to multilayer systems with site dilution.

With that in mind, in this work we present a Monte Carlo study of the magnetic properties of a spin-$1/2$ Ising system composed of two kinds of non-equivalent planes, A and B, stacked alternately. All intralayer interactions are ferromagnetic while the interlayer interactions are antiferromagnetic. We also consider the presence of site dilution in one of the kinds of planes. The simulations are performed with the Wolff algorithm Wolff (1989) and with the aid of a reweighting multiple histogram technique Ferrenberg and Swendsen (1988, 1989). The model is presented in Sec. II, the simulation and data analysis methods are discussed in Sec. III, and the results are presented and discussed in Sec. IV.

II Model and observables

The multilayer system we study consists of a simple cubic crystalline lattice such that non-equivalent monolayers (A and B) are stacked alternately (see Fig. 1). The A planes are composed exclusively of A-type atoms while the B planes have B-type atoms as well as non-magnetic impurities. The Hamiltonian describing our system is of the Ising type with spin $1/2$ and can be written as

[TABLE]

where the sums run over nearest neighbors, $\beta\equiv(k_{B}T)^{-1}$, $T$ is the temperature, $k_{B}$ is the Boltzmann constant, the spin variables $s_{i}$ assume the values $\pm 1$, the occupation variables $\epsilon_{i}$ are uncorrelated quenched random variables which take on the values $\epsilon_{i}=1$ with probability $p$ (spin concentration) or $\epsilon_{i}=0$ with probability $1-p$ (spin dilution). The couplings are $K_{AA}>0$ for an AA pair, $K_{BB}>0$ for a BB pair, and $K_{AB}<0$ for an AB pair. The corresponding exchange integrals (see Fig. 1) are given by $J_{\gamma\delta}=\beta^{-1}K_{\gamma\delta}$, where $\gamma=A,B$ and $\delta=A,B$.

For the purpose of the numerical analysis to be discussed in Sec. III, we define some observables to be measured in our simulation. Namely the dimensionless extensive energy is given by $E\equiv\mathcal{H}/J_{BB}$, and the magnetizations of A-type atoms and B-type atoms are, respectively

[TABLE]

where $N_{A}=L^{3}/2$ is the total numbers of A-type atoms in the system and $N_{B}=pL^{3}/2$ is the number of B-type atoms. The total magnetization of the system is

[TABLE]

We also define the magnetic susceptibilities

[TABLE]

where $\langle\cdots\rangle$ is the thermal average for a single disorder configuration whereas the over-line denotes the subsequent average over disorder configurations, $K=\beta J_{BB}$ is the inverse dimensionless temperature, and $\gamma=A,B,{\mbox{\scriptsize{tot}}}$. The total number of atoms in the system is $N_{{\mbox{\scriptsize{tot}}}}=N_{A}+N_{B}$.

III Monte Carlo methods

III.1 Simulational details

We employed the Wolf single-cluster algorithm Wolff (1989) for the MC analysis of Hamiltonian (1) on cubic lattices of size $L^{3}$ with periodic boundary conditions. We performed simulations for linear sizes $L$ from 10 to 100 and for a range of values of the Hamiltonian parameters: $0.0<p\leq 1.0$, $0.0<J_{AA}/J_{BB}\leq 1.0$, and $-1.0\leq J_{AB}/J_{BB}<0.0$. All random numbers were generated using the Mersenne Twister pseudo-random number generator Matsumoto and Nishimura (1998).

For each set of values chosen for the parameters above, we performed simulations in a range of temperatures close to either the critical point or the compensation point. To determine $T_{comp}$ we typically divided the range in 5 to 10 equally spaced temperatures whereas for $T_{c}$ we used from 8 to 17 equally spaced temperatures. In a single simulation for the largest systems considered, i. e., $L=100$, we performed up to $2\times 10^{5}$ steps (Wolff single-cluster updates) and discarded up to $7\times 10^{4}$ steps to account for equilibration. The total number of steps after equilibration was always at least 2000 times the relevant integrated autocorrelation time Newman and Barkema (1999).

To analyze the data generated in the simulations for a particular range of temperatures, we use the multiple-histogram method Ferrenberg and Swendsen (1988, 1989) to compute the observables defined in the previous section at any temperature inside this range. The thermal error associated with those observables is estimated via the blocking method Newman and Barkema (1999), in which we divide the data from each simulation in blocks and repeat the multiple-histogram procedure for each block. The errors are the standard deviation of the values obtained for a given observable for different blocks.

The process is repeated for $N_{s}$ samples of quenched disorder to obtain the final estimate of our observable. We chose $10\leq N_{s}\leq 50$, such that the error due to disorder was approximately the same as the thermal error obtained for each disorder configuration. Finally, we sum both thermal and disorder errors for an estimate of the total error.

III.2 Determining $T_{c}$

In this work we wish to determine the critical point accurately. This is accomplished by means of a finite-size scaling analysis Yeomans (1992), in which we examine the size dependency of certain observables measured for finite systems of several sizes and extrapolate these results to the thermodynamic limit, i. e., $L\rightarrow\infty$. In this approach, the singular part of the free energy density for a system of size $L$, near the critical point, is given by the scaling form

[TABLE]

where $t$ is the reduced temperature, $t=(T-T_{c})/T_{c}$, $T_{c}$ is the critical temperature of the infinite system, and $h$ is the external magnetic field given in units of $k_{B}T$. The critical exponents $\alpha$, $\beta$, $\gamma$ and $\nu$ are the traditional ones associated, respectively, with the specific heat, magnetization, magnetic susceptibility, and correlation length. Various thermodynamic properties can be determined from Eq. (6) and other observables have their corresponding scaling forms, e. g., for the magnetic susceptibility at null magnetic field we have

[TABLE]

where $x_{t}\equiv tL^{1/\nu}$ is the temperature scaling variable.

Eq. (7) provides a powerful method to determine the critical point. It is clear from this scaling law that $\chi_{\mbox{\scriptsize{tot}}}$ diverges at the critical point only in the thermodynamic limit, whereas for a finite system size $L$, $\chi_{\mbox{\scriptsize{tot}}}$ does not diverge but has a maximum at a pseudo-critical temperature $T_{c}(L)$, which asymptotically approaches the real $T_{c}$ as $L$ increases. The maximum occurs when

[TABLE]

which yields the following relation

[TABLE]

where $A$ is a constant, $T_{c}$ is the critical temperature and $\nu$ is the critical exponent associated with the correlation length.

The finite-size scaling method is applicable to different quantities, such as specific heat and other thermodynamic derivatives Ferrenberg and Landau (1991b). We expect the results obtained from the scaling behavior of these other quantities to be consistent, as we were able to verify in preliminary simulations. In this study, however, we focused only on the peak temperatures of the magnetic susceptibilities, defined in Eq. (5), for these peak temperatures occurred fairly close to one another and were the sharpest peaks from all the quantities initially considered.

To determine the pseudo-critical temperatures $T_{c}(L)$, we perform simulations close to the peaks of the susceptibilities and use the multiple-histogram method to obtain $\chi_{A}$, $\chi_{B}$, and $\chi_{\mbox{\scriptsize{tot}}}$ as continuous functions of $T$, as shown in Fig. 2 for the total magnetic susceptibility of a system with $p=0.60$, $J_{AA}/J_{BB}=0.80$, $J_{AB}/J_{BB}=-0.50$, and $L$ from 10 to 100. The location of the peak temperatures is automated using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method Broyden (1970) for each of the $N_{s}$ disorder configurations and the errors are estimated as discussed in Sec. III.1.

Finally, to obtain $T_{c}$, we perform least-square fits with Eq. (9) using the estimates for $T_{c}(L)$ as input. This equation has three free parameters to be adjusted in the fitting process and requires great statistical resolution in order to produce stable and reliable estimates of the parameters. Since the precise estimation of all of these parameters would lead to a drastic increase in computational work, and since for the present work we are not interested in a precise value for the exponent $\nu$, we employ the same procedure presented in Ref. Diaz and Branco, 2017, in which we set a fixed value for the exponent $\nu$ and perform fits with two free parameters, instead of three. These fits are made, for a fixed value of $\nu$, for system sizes not smaller than $L_{{\mbox{\scriptsize{min}}}}$ and the value of $L_{\mbox{\scriptsize{min}}}$ that gives the best fit is located, i. e., the one that minimizes the reduced weighted sum of errors, $\chi^{2}/n_{DOF}$, where $n_{DOF}$ is the number of degrees of freedom. Next, we keep changing the values of $\nu$ and $L_{\mbox{\scriptsize{min}}}$ iteratively until we locate the set of values that globally minimizes $\chi^{2}/n_{DOF}$ and use these values to determine our best estimate of $T_{c}$. This procedure effectively linearizes the fit, although it does not allow for an individual error estimate for the exponent $\nu$.

In Fig. 3 we show examples of these fits of the pseudo-critical temperatures obtained from the maxima of the magnetic susceptibilities for $p=0.50$, $J_{AA}/J_{BB}=0.80$, and $J_{AB}/J_{BB}=-0.50$. In this figure we show the fits using the values of $L_{\mbox{\scriptsize{min}}}$ and $\nu$ that minimize $\chi^{2}/n_{DOF}$. We note that this method gives a very small statistical error for $T_{c}$, even negligible in some cases, but it is important to point out that this error is underestimated when compared to the actual error, obtained through a true non-linear fit. In order to obtain a more realistic error bar, we follow the criterion used in Ref. Diaz and Branco, 2017, in which the values obtained from fits that give $\chi^{2}/n_{DOF}$ up to $20\%$ larger than the minimum are considered in the statistical analysis.

It is also clear from Fig. 3 that the values of $\nu$ obtained through this process are imprecise and unreliable, however, it is worth stressing that it is not our goal to obtain a precise description of the critical behavior for the model. Therefore, the values of $1/\nu$ presented in Fig. 3 serve only as an “effective exponent” used to achieve a good estimate of $T_{c}$.

III.3 Determining $T_{comp}$

The compensation point is the temperature $T_{comp}$ where ${\langle m_{\mbox{\scriptsize{tot}}}\rangle}=0$ while $\langle m_{A}\rangle,{\langle m_{B}\rangle}\neq 0$, as seen in Figs. 4(a) and 5(a). In order to estimate that temperature we perform simulations for a range of temperatures around $T_{comp}$ and obtain the ${\langle m_{\mbox{\scriptsize{tot}}}\rangle}$ values as a continuous function of $T$ using the multiple-histogram method, as show in Fig. 6 for the total magnetization of a system with $p=0.80$, $J_{AA}/J_{BB}=0.65$, $J_{AB}/J_{BB}=-0.80$, and several system sizes $L$ from 20 to 100. The root of $\langle m_{\mbox{\scriptsize{tot}}}(T)\rangle$ is then located using Brent’s method Brent (1973). The process is repeated for $N_{s}$ configurations and the final value and associated error of $T_{comp}$ are determined as discussed in Sec. III.1.

To obtain a final estimate of $T_{comp}$, it is necessary to combine the estimates for different system sizes. In Fig. 5(a) we see that the different $T_{comp}(L)$ are fairly close to one another. However, it is clear from Fig. 6 that the smaller lattices provide somewhat inconsistent results. Fig. 7 shows the size dependence of the compensation temperature estimates obtained from the same data in Fig. 6. We can see that, as $L$ increases, the compensation temperature approaches a fixed value. As the compensation effect is not a critical phenomenon, we have no a priori reason to expect a particular expression for the dependence of $T_{comp}$ on $L$. Based on the form of the curve, however, we propose a power law behavior similar to the one obeyed by the critical temperature:

[TABLE]

where $a$, $b$, and $\varepsilon$ are parameters to be determined in the fitting process. Also similar to the critical temperature, we lack the statistical resolution to determine the three parameters independently; thus, to obtain $T_{comp}$, we resort to the same fitting method described in Sec. III.2 to determine $T_{c}$.

An alternative procedure is employed: we also fit our data to

[TABLE]

for $L>L_{{\mbox{\scriptsize{min}}}}$, which corresponds to averaging the different compensation temperatures considering only the values of $L$ after the $T_{comp}(L)$ curve has approximately converged. The value of $L_{{\mbox{\scriptsize{min}}}}$ is also determined by minimizing the $\chi^{2}/n_{DOF}$ of the fit. The latter method corresponds to the same one used in Ref. Diaz and Branco, 2017.

Table 1 shows the results of the fits presented in Fig. 7 with Eqs. (10) and (11) for different values of $L_{{\mbox{\scriptsize{min}}}}$. We see that values of $a$ obtained from each method are all consistent, irrespective of $\varepsilon$ or $L_{{\mbox{\scriptsize{min}}}}$. At first glance, we could be tempted to choose Eq. (10) with $L_{{\mbox{\scriptsize{min}}}}=40$ as the best fit, based on the value of $\chi^{2}/n_{DOF}$, and use it to obtain the final estimate of $T_{comp}$. Nonetheless, we note that the fit with Eq. (10) cannot be used to determine $T_{comp}$ for any case from Tab. 1 where $L_{{\mbox{\scriptsize{min}}}}\geq 40$, given that $\varepsilon<0$ contradicts our assumption that $T_{comp}(L)$ approaches $T_{comp}$ asymptotically and, in these cases, the limit used to determine the compensation point, i. e., $T_{comp}=\lim_{L\rightarrow\infty}T_{comp}(L)$, does not exist. The fact that these particular values of $a$ are so close to the other $T_{comp}$ estimates is a consequence of the small values obtained for parameter $b$. In fact, in these cases, the term $bL^{-\varepsilon}$ is less than $10^{-4}$, at least for $L\leq 100$. This means that the dependence with $L$ of Eq. (10) becomes irrelevant to the fitting process, i. e., these results become equivalent to a fit using Eq. (11). On the other hand, if we impose $\varepsilon>0$ we obtain the best fits for $\varepsilon$ very close to zero ($\varepsilon\approx 10^{-6}$) for $L_{{\mbox{\scriptsize{min}}}}\geq 40$, which is consistent with a fit using Eq. (11). Yet, the fit with such a small value of $\varepsilon$ produces imprecise estimates for $a$ and $b$, with errors that reach up to three orders of magnitude higher than the parameters themselves.

These results are consistent with the fact that the compensation phenomenon is not in any way related to criticality. Hence, the observables near the compensation point are not supposed to exhibit a power law scaling behavior. Consequently, the method using Eq. (11) turns out to be more adequate to determine $T_{comp}$, besides being simpler and more robust. To estimate the final error bars we follow the same procedure used in Ref. Diaz and Branco, 2017, i. e., we combine the error obtained in the fitting process with both the errors obtained for $N_{s}$ samples and for a single sample via the blocking method.

IV Results and Discussion

The system studied in this work has three parameters in its Hamiltonian, namely the concentration of magnetic sites $p$, and the ratios $J_{AA}/J_{BB}$ and $J_{AB}/J_{BB}$, which represent, respectively, the asymmetry between type-A and type-B intraplanar couplings and the relative strength of interplanar coupling in relation to the type-B intraplanar coupling. All these parameters play important roles in determining the behavior of the system. Our goal is to outline the contribution of each parameter to the presence or absence of the compensation phenomenon. To that end we map out the regions of the parameter space for which the system has a compensation point, as seen in Figs. 4(a) and 5(a), and the regions for which the compensation effect does not take place, as seen in Figs. 4(b) and 5(b).

We start our analysis by fixing the values of $J_{AA}/J_{BB}$ and $J_{AB}/J_{BB}$ and looking at the dependence of both the critical and compensation temperatures with $p$, as seen in Fig. 8, where we plot $T_{c}$ and $T_{comp}$ as functions of $p$ for $J_{AB}/J_{BB}=-1.00$ and for both $J_{AA}/J_{BB}=0.01$ and $J_{AA}/J_{BB}=0.50$. The solid lines are either cubic spline interpolations or linear extrapolations just to guide the eye. The vertical dashed lines mark the characteristic concentration $p^{\ast}$ where the critical temperature and compensation temperature curves meet and below which there is no compensation.

The behavior depicted in Fig. 8 for the multilayer (3d) system is qualitatively the same as displayed by the bilayer (2d) system for both pair approximation Balcerzak and Szałowski (2014) and Monte Carlo Diaz and Branco (2017). This can be seen by comparing our Fig. 8 to Fig. 2 in reference Balcerzak and Szałowski, 2014 and Fig. 7 in reference Diaz and Branco, 2017, since all where made using the same values for the Hamiltonian parameters. As expected, the critical temperatures obtained for the 3d system are consistently higher than those reported for the 2d one using the same approximation Diaz and Branco (2017). In all cases, though, we see that $p^{\ast}$ is higher for $J_{AA}/J_{BB}=0.50$ than it is for $J_{AA}/J_{BB}=0.01$, which indicates that $p^{\ast}$ should increase as the strength of the A intraplanar coupling increases. This is expected, since, for high $J_{AA}/J_{BB}$, the magnetization of planes B will never be the same (in absolute value) as the magnetization of planes A for small $p$. We are able to confirm this trend, for both cases of a weak and a strong interplanar coupling, in Fig. 9, where we show $p^{\ast}$ as a function of the ratio $J_{AA}/J_{BB}$ for $J_{AB}/J_{BB}=-0.01$ and $J_{AB}/J_{BB}=-1.00$. In both cases we see that $p^{\ast}$ increases monotonically as $J_{AA}/J_{BB}$ increases.

In Fig. 10 we have the behavior of $p^{\ast}$ as a function of the ratio $J_{AB}/J_{BB}$ for several values of $J_{AA}/J_{BB}$. Curiously, $p^{\ast}$ is not monotonic as is the case for the 2d system Diaz and Branco (2017). Figure 9 in Ref. Diaz and Branco, 2017 shows the same $p^{\ast}\times J_{AB}/J_{BB}$ diagram for the bilayer system for the particular case of $J_{AA}/J_{BB}=0.50$ and, comparing the latter figure to the 3d results we see that the values of $p^{\ast}$ for the 3d system are consistently lower than those for the 2d one. However, this difference gets smaller as $|J_{AB}/J_{BB}|$ gets weaker ($p^{\ast}$ for the multilayer is $\approx 8\%$ lower than value for the bilayer for $J_{AB}/J_{BB}=-1.00$ and it is $\approx 3\%$ lower for $J_{AB}/J_{BB}=-0.01$), which is consistent with the fact that as $J_{AB}/J_{BB}\rightarrow 0$ the behaviors of both multilayer and bilayer systems cross over to the 2d behavior of non-interacting planes A and B. We also notice that the values of $p^{\ast}$ are less dependent on the values $J_{AB}/J_{BB}$ for 3d than for 2d. For instance, for $J_{AA}/J_{BB}=0.50$, the percentile increase in $p^{\ast}$ as $|J_{AB}/J_{BB}|$ increases from $0.01$ to $1.00$ is $\approx 13\%$ for the bilayer, whereas for the multilayer this increase is only $\approx 7.6\%$. Even if we consider the maximum percentile increase (which for 2d remains $\approx 13\%$ for its monotonic behavior), we have, for the multilayer, at most an increase of $\approx 10\%$ from the minimum at about $J_{AB}/J_{BB}=-0.20$ to the maximum at $J_{AB}/J_{BB}=-1.00$.

Next, for fixed values of $p$ and $J_{AB}/J_{BB}$, we can analyze the influence of the ratio $J_{AA}/J_{BB}$ in the behavior of our multilayer. In Fig. 11 we plot $T_{c}$ and $T_{comp}$ as functions of $J_{AA}/J_{BB}$ for $p=0.70$ (Fig. 11(a)) and $p=0.90$ (Fig. 11(b)), and fixed $J_{AB}/J_{BB}=-0.50$ and $J_{AB}/J_{BB}=-1.00$ for both cases. The vertical dashed lines mark the value of $J_{AA}/J_{BB}$ above which there is no compensation for each case.

In order to contrast the behaviors of 3d and 2d systems, we can compare Fig. 11(a) with Fig. 10 in Ref. Diaz and Branco, 2017, as well as Fig. 11(b) with Fig. 11 in the same reference, since in both pairs of figures we have the same quantities calculated for the same values of fixed parameters. For instance, for $(p,J_{AB}/J_{BB})=(0.70,-1.00)$, the value of $J_{AA}/J_{BB}$ where $T_{c}=T_{comp}$ is $\approx 77\%$ higher for the 3d system, whereas, for $(p,J_{AB}/J_{BB})=(0.90,-0.50)$, this difference drops to $\approx 10\%$. Still, in both cases the region of the diagram for which the system has a compensation point is larger for the multilayer.

Also, in both cases we have the critical temperatures for the 3d system consistently higher than those obtained with MC simulations for the bilayer, as expected. For $(p,J_{AB}/J_{BB})=(0.70,-1.00)$, the critical temperature for the 3d system is almost twice as high as the 2d value at $J_{AA}/J_{BB}=0.0$, while at $J_{AA}/J_{BB}=1.0$, the 3d value is only $37\%$ higher. For $(p,J_{AB}/J_{BB})=(0.90,-0.50)$, the 3d system $T_{c}$ is $26\%$ higher than the 2d one at $J_{AA}/J_{BB}=0.0$ and $25\%$ higher at $J_{AA}/J_{BB}=1.0$. Still consistent with the 3d to 2d crossover, the difference between the 3d and 2d critical temperatures is bigger for a stronger interplanar coupling. However, this difference does not get smaller as $J_{AA}/J_{BB}$ gets weaker, which agrees with the fact that there is no crossover from 3d to 2d as $J_{AA}/J_{BB}\rightarrow 0$. As a matter of fact, we are able to set $J_{AA}/J_{BB}=0$ and the system is still a multilayer as long as $p\neq 0$ and $J_{AB}/J_{BB}\neq 0$.

We can also compare our multilayer MC results with the bilayer PA behavior reported in Ref. Balcerzak and Szałowski, 2014. This is made by contrasting Fig. 11(a) with Fig. 6 in the latter reference, as well as Fig. 11(b) with Fig. 7 in the same reference. After correcting for the difference in temperature scale (due to the PA using $s_{i}\pm 1/2$ and MC using $s_{i}\pm 1$), as discussed in Ref. Diaz and Branco, 2017, for $(p,J_{AB}/J_{BB})=(0.70,-1.00)$ we have the critical temperature for the 3d system in MC $\approx 90\%$ higher than for the 2d system in the PA at $J_{AA}/J_{BB}=0.0$ and $12\%$ higher at $J_{AA}/J_{BB}=1.0$. For $(p,J_{AB}/J_{BB})=(0.90,-0.50)$ the critical temperature for the 3d system in MC is $\approx 1.3\%$ lower than for the 2d system in the PA at $J_{AA}/J_{BB}=0.0$ and $4.7\%$ higher at $J_{AA}/J_{BB}=1.0$. Again we see that the agreement is better for weaker interplanar couplings as in this case both systems are closer to a 2d behavior. It is also worth mentioning that the PA results for the bilayer are much closer to the MC results for the multilayer than to the MC results for the bilayer. This is consistent with the fact that mean-field-like approximations tend to work better in higher dimensions.

In Fig. 12 we have $T_{c}$ and $T_{comp}$ as functions of $J_{AB}/J_{BB}$ for $(p,J_{AA}/J_{BB})=(0.60,0.25),(0.60,0.30)$ (12(a)) and $(p,J_{AA}/J_{BB})=(0.70,0.46),(0.70,0.50)$ (12(b)). The vertical dashed lines mark the value of $J_{AB}/J_{BB}$ below which there is no compensation for each case. Unfortunately no direct comparison can be made between this figure and the bilayer MC results, since the analogous figure for the bilayer (Fig. 12 in Ref. Diaz and Branco, 2017) does not have the same parameters as ours. The reason we did not use the same parameter for our figure 12 as those used for Fig. 12 in Ref. Diaz and Branco, 2017 is that both sets of parameters in the latter figure ($(p,J_{AA}/J_{BB})=(0.7,0.3)$ and $(0.9,0.8)$) are in a ferrimagnetic phase where there is always compensation for $-1.0\leq J_{AB}/J_{BB}<0.0$. Qualitatively, though, it is safe to say that the $T_{c}$ and $T_{comp}$ curves do not get further apart as $J_{AB}/J_{BB}\rightarrow 0$ in the multilayer system as it happens in the bilayer, for both MC Diaz and Branco (2017) and PA Balcerzak and Szałowski (2014). In fact, for fixed $p$ and $J_{AA}/J_{BB}$ we seem to have $T_{c}$ and $T_{comp}$ always very close to one another, and the difference $T_{c}-T_{comp}$ dos not increase monotonically as $J_{AB}/J_{BB}$ increases, having a maximum value somewhere between the point where both curves meet and $J_{AB}/J_{BB}=0.0$.

Finally, as it follows from the analyzes presented above, we can divide the parameter space of our Hamiltonian in two distinct regions of interest. One is a ferrimagnetic phase for which there is no compensation at any temperature and the second is a ferrimagnetic phase where there is a compensation point at a certain temperature $T_{comp}$. We present this results in Fig. 13, where we plot the phase diagrams for concentrations $p=0.5,0.6,0.7,0.8$, and $0.9$. For each concentration, the line marks the separation between a ferrimagnetic phase with compensation (to the left) and a ferrimagnetic phase without compensation (to the right).

These diagrams show that the compensation phenomenon is favored by greater values of $p$, as it is clear that as $p$ approaches $1$, as long as $p\neq 1$, the area occupied by the ferrimagnetic phase with compensation greatly increases. It is also clear that the compensation phenomenon is only present if there is intraplanar coupling asymmetry and, as $p$ decreases, it is necessary to increase the asymmetry for the phenomenon to occur. We also notice that, as $p$ increases, the line separating the phases becomes more vertical, i. e., the presence or absence of compensation becomes less sensitive to the value of $J_{AB}/J_{BB}$. Nonetheless, on all cases presented it is possible to intersect each curve with an actual vertical line in two points, although this is easier to see for smaller values $p$. Therefore, there are sets of values for the parameters for which we can be in a phase without compensation, decrease $|J_{AB}/J_{BB}|$ to get to the phase with compensation, but if we keep decreasing $|J_{AB}/J_{BB}|$ the system goes back to the phase without compensation. This is in agreement with the discussion presented above for Fig. 12 and it is not seen on the behavior of the bilayer for either AP Balcerzak and Szałowski (2014) or MC Diaz and Branco (2017).

It is possible to compare Fig. 13 with the phase diagram for the same system and same parameters obtained with the pair approximation, presented in Fig. 2(a) of Ref. Szałowski and Balcerzak, 2014. We first notice that the area occupied by the ferrimagnetic phase with compensation (to the left of the curve) is consistently greater for MC than for the AP. The lines also become more vertical as $p$ increases in the MC approach and, although the reentrant behavior discussed in the last paragraph is also present in the PA, it is only visible at $p=0.5$ in Fig. 2(a) of Ref. Szałowski and Balcerzak, 2014, while it is clearly present in MC for all values of $p$ presented in Fig. 13. The PA results also indicate that there is no compensation whatsoever for $p=0.6$ and $J_{AB}/J_{BB}$ below $\approx-0.86$ and for $p=0.5$ and $J_{AB}/J_{BB}$ below $\approx-0.36$, whereas the MC results show we always have a phase with compensation for $-1.0\leq J_{AB}/J_{BB}<0.0$ and the concentrations considered in Fig. 13.

V Conclusion

In this work we have applied Monte Carlo simulations to study the magnetic behavior of an Ising multilayer composed of alternated non-equivalent planes of two types, A and B. Both A and B intralayer couplings are ferromagnetic while the interlayer couplings are antiferromagnetic. Additionally, only the B layers are site-diluted. The simulations were performed using the Wolff algorithm and the data were analyzed with multiple histogram reweighting and finite-size scaling methods. Our main goal is to obtain the conditions of existence of a compensation point for the system, i. e., a temperature $T_{comp}$ where the total magnetization is zero below the critical point $T_{c}$.

We studied simple cubic lattices with linear sizes $L\leq 100$ to allow for a precise evaluation of both $T_{c}$ and $T_{comp}$ for a wide range of values of each of the Hamiltonian parameters. The results for the multilayer are compared with those for the bilayer reported in both Monte Carlo Diaz and Branco (2017) and pair approximation Balcerzak and Szałowski (2014) approaches. We see that the compensation phenomenon in the multilayer is favored by small but non-null dilutions and by large intralayer coupling asymmetry, as it is also the case for the bilayer. Similarly, the effect is favored by a weak interlayer coupling, although the sensitivity to this parameter is more pronounced for the 2d system than it is for the multilayer. In addition, we notice that the behavior of the multilayer gets closer to the bilayer as the interlayer coupling gets weaker. This agrees with the crossover that happens at that limit, where the multilayer becomes a set of non-interacting two-dimensional systems.

A summary of the results is then depicted in a convenient way on $J_{AB}\times J_{AA}$ diagrams for several values of site concentration. These diagrams are compared with the PA results from Ref. Szałowski and Balcerzak (2014) for the same model and we emphasize that the MC and PA results are considerably different, both quantitatively and qualitatively.

Work is now underway to accurately determine the critical exponents of the model, as well as to extend the analysis to consider continuous spin symmetry, i. e., Heisenberg spins.

Acknowledgements.

This work has been partially supported by Brazilian Agency CNPq.