Monte Carlo study of the magnetic properties and magnetocaloric effect of an AFM/FM BiFeO₃/Co bilayer

BiFeO₃ belongs to the family of multiferroic materials and also encompasses ceramic materials. It can play a significant role in the application of magnetoelectric read heads [1], magnetoelectric transformers [2], memory devices [3], and tunable microwave devices [4]. In practical applications, BiFeO₃-based components, such as the BiFeO₃/Co bilayer, which can be synthesized by rapid hot pressing [5], solid-state reaction methods [6], the Pechini method [7], and so on, have attracted much attention. In addition, much effort has been devoted to the study of the magnetic, as well as electrical properties of the BiFeO₃/Co bilayer based on diverse methods from theoretical and experimental aspects [8–10].

In recent years, research on the magnetocaloric effect (MCE), which was first proposed in 1881 by Warburg [11], is attractive because of its potential application in sustainable energy-related emerging technologies, for instance, magnetic refrigeration [12–15]. As far as we know, MCE is an intrinsic and exciting nature of the magnetic system, which can be defined by the fact that the temperature of a ferromagnet varies with the strength of the magnetic field in an adiabatic process [16, 17]. From a physical point of view, MCE is directly related to the coupling between magnetic atoms in the magnetic sub-lattice and the applied magnetic field [18]. It is worth mentioning that isothermal entropy change (−ΔS), adiabatic temperature change (ΔT_ad), and relative cooling power (RCP) are representative characteristics that can better describe MCE, and these parameters are necessary to obtain a high magnetic refrigeration efficiency [19]. There are many investigations on the experimental and theoretical study of MCE. Experimentally, there are many reports of MCE on diverse nanoparticles or compounds. To name a few, the Ln_0.67Sr_0.33MnO₃ (Ln = La, Pr, and Nd) nanoparticles were fabricated by the sol-gel method, and much work was devoted to studying their structural and magnetic properties [20]. Another conclusion is that all the samples showed considerable MCE in a wide temperature range. Some research results indicated that the physical properties of thin films are different from their bulk counterparts [21, 22]. A report suggested that MCE was successively reached at a low magnetic field [23], and the study also paid attention to the MCE of Ni-Mn-Ga thin films. Further, some other studies focused on the LaFe_13−x Si_x compounds [24], Heulser alloys Ni-Mn-X with (X = Ga, Sn,...) [25], [Gd/Tb] multilayers [26], and so on. The experimental research of BiFeO₃-based ceramics shows that the external magnetic and electric field have a significant influence on the magnetic properties of the system [27, 28]. In addition, the atom doping concentration can modify the MCE of the BiFeO₃-based system [29].

On the other hand, more efforts have been devoted to the theoretical investigations of MCE for nanomaterials, such as nanoparticles, nano-compounds, and thin film by various methods. A study was implemented by Hu et al to study the MCE in anisotropic zinc ferrite nanoparticles based on Monte Carlo (MC) simulation [30]. The magnetic entropy change and adiabatic temperature change were predicted in detail, and the emphasis was on the effect of anisotropy on these properties. In order to foster favorable materials for magnetic refrigeration devices, a study was carried out on the magnetic properties and MCE of the MnFe₄Si₃ system based on the MC method [31]. The structure was composed of a hexagonal lattice with Fe (S = 3/2) and Mn (S = 3/2). The magnetic entropy change, adiabatic temperature change, RCP, as well as magnetization, were obtained for different parameters and discussed in detail. In addition, the second-order phase transition was predicted at the Curie temperature. The results in this study provided good proof that MnFe₄Si₃ can be a suitable material for magnetic refrigeration. At present, many researchers are dedicated to the investigation of MCE for thin films by various methods. Giant magnetic entropy change is mainly rooted in the magnetic state change, which was concluded from a study of magnetic properties and MCE of multifunctional α-FeRh thin films [32]. MC simulation, as one of the most conventional and successful methods, has also been extended to study the MCE of thin films. An investigation of MCE for La_1−x Ca_x MnO₃ thin films based on this method revealed that the external magnetic field has an important impact on MCE for such thin films [33]. The surface effect on the magnetocaloric properties of perovskite ferromagnetic thin films was studied by Elyacoubi et al based on MC simulation [34]. Some results were discussed in detail, such as the maximal entropy change and that RCP increases when enlarging the external magnetic field, while RCP decreases when the value of the surface exchange coupling increases. In the field of theoretical calculation of magnetic nanomaterials, the mixed-spin Ising model is a well-known model which can help in understanding the magnetic nature of various magnetic nano-systems and could be taken into account to investigate the MCE. Albayrak studied the entropy and isothermal entropy change of the mixed-spin(1,1/2) Ising model [35]. The influence of the crystal field and external magnetic field on the variation of the peak value of isothermal entropy change was discussed emphatically. Apart from the research mentioned above, there are many other investigations on the MCE of thin films, for instance AlCMn₃ [36], NdMnO₃ [37], Ho₃Pd₂[38], Sr₂CrMoO₆ [39], Sr₂CrIrO₆ [40], and SmFe_1−x Mn_x O₃ [41]. All these studies reveal that exchange coupling and an applied external field play a significant role in the variation of the MCE, and this is exactly what our present research focuses on.

In our previous works, the MC simulation and the Ising mixed-spin model were significantly applied to investigate the magnetic, electric, as well as thermal dynamic, properties of various nanosystems [42–48]. We have ever devoted ourselves to the study of the dielectric properties of the BiFeO₃/FM bilayer system [49]. It should be mentioned that there have been many experimental [50, 51] and theoretical [52, 53] works which have studied the MCE behavior of multiferroic materials, and a study has released that BiFeO₃ has great potential applications because of its higher Curie temperature (T_C = 1100 K) and Nell temperature (T_N = 650 K) [54, 55]. However, few theoretical works have been conducted on the MCE in an AFM/FM BiFeO₃/FM bilayer system. This further stimulated our interest to conduct such a study to discuss, in depth, the effect of exchange coupling, the external magnetic field on MCE, and the magnetic properties of the BiFeO₃/FM bilayer system. The rest of this paper is arranged as follows: in the next section, the model and the method applied to this study are presented. In the third section, we show the detailed results of the simulation and provide some discussion. The last section includes a brief summary.

An AFM/FM BiFeO₃/Co bilayer with mixed-spin (5/2, 3/2) consisting of 2N² × L spins on a three-dimensional honeycomb lattice, as plotted in figure 1, is constructed. N and L are defined as the number of spins in each layer and the thickness of the bilayer system, respectively. For such a bilayer system, we adopt L = 2.

Zoom In Zoom Out Reset image size

The Hamiltonian of the system is considered as follows:

$$\begin{eqnarray}\begin{array}{rcl}H & = & -{J}_{{ab}}\sum _{\langle i,j\rangle }{\sigma }_{{ia}}^{z}{\sigma }_{{jb}}^{z}-{J}_{{cd}}\sum _{\langle m,n\rangle }{S}_{{mc}}^{z}{S}_{{nd}}^{z}-{J}_{z}\\ & & \times \sum _{\langle i(j),m(n)\rangle }{\sigma }_{{ia}({jb})}^{z}{S}_{{mc}({nd})}^{z}\\ & & -h(\sum _{i(j)}{\sigma }_{{ia}({jb})}^{z}+\sum _{m(n)}{S}_{{mc}({nd})}^{z}),\end{array}\end{eqnarray} \tag{ 1 }$$

where the first item J_ab (<0) represents the intralayer exchange coupling between sublattices a and b in the upper layer, and J_cd (>0) refers to the intralayer exchange coupling between sublattices c and d in the lower layer. J_z denotes the interlayer exchange coupling which is positive between sublattices a and c, but negative between sublattices b and d. In the whole simulation process, J_ab is considered to be the unit of energy and temperature, thus ∣J_ab ∣ takes the value of 1. 〈i, j〉, 〈m, n〉, and 〈i, m〉 stand for the summation over all the first nearest neighbor sites at sublattices a, b, sublattices c, d, and between sublattices a(b) and c(d), respectively. ${\sigma }_{{ia}}^{z}$,${\sigma }_{{jb}}^{z}$ are the z-component spins of spin-5/2 with the values ±5/2, ±3/2, ±1/2 and ${S}_{{mc}}^{z}$, ${S}_{{nd}}^{z}$ are the z-component spins of spin-3/2 with the values ±3/2, ±1/2. The last item h is the applied external magnetic field.

The method used for this simulation to calculate the BiFeO₃/Co bilayer is MC simulation based upon the standard Metropolis algorithm [56]. For the present bilayer, the periodic boundary condition is applied in the xy plane, and the free boundary condition is considered in the z-direction. In order to determine the finite-size effect caused by the number of each layer 2 × N² for such a bilayer system, we present in figure 2 the temperature dependence of susceptibility χ when 2 × N² is changed from 2 × 5² to 2 × 50² with fixed ∣J_z ∣/∣J_ab ∣ = 0.8, J_cd /∣J_ab ∣ = 1.2, T/∣J_ab ∣ = 0.05. From the figure, we can see every χ curve has a peak that corresponds to a temperature defined as critical temperature (T_C). It is an important parameter for better understanding magnetic phase transition from order to disorder for magnetic nanostructures [28, 34, 35]. It is found that the χ curves intersect each other at the same T_C when 2 × N² is larger than 2 × 20², whereas the T_C appears earlier for 2 × N² smaller than 2 × 20², suggesting that the system can exhibit an obvious finite-size effect for small size. Consequently, in order to save computation time, we choose 2 × N² = 2 × 20² for the following simulation. In order to balance the system, the first 2 × 10⁵ MC steps (MCS) among the total was abandoned in our simulation and we used the rest 3 × 10⁵ MCS to calculate the average values of the thermodynamic quantities as follows:

Zoom In Zoom Out Reset image size

The magnetization of sublattices a, b, c, and d are calculated as:

$$\begin{eqnarray}&&{M}_{a}=\displaystyle \frac{1}{{N}^{2}}\sum _{i}{\sigma }_{{ia}}^{z},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{M}_{b}=\displaystyle \frac{1}{{N}^{2}}\sum _{j}{\sigma }_{{jb}}^{z},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{M}_{c}=\displaystyle \frac{1}{{N}^{2}}\sum _{m}{S}_{{mc}}^{z},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{M}_{d}=\displaystyle \frac{1}{{N}^{2}}\sum _{n}{S}_{{nd}}^{z}.\end{eqnarray} \tag{ 5 }$$

Thus, the total magnetization of the system is shown as:

$$\begin{eqnarray}&&M=\displaystyle \frac{1}{4}({M}_{a}+{M}_{b}+{M}_{c}+{M}_{d}).\end{eqnarray} \tag{ 6 }$$

The total susceptibility of the bilayer system is:

$$\begin{eqnarray}&&\chi =\beta (\langle {M}^{2}\rangle -\langle M{\rangle }^{2}),\end{eqnarray} \tag{ 7 }$$

where $\beta =\tfrac{1}{{k}_{{\rm{B}}}T}$, k_B is the Boltzmann constant and T represents the absolute temperature. In this study, k_B is assigned a value of 1 just to simplify the calculation procedure.

The characterization of the MCE can be described by the following formulas:

The magnetic entropy per site is obtained by:

$$\begin{eqnarray}&&{S}_{m}(T,h)={\int }_{0}^{T}\displaystyle \frac{{C}_{m}(T,h)}{T}{\rm{d}}T.\end{eqnarray} \tag{ 8 }$$

The magnetic entropy change is presented by:

$$\begin{eqnarray}&&{\rm{\Delta }}{S}_{m}(T,h)={\int }_{0}^{h}{\left(\displaystyle \frac{\partial M}{\partial T}\right)}_{{h}_{i}}{\rm{d}}{h}_{i},\end{eqnarray} \tag{ 9 }$$

where ${\left(\tfrac{\partial M}{\partial T}\right)}_{{h}_{i}}$ denotes themal magnetization at a fixed magnetic field h_i .

The adiabatic temperature change is deduced by:

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{{ad}}(T,h)=-T\displaystyle \frac{{\rm{\Delta }}{S}_{m}(T,h)}{C(T,h)}.\end{eqnarray} \tag{ 10 }$$

The C appearing in the above equation is the representation of the total magnetic specific heat. The lattice term, electronic term, and magnetic term are taken into account. On the strength of the Debye model [57], the approximated value of it can be described as:

$$\begin{eqnarray}&&C=7.378T+0.115{T}^{3}+{C}_{m}.\end{eqnarray} \tag{ 11 }$$

Finally, the RCP of the system is given by:

$$\begin{eqnarray}&&\mathrm{RCP}={\int }_{{T}_{1}}^{{T}_{2}}{\rm{\Delta }}{S}_{m}(T){\rm{d}}T,\end{eqnarray} \tag{ 12 }$$

T₁ and T₂ are the cold and hot temperatures corresponding to both ends of the 1/2-maximum value of ${\rm{\Delta }}{S}_{m}^{\max }$, respectively.

In this section, some significant results of the effects of multifarious physical parameters, for instance, exchange couplings (J_cd , ∣J_z ∣) and external magnetic field h on the magnetic properties, the MCE as well as the hysteresis loops of a mixed-spin (5/2, 3/2) Ising bilayer, are presented in detail and discussed.

3.1. Magnetization and susceptibility

Figure 3 shows the magnetization (M), susceptibility (χ), and the critical temperature (T_C) of the bilayer system under the influence of the exchange interaction J_cd with other fixed parameters ∣J_z ∣/∣J_ab ∣ = 0.8, h/∣J_ab ∣ = 0.5. Figure 3(a) shows the variation of magnetization with J_cd . Just one saturation is observed, so a conclusion can be drawn that M is not sensitive to the variation of interlayer exchange interaction J_cd . The M curve holds a constant when J_cd /∣J_ab ∣ = 2.0 for T is lower (0 ∼ 1), while with the decreasing of the value of J_cd /∣J_ab ∣, this phenomenon gradually becomes inconspicuous. In the perspective of theory, the magnetization is presented under the interaction of exchange coupling, temperature, and external magnetic field. A larger exchange coupling value is more advantageous in competition with temperature, which causes the spins in the sublattice to align in its direction, thus showing a stable magnetization externally. However, when the exchange coupling decreases or the temperature increases, the spins have more chance to flit, so the magnetization changes for the response. It gradually decreases with the increase in temperature, and at the same temperature, it decreases with the decrease of J_cd /∣J_ab ∣. A similar phenomenon can also be seen in [38, 58]. The total susceptibility of the system varying with J_cd is displayed in figure 3(b). The six χ curves have the same shape. As the temperature rises, χ first increases and then gradually decreases after reaching a peak. The increase of J_cd /∣J_ab ∣ can cause the position of the peak to shift to a higher temperature region. Physically, with the increase of J_cd /∣J_ab ∣, more thermal energy is needed to overcome the forces from exchange coupling to break the magnetic order, thus T_C increases correspondingly. For a deeper understanding of the variation of T_C with J_cd , the T_C curve is plotted in the (T_C − J_cd /∣J_ab ∣) plane. In figure 3(c), a smooth curve can be seen. With the increase of J_cd /∣J_ab ∣, the curve gradually extends upward. In addition, the T_C curve splits the entire phase plane into two parts, with the ferromagnetic phase below T_C and the paramagnetic phase above T_C. A similar phenomenon can be seen in [47, 48].

Zoom In Zoom Out Reset image size

Similarly, as shown in figure 4, the influence of intralayer exchange coupling ∣J_z ∣/∣J_ab ∣ on M, χ, T_C when J_cd /∣J_ab ∣ =1.2, h/∣J_ab ∣ = 0.5 is studied. The M curve plotted in figure 4(a) is slightly different from that in figure 3(a). There is just one saturation as well, but as the temperature rises, the six M curves finally coincide with each other, which is different from the phenomenon in figure 3(a). Figure 4(b) shows the change of T_C with ∣J_z ∣. Obviously, the peak of the curves shifts to the right as ∣J_z ∣ increases. The curve is shown in figure 4(c), depicted in the (T_C − ∣J_z ∣/∣J_ab ∣) plane. It is clear that as ∣J_z ∣/∣J_ab ∣ increases, the T_C goes up in response. The phase plane is also divided into two parts by T_C—the lower and upper parts correspond to the magnetic phase and paramagnetic phase, respectively. Taking interlayer exchange coupling J_cd and intralayer exchange coupling ∣J_z ∣ into consideration, we can conclude that both have a similar impact on χ and T_C. Similarly, these results are in good accordance with that in [46, 48, 58, 59].

Zoom In Zoom Out Reset image size

3.2. Magnetocaloric effect

Figure 5 shows the magnetic entropy change (−ΔS) and adiabatic temperature change (ΔT_ad) of the bilayer system for different values of J_cd /∣J_ab ∣ for fixed values of Δh/∣J_ab ∣ = 0.5, ∣J_z ∣/∣J_ab ∣ = 0.8. From figure 5(a), a phenomenon can be seen whereby −ΔS initially increases with the increase of T/∣J_ab ∣, and then decreases with the continuous increase of T/∣J_ab ∣ after reaching the maximum value. We have to mention that the temperature corresponding to the maximum is near the critical temperature. This is consistent with that shown in [14, 37, 40], which present that the maximum of the −ΔS appears in the vicinity of the critical temperature. Thus, the prominent increase in the variation of −ΔS results from the heat that the system releases during the phase transition. Turning our attention to figure 5(b), similar curves to those shown in figure 4(a) are seen again. The biggest difference between them is the order of magnitude. When J_cd /∣J_ab ∣ takes the value of 0.2, the −ΔS is observed to be 0.033 75, while the ΔT_ad is observed to be 0.004 25. Moreover, the ΔT_ad curve tends to decline with the increase of T/∣J_ab ∣ more obviously than the −ΔS curve.

Zoom In Zoom Out Reset image size

The −ΔS and ΔT_ad of the bilayer system under the effect of ∣J_z ∣/∣J_ab ∣ for fixed J_cd /∣J_ab ∣ = 1.2, Δh/∣J_ab ∣ = 0.5 is shown in figure 6. In figure 6(a), each curve has a maximum, and the temperature corresponding to it is in the vicinity of T_C. The location of the maximum shifts to the high temperature region as ∣J_z ∣ increases. The value of it can be seen as 0.026 when ∣J_z ∣/∣J_ab ∣ = 0.2. Figure 6(b) presents the maximum moves to the high temperature region with the increase of ∣J_z ∣/∣J_ab ∣, and it can reach 0.0028 for ∣J_z ∣/∣J_ab ∣ = 0.2.

Zoom In Zoom Out Reset image size

All the results regarding the properties of the bilayer system are simulated in the presence of an external magnetic field, so it is of great necessity to investigate the dependence of −ΔS and ΔT_ad on the Δh, which is displayed in figure 7. Figure 7(a) gives the −ΔS as a function of Δh and T. The other parameters are selected as J_cd /∣J_ab ∣ = 1.2, ∣J_z ∣/∣J_ab ∣ = 0.8, respectively, and this figure shows us an intuitive phenomenon that when T changes from 6.0 to 9.8, there exists the maximum and it increases as Δh changes from 3.25 to 4.0. Figure 7(b) presents the ΔT_ad as a function of Δh and T. The maximum of the curves can be seen when T changes from 4.5 to 8.8, and it increases with the enlargement of Δh. Some comparisons may be implemented between these results and other studies, which investigated the influence of an external field on the −ΔS and ΔT_ad of the system [14, 31, 34, 35, 37, 38].

Zoom In Zoom Out Reset image size

The relative cooling power (RCP) is a significant physical parameter to elucidate the character of the MCE. As given in figure 8, the Δh dependence of RCP corresponds to different values of exchange couplings—figure 8(a) for J_cd /∣J_ab ∣ = 0.6, 1.2, and figure 8(b) for ∣J_z ∣/∣J_ab ∣ = 0.2, 0.8. It can be seen from figure 8(a) that RCP increases monotonously and smoothly as the field Δh increases, and can reach values of 2.00 and 1.75 for J_cd /∣J_ab ∣ = 0.6 and 1.2, respectively. Similar results can be seen in [31, 36, 38], and RCP decreases gradually with the increase of J_cd /∣J_ab ∣ for a certain Δh. A similar phenomenon in figure 8(a) can also be seen in figure 8(b)—that is, RCP goes smoothly with the increase of Δh, and the large ∣J_z ∣/∣J_ab ∣ corresponds to the small RCP. Therefore, a conclusion can be drawn that weakening the exchange coupling and heightening the external magnetic field can contribute to magnifying the RCP of the system and make the material a better candidate for a magnetic refrigerant. This conclusion is in accordance with [34, 37, 40].

Zoom In Zoom Out Reset image size

3.3. Hysteresis loops

In this subsection, we discuss the hysteresis loops of the system under the influence of exchange couplings, and the temperature dependence of the hysteresis loop is also presented. Figure 9 presents the magnetization as a function of the external magnetic field h for different values of J_cd /∣J_ab ∣ with ∣J_z ∣/∣J_ab ∣ = 0.8, T = 0.05. Only one loop can be seen and the area of the loops gets wider with J_cd /∣J_ab ∣ changing from 0.1 to 2.1. It is reasonable that the coercive field and the residual magnetization augment with the increase of exchange coupling. This action originates from the enlargement of the exchange interaction energy. An analysis of these results shows that the present AFM/FM BiFeO₃/Co can exhibit some properties similar to those found in ferromagnetic materials [60, 61].

Zoom In Zoom Out Reset image size

Analogically, figure 10 displays the magnetization as a function of the external magnetic field h for various values of ∣J_z ∣/∣J_ab ∣ with J_cd /∣J_ab ∣ = 1.2, T = 0.05. Similarly, the increase of ∣J_z ∣/∣J_ab ∣ is not advantageous for the number of hysteresis loops, but it is beneficial for the increase of the area for the loops.

Zoom In Zoom Out Reset image size

Finally, the impact of temperature on the hysteresis loop with the parameter J_cd /∣J_ab ∣ = 1.2, ∣J_z ∣/∣J_ab ∣ = 0.8 is shown in figure 11. From the figure, we can see that the remaining magnetization and the coercive field decrease gradually as T changes from a lower value to a larger; eventually, the loop vanishes at T = 6.5, corresponding to the critical temperature. Similar phenomena have also been found in the investigations of nanoislands [62, 63], graphene-like quantum dots [64], ladder-type graphene nanoribbons [65], a graphene monolayer [66] and triple-layer [67] as well as other magnetic nanosystems [68–72].

Zoom In Zoom Out Reset image size

The magnetic properties, magnetocaloric effect as well as hysteresis behaviors have been investigated and predicted for a mixed-spin (5/2, 3/2) AFM/FM Ising BiFeO₃/Co bilayer applying the Monte Carlo method. The results show that either exchange coupling or an external magnetic field has an important influence on the magnetic properties and the magnetocaloric effects. The presence of J_cd , ∣J_z ∣ and h causes the system to undergo a rough process to accomplish phase transition, while the decrease of J_cd , ∣J_z ∣, and the increase of h has a significant advantage in increasing the values of Δh, ΔT_ad, and RCP. Furthermore, the larger J_cd and ∣J_z ∣ are beneficial for the enlargement of the area for the hysteresis loops. These results indicate that the BiFeO₃/Co bilayer has appropriate properties to be a suitable candidate for use as a magnetic refrigerant material. We hope the study that we have implemented can provide direction for experimental research and theoretical reference for other works. It can be remarked that, although much effort has gone into the magnetic and electric properties of the BiFeO₃/Co bilayer, there is still much more to do, such as dynamic magnetic properties as well as the effect of atom doping, which will be investigated in detail in the future.