Modeling temperature dependent exchange bias in systems with magnetoelectric chromia

The requirement of ultra-low power dissipation in spintronic devices triggers the desire to move to voltage controlled switching elements in lieu of current controlled ones. Magnetoelectric (ME) chromia has drawn much attention in this context due to its salient feature of switching the exchange bias (EB)¹ of the adjacent ferromagnet (FM) with applied voltage. The linear ME effect present in chromia manifests the EB by developing a net boundary magnetization^2,3 (BM). Phenomenologically, simultaneous application of electric (E) and magnetic field (H) below Néel temperature (T_N) favors one spin configuration over the opposite one due to the lift in degeneracy of the two configurations.^4,5 The applied field product EH determines the directionality of the BM and the EB of the adjacent FM.⁵ The manipulation of the directionality has been utilized to achieve electric field dependent switching of the EB in systems with bulk and thin film chromia.^6–8 

Evidently, quantitative understanding of the EB is essential to assess the possibility of potential application of systems built with chromia. Due to lack of the absolute measurement of the in situ BM property, the traditional Meiklejohn and Bean⁹ (MB) model of EB, with first-principles calculation, is inadequate to understand the experimentally observed low magnitude of the EB. Additionally, the non-monotonous temperature dependence of the experimentally observed EB is difficult to perceive through mere interfacial exchange between chromia and the adjacent FM. A fundamental approach to circumvent this issue of modeling the EB is to investigate domain formation in the AF and incorporate the domain wall (DW) energy in the EB calculation proposed by Mauri.¹⁰ Mauri’s DW model is reported to predict EB magnitude much less than that predicted by MB model.^10,11 Therefore, implementation of the model in estimating the EB in systems with ME chromia poses a plausible solution to explaining observed low EB in experimental systems.

In this work, Monte Carlo (MC) simulation has been deployed to study the temperature dependent EB motivated by Mauri’s DW model. Also, the capability of the ME effect in determining the EB at any temperature below T_N has been studied through estimation of the temperature dependent energy terms.

Cr₂O₃ has the structure of alpha-corundum with hexagonal closed packed arrangement of O^2- ions.¹² Cr³⁺ ions occupy the octahedral sites with 2/3 occupancy in a buckled arrangement between O^2- layers. Ideally, spins are arranged within the buckled Cr³⁺ layer according to the maximum order parameter (OP). FIG. 1(a) demonstrates the maximum positive OP arrangement (spin order: up–zero–down–up in the z direction). This configuration is said to have positive BM. The maximum negative OP refers to the opposite arrangement of spins. The choice between the two configurations depends on the sign of the applied ME energy, $αEH$⁠. Here, α refers to the ME susceptibility which is inherent to Cr₂O₃ and possesses interesting temperature dependence.¹³ In the absence of the ME energy, Cr³⁺ ions in a buckled layer take their equilibrium positions with nominal displacement. With the applied ME energy, the displacement between Cr³⁺ ions with opposite spins within the buckled layer changes and the surface develops an uncompensated spin arrangement with all spins pointing up (e.g. FIG. 1(a)) or down depending on the OP set by the sign of the ME energy. This is obviously an idealistic picture with no temperature effect taken into consideration. In situ switching of the uncompensated surface spins demonstrated experimentally¹⁴ suggests that the choice of the surface configuration presented in FIG. 1(a) is reasonable. FIG. 1(b) shows the cartoon model of formation of a 180° DW from the surface (z = 0) towards deep into the sample with OP reversed at the end of the DW. Details of this model will be discussed in Section III.

Continuous spin MC simulation has been performed throughout the work to calculate magnetic configuration and system energy.¹⁵ The Hamiltonian utilizing Heisenberg exchange and anisotropy is as following.

Here, J_i,j represents the exchange between spin i with angle $θi$ and spin j with angle $θj$⁠. The sum over index j extends to Cr atoms up to 5^th nearest neighbor with exchange parameters reported by Shi et al.¹⁶ The exchange parameters have been scaled to match T_N in bulk magnetization calculation. K_u denotes the uniaxial anisotropy per Cr³⁺ spin. Bulk anisotropy value (∼2x10⁵ ergs/cm³) has been assumed for K_u.¹⁷ A total of 50000 MC steps has been performed per Cr³⁺ spin with 5000 steps to equilibrate initially.

Change in magnetism at the surface from bulk is always interesting to analyze due to the change in coordination. Bulk sublattice magnetization has been calculated by applying periodic boundary conditions in all directions. Surface magnetization for up spin sites has been estimated by removing the down spins at the surface buckled Cr³⁺ layer as shown in FIG. 1(a). 24x24x36 Cr³⁺ sites have been considered for the MC simulation. Magnetization normalized to the zero temperature value for up spin sites in bulk and surface configuration is shown in FIG. 2. Magnetization vanishes at T_N = 307 K. The steep drop in the surface magnetization with temperature may be attributed to the loss by surface sites of half their neighbors, especially the first two nearest neighbors which have the most prominent exchange. Our MC simulation results can be compared with the mean-field calculation by Wysocki et al.¹⁸ It can be seen that our bulk sublattice magnetization at high temperature approaches T_N more sharply and low temperature magnetization is reduced. These are direct consequences of the Bloch T^3/2 law and collective spin oscillation near T_N, neither of which are captured by mean-field theory. It is also worth noting that the previous results¹⁸ neglect anisotropy.

The surface magnetization presented in FIG. 2 does not represent the BM useful for estimating the EB. Hence, it should not be used in determining the EB with MB model. For instance, using the calculated surface magnetization, a 30 meV AF exchange¹⁹ per Co atom attached to the surface would yield an EB of (Co 0.6 nm/Pd 1.0 nm)₃ according to MB model more than two orders of magnitude higher at low temperature than experiment.⁶ This intriguing issue of calculating EB in systems with Cr₂O₃ thus requires exploring another EB mechanism, e.g. DW model.

In this work, a 180° DW in Cr₂O₃ has been modeled by forcing opposite OP at the edges of the sample. The surface (z = 0) has been pinned with positive OP and the other end (z = t_sim) with negative OP as shown in FIG. 1(b). For the DW calculation, 24x24x130 Cr³⁺ sites have been considered. At each temperature, MC simulation has been initiated with constant change in spin angle from positive OP to negative OP (FIG. 1(b)). FIG. 3(a) depicts the change in cosine of the average angle of up spin sites, $<θ(z)>$ with distance from the interface. An inference drawn from the slopes of the lines in FIG. 3(a) is that the DW becomes narrower with increasing temperature. While the reduction in DW width with temperature is unusual in a FM, it has been experimentally observed in an AF near T_N.²⁰ It is worth mentioning that due to the fixed simulation thickness (t_sim) and perfectly symmetric edges of the sample (at z = 0, t_sim), the domains are shifted to the middle of the sample. In reality, we anticipate that domains are formed close to the interface.

The DW energy has been calculated from the energy difference between configurations with (opposite OP at the edges) and without (similar OP at the edges) a 180° DW. We have found that the invariance of the DW energy with simulation thickness requires that t_sim ≥ 30 nm as used in our computation. To maintain consistency with Mauri’s model, we will refer the DW energy, U_dom as half the energy per unit area of the 180° DW. The DW energy per unit area has been presented in FIG. 3(b). The waviness of the temperature dependent DW energy arises from the thermal fluctuations inherent to the simulation scheme. There is a noticeable peak at elevated temperature, which arises from an increase in exchange energy within the DW as previously reported in a FM.²¹ FIG. 3(b) also shows the exchange bias, H_E calculated using Mauri’s model H_E = U_dom/(M_FM t_FM). Here, M_FM and t_FM represent the saturation magnetization and the thickness of the adjacent FM respectively. For EB calculation, we have chosen Cr₂O₃/Pd 0.5 nm/(Co 0.6 nm/Pd 1.0 nm)₃ as the AF/FM system, as reported by He et al.⁶ The M_FM value for Co/Pd superlattice has been taken from the work by Engel et al.²² The peak value of H_E is remarkably smaller compared to that predicted by MB model (>10⁴ Oe) and the unusual presence of a peak near T_N qualitatively resembles the experiment.⁶ 

The EB estimation presented herein does not include the magnetoelectric effect explicitly. It has been assumed that the OP is already set by the ME energy at low temperature. To gain insight on how effective the ME effect is in setting the OP and hence determining the EB, comparative study of temperature dependent relevant energy terms is essential. According to Borisov et al., applied ME energy must overcome the Zeeman energy of the surface spins and the exchange energy, U_ex at the AF/FM interface to successfully switch the EB.⁵ In case of the DW model, U_ex will be replaced by the DW energy at the interface. Temperature dependent ME energy per unit area is given by, $UME=α(T)EHtAF$⁠. Here, applied field values are E = 16.2 kV/cm, H = 3 kOe, and Cr₂O₃ thickness, t_AF = 0.5 mm.⁶ The ME susceptibility, $α(T)$ for the 0.5 mm thick Cr₂O₃ has been taken from experimental measurement.²³ The Zeeman energy per unit area of the surface spins has been calculated as U_Zee = m_surf(T)H. The surface moment per unit area, m_surf(T) has been estimated from FIG. 2. FIG. 4 shows the calculated temperature dependent energies. According to our DW energy calculation, below 190 K, U_ME cannot exceed U_Zee and U_dom, making the effectiveness of ME energy in governing EB below this temperature dubious. Below 80 K, U_ME cannot even support a minimal surface Zeeman energy due to $α(T)$ becoming negative. Hence, we anticipate that ME effect cannot cause the EB in the entire temperature range below T_N. At low temperature, $α(T)$ droop makes ME effect ineffective in switching the EB with applied voltage. Low temperature EB might be arising due to a small residual uncompensated moment developed at Cr₂O₃ surface.

Magnetic properties of magnetoelectric chromia have been analyzed through MC simulation. Magnetization with configurational change and DW formation have been investigated. Calculated DW energy has been utilized to estimate the EB according to Mauri’s model. The calculated EB has been found to qualitatively resemble the experimentally observed behavior at high temperature and its magnitude is much smaller than the MB model. Hence, this work demonstrates the significance of analyzing domain formation in the chromia film for understanding the EB. In addition to that, important insight on the effectiveness of the ME effect in governing the EB at low temperature has been presented through energy arguments. Our analysis shows that the temperature dependent magnetoelectric susceptibility puts a limit on the capability of the applied ME energy in setting the exchange bias in the low temperature regime.

The authors would like to thank C. Binek for valuable information pertinent to the work. This work was supported by C-SPIN, one of six centers of STARnet, a Semiconductor Research Corporation program sponsored by MARCO and DARPA.