Electric-field manipulation of magnetization in an insulating dilute ferromagnet through piezoelectromagnetic coupling

Repetitive switching of magnetization in ferromagnetic materials is a key process in magnetic recording and information storage. Traditionally, this energy-intensive process is achieved by applying an external magnetic field or spin-polarized currents. An emerging alternative, however, is the use of an electric field to control magnetic anisotropy, which could offer a more energy-efficient solution compared to conventional methods. Different approaches have been proposed in the past to observe the magnetoelectric (ME) effect, in which coupling between magnetism and electricity is present, such as the gate-controlled hole concentration and the Curie temperature in dilute ferromagnetic semiconductors^1,2,3,4, voltage induced structural and magnetic phase transition^5,6, piezoelectric based modification of magnetism^7,8,9 or magnetoelectric effect in multiferroics^10,11. The ME effect has been observed in ferromagnetic⁷, antiferromagnetic^12,13, superparamagnetic⁹, and paramagnetic⁸ materials.

A first step to accomplish the repetitive switching of magnetization is the demonstration of the control of the key magnetic properties, such as the coercive field, saturation magnetization, or remanent magnetization via electric fields^3,14,15. A strong coupling between an electric field and magnetism was shown to exist and theoretically quantified in piezoelectric wurtzite (wz) Ga_1−xMn_xN in the paramagnetic phase, that is for Mn concentrations x ≲ 2.5% and temperatures T ≥ 2 K (ref. ⁸). Those samples were grown by metalorganic vapour-phase epitaxy (MOVPE)¹⁶. Now aided by a carefully elaborated molecular beam epitaxy (MBE) technique, we incorporate substitutionally up to 10% of Mn ions into GaN host^17,18. In those semi-insulating layers, Mn³⁺ ions are coupled by a short-range ferromagnetic (FM) superexchange^16,19,20,21, and FM ordering is realized in a percolation-like fashion^22,23,24.

Our combined experimental and theoretical studies of the ME effect in the ferromagnetic (Ga,Mn)N address two classes of questions. First, we demonstrate that the sign and magnitude of magnetization generated by an electric field, M_E, which we find by superconducting quantum interference device (SQUID) magnetometry, result from a change of uniaxial magnetic anisotropy driven by the inverse piezoelectric effect. Our conclusion is supported by two complementary theoretical approaches: (i) a macrospin model, in which changes of magnetic anisotropy are described in terms of magnetizing or demagnetizing anisotropy fields, depending on the direction of the external magnetic field with respect to the easy axis; (ii) an atomistic Landau-Lifshitz-Gilbert approach for interacting classical spins^25,26,27 with input parameters obtained by fitting our magnetization M(H) data. Our results, on the one hand imply that the ME effect can be understood and designed in well characterized magnetic systems. On the other, we find that a short time scale inherent to the LLG approach leads to an overestimation in the theoretical magnitudes of the remanent magnetization and magnetoelectric effect in dilute ferromagnets in the magnetic fields H → 0.

Second, the percolating ferromagnetism means the presence of an infinite FM cluster at T ≤ T_C. This cluster encompasses all spins at T = 0 but only about 20% of spins at T approaching T_C (ref. ²²). Thus, at any non-zero temperature, a large portion of the Mn ions resides in finite clusters of different sizes. Such a magnetic system is expected to exhibit a broad spectrum of relaxation times and, in the hysteretic regime, constitutes a glass far from thermal equilibrium. Thus, a possibility of perturbing ferromagnetic (Ga,Mn)N by an electric field offers a unique opportunity to test theoretical predictions concerning dynamics of driven non-equilibrium glass systems^28,29,30.

As we show, our data provide a strong support for persisting theoretical proposals arguing that a non-thermal and non-dissipative drive, bringing a glass beyond the linear response regime, can be described by enlarged effective temperature³⁰. Interestingly, this interpretation appears explaining surprising magnetization behaviour observed in dilute ferromagnetic semiconductor (Bi,Sb,Cr)₂Te₃ in the regime of the quantum anomalous Hall (QAH) effect³¹. Tapping of that QHE ferromagnet by gate voltage modulation resulted, according to local SQUID measurements, in magnetization shift towards the direction of the external magnetic field in the hysteretic regime³¹, as reported here for (Ga,Mn)N. As short-range ferromagnetic superexchange dominates also in topological ferromagnetic insulators³², we expect similar glassy-like relaxation towards thermal equilibrium under a time-dependent drive in both systems.

It is worth noting, however, that our conclusions should not be extended towards experiments, in which magnetization switching in dilute ferromagnetic semiconductors was revealed by resistance measurements^33,34, as an important additional ingredient there is Joule heating that may lead to surprising bistable behaviours³³. Still another physics may show up in mesoscopic QAH samples, in which resistance switching was found also in magnetic fields far above the coercivity field^35,36. In such submicron wires, a highly probable process of backscattering between edge states is controlled by (i) orientation of local magnetization and (ii) spatial arrangement of localized in-gap carriers that are frozen far from thermal equilibrium under the presence of the Hall electric field. Resistance switching dynamics was found to accelerate with the magnetic field in the ferromagnetic QAH case³⁶, as found earlier in mesoscopic antiferromagnetic spin-glasses³⁷. It is to be seen on whether the AC source-drain voltage in QAH nanostructures³⁶ just allowed tracing magnetization dynamics or rather had an effect on the dynamics and effective temperature of the electron glass.

In our work we demonstrate a clear control of magnetic anisotropy of piezoelectric and ferromagnetic (Ga,Mn)N layers by application of an electric field. We observe a reduction of the width of a hysteresis curve from 120 Oe to about 10 Oe as the magnitude of the applied AC voltage V_AC increases from 0.2 V to 2 V and a non-reversible switching of magnetization toward the equilibrium direction for magnetic fields close to the coercive field, occurring under a stepped change of the electric field either elongating or shortening the lattice parameter c.

Samples

Two samples with different thicknesses of (Ga,Mn)N layer d_(Ga,Mn)N and similar concentrations of Mn ions x are investigated here, namely, d_(Ga,Mn)N = 10 nm, x = 5.6% (Sample A) and d_(Ga,Mn)N = 160 nm, x = 6% (Sample B). They were grown along the c axis on a thick GaN buffer followed by an n⁺-GaN:Si layer (n ~ 3 × 10¹⁸ cm⁻³), which serves as a back-contact contact layer, and characterized according to protocols elaborated previously^17,18. The experimental findings^19,20 and the theoretical modeling³⁸ agrees that for x < 10% Curie temperature T_C of Ga_1−xMn_xN obeys the power law T_C  ∝ x^2.2, which well approximates the percolation exponential formula for small x values. The layers studied here do fall into this dependency exhibiting 7 ≲ T_C ≲ 8 K for x ≃ 6%.

In Fig. 1 we compare M(H) data of Sample A acquired at T = 2 K for both in-plane (H⊥c) and out-of-plane (H∥c) configuration for which the demagnetization field was subtracted from the external field value. We see that at that temperature T ≈ T_C/4 and the time scale of the experiment, magnetic hysteresis involves only a part of spins, substantiating expectations of the percolation model^22,23,24.

a, b Full points: experimental magnetization M of Ga_1−xMn_xN with x = 5.6% measured at T = 2 K for two perpendicular magnetic field orientations, in-plane (H⊥c) and out-of-plane (H∥c). Shape anisotropy (demagnetization) effects are removed from experimental data. Inset in (b), Close-up on the weak magnetic field regime (H in Oe). Lines in (a) and (b): magnetization calculated for magnetic field sweep direction shown by arrow using atomistic Landau-Lifshitz-Gilbert equations with Langevin dynamics for 9464 classical spins with single-ion magnetocrystalline anisotropy and Heisenberg ferromagnetic interactions between spins (up to 18th neighbours), as described in the theoretical subsection.

Full size image

The data indicate also that magnetocrystalline anisotropy is close to zero for the Mn content and strain in question. In this case, shape anisotropy makes that the easy axis is in-plane in ferromagnetic thin films, as found previously^18,20 and here. However, in samples with smaller x values, uniaxial magnetocrystalline anisotropy is clearly non-zero and its sign corresponds to the in-plane easy axis^8,39. Those data imply a variation with x of the parameter (xi =c/a-sqrt{8/3}), where c and a are lattice parameters of the wz unit cell. The parameter ξ governs the trigonal deformation and the corresponding uniaxial magnetic anisotropy in magnetically doped wz crystals like Mn-doped GaN^8,39,40 or Co-doped ZnO⁴¹.

To substantiate those findings, we show in Fig. 2a ξ values obtained by x-ray diffraction for a series of Ga_1−xMn_xN samples deposited by MOVPE and MBE on a GaN buffer clamping the a value. We also present results of our first-principles calculations with the structure optimization (see “Methods”) by plotting the computed values of ξ and the total energy difference Δε for (Ga,Mn)N with a Mn magnetic moment lying in a plane (along the a-axis) and a Mn magnetic moment aligned with the c-axis as a function of x. As seen in Fig. 2a, the energy difference changes the sign from negative to positive for the Mn concentration of about 2.5%, which correlates but is slightly lower than the computed value corresponding to the sign change of ξ. Those experimental and theoretical data imply that as a function of Mn content, and in the absence of the shape anisotropy, the Mn magnetic moment will have a tendency to reorient from the in-plane to perpendicular direction. Our two samples are close to that reorientation transition with ξ = 0.0029(3) (Sample A) and  − 0.0005(3) (Sample B), in agreement with magnetization data in Fig. 1.

a Left axis: dependence of the parameter (xi =c/a-sqrt{8/3}) on the Mn concentration x measured experimentally at 300 K (stars) essentially identical to evaluated values for 4.2 K (full squares) compared to determined by ab initio computations (open red circles). Right axis: total energy difference between the in plane (along the a-axis) and parallel to the c-axis configurations of the Mn magnetic moment in GaN (open blue squares). Labels indicate ξ values for Samples A, B and MOVPE film from ref. ⁸. b, c In Ga_1−xMn_xN, the uniaxial magnetic anisotropy originates primarily from the trigonal deformation of the tetrahedron formed by anion atoms next to Mn. The application of an external electric field E elongates lattice parameter c by the inverse piezoelectric effect and controls uniaxial magnetocrystalline anisotropy.

Full size image

Experimental

As presented in “Methods”, for magnetoelectric investigations Ga_1−xMn_xN layer resides between the n⁺-GaN:Si bottom contact and a top HfO₂ film covered by a metal gate. Such an arrangement allows for blocking gate currents resulting from a possible conductance via threading dislocations in the (Ga,Mn)N epilayer. In the case of the DC voltage applied to the structure, the electric field on the HfO₂ and (Ga,Mn)N films is distributed according to relative resistances of the two layers, so that a rather small electric field is expected in (Ga,Mn)N samples as, presumably, R_HfO2 ≫ R_(Ga,Mn)N. In contrast, for resistances and frequencies in question (f = 17 Hz), the AC voltage is distributed according to relative capacitances, meaning that the application of the AC voltage of V_AC = 1 V to the investigated structures (see “Methods”) corresponds to the AC electric field E_AC of about 0.19 MV/cm and 0.04 MV/cm to the (Ga,Mn)N layer of Sample A and B, respectively. Given the structure parameters and wiring layout, we do not expect temperature changes by capacitor charging currents at experimental temperature of 2 K.

The measurement system, depicted in Methods, utilizes a SQUID magnetometer and phase sensitive detection. An AC voltage V_AC is applied between the metal gate and the GaN:Si contact layer. The resulting AC magnetization signal of the sample is picked up by SQUID coils and detected by a lock-in amplifier. Alternatively, we look for magnetization alterations generated by stepped DC voltage changes. As shown in panels b and c of Fig. 2, in our piezoelectric structure with the clamped value of a by the thick buffer layer, an electric field E alters only the lattice parameter c, according to Δc/c = d₃₃E, where d₃₃ = 2.8 pm/V, as obtained for an epitaxial GaN film with the same orientation of the c vector as in our samples, and also with a clamped value of a by a thick GaN buffer layer⁴².

Experimental results and macrospin model

As mentioned in the introduction, our data reveal the existence of two distinct magnetoelectric phenomena. First of them, discussed in this subsection, corresponds to the regime, where the magnetoelectric response M_E is linear in the electric field E. Given a rather unique combination of ferromagnetism and piezoelectricity offered by single crystals of wurtzite (Ga,Mn)N, the data in the linear range makes possible to test quantitatively our understanding of magnetoelectricity in dilute ferromagnets. The second phenomenon, presented in the last subsection of Results, is specified by strong nonlinearities and irreversibilities of M_E. As we show, the data in that regime give access to the physics of driven glassy systems far from thermal equilibrium.

We recall that in our experiments the electric field is applied along the c axis of the wurtzite structure, and we measure magnetization changes in the direction of the applied magnetic field H. We have found that in high magnetic fields, M_E(H) is linear in E up to the highest electric field supported by our devices. In contrast, in the hysteresis regime, nonlinearties are already apparent in the electric fields lower by a factor of ten. Thus, being interested here in the linear regime, we discuss M_E values normalized to V_AC = 0.5 V (Sample A). Such a low value of V_AC has been employed in low magnetic fields H ≤ 1 kOe but V_AC has been risen up to 5 V for H ≥ 10 kOe.

According to the data displayed in Figs. 3 and 4 for Sample A at 2 K, the electric field decreases magnetization for the magnetic field H⊥c, but increases M for H∥c. Interestingly, comparing M_E(H) in Fig. 3 to M(H) in Fig. 1, we see that M_E is roughly proportional to ∂M/∂H in high magnetic fields but according to Fig. 4, M_E becomes proportional to M in weak magnetic fields. The presented data implies that the converse magneto-electric coupling coefficient, obtained from the relation α_C = 4πM_E/E, attains 7.3 Oe cm/MV for wz-Ga_0.94Mn_0.06N at 2 K and in 0.5 kOe, as E_(Ga,Mn)N = 0.095 MV/cm for 0.5 V applied to Sample A. The highest value of M_E we have detected is 0.23 emu/cm³ at 5 V, to be compared to the magnetization saturation value M_sat = 80 emu/cm³.

Areal density of the magnetic moment M_E/A normalized to the voltage V_AC = 0.5 V in Ga_1−xMn_xN layer with x = 5.6% at T = 2 K, i.e., in the ferromagnetic phase (blue solid points). The data are plotted as a function of the square root of the magnetic field to show simultaneously weak and strong field regions. Measurements are taken for the in-plane (H⊥c) and out-of-plane (H∣∣c) configurations. a, c Open red points: theoretical results of magnetoelectric signal M_E obtained from the macroscopic spin model (Eq. (1)) with one fitting parameter b = 0.04. b, d Open green stars: the magnetoelectric signal computed by using the atomistic Landau-Lifshitz-Gilbert (LLG) equations with Langevin dynamics.

Full size image

a Areal density of the in-plane magnetic moment amplitude M_E/A of Ga_1−xMn_xN layer with x = 5.6% induced by AC voltage V_AC = 1V at T = 2 K. The measurements are taken in the in-plane (H⊥c) and out-of-plane (H∥c) configuration (full circles and triangles, respectively). These results are compared to directly measured magnetizations  − M_⊥(H) and M_∥(H) (grey lines, relative units), measured for Sample B before processing. The virgin magnetization data obtained after cooling in H ≈ 0.15 Oe are shown by empty symbols. b Areal density of the in-plane magnetoelectric signal M_E/A at H = 0 as a function of temperature T (blue circles), compared also to experimental value M(T) (grey line, relative units). c Magnetoelectric response to V_AC = 1 V measured in the presence of the DC bias V_DC = ±28 V applied in the magnetic field beyond the hysteretic regime.

Full size image

We describe M_E(H) in this easy-plane piezoelectric ferromagnet in terms of an effective anisotropy magnetic field H_E affected by the electric field. If H_E ≪ H,

where j refers to the two explored configurations of the magnetic field; χ_jj = ∂M_j/∂H_j and, since uniaxial anisotropy is proportional to ξ(E),

where  ± corresponds to H∥c and H⊥c, respectively, and b is a dimensionless fitting parameter. The above equations define the form of the magnetoelectric coefficient⁴³α_C = M_E/E for the case under consideration.

We see in Fig. 3 that our macrospin model with experimental values of χ_jj(H_j) and b = 0.04 describes quite satisfactorily the field dependence for the two orientations of the magnetic field down to 4 kOe at V_AC = 0.5 V.

While in high magnetic fields M_E(H_j) ∝ χ_jj(H), in weak magnetic fields, M_E(H_j) becomes proportional to M_j(H_j), as shown in Fig. 4. This observation supports experimentally the expectation that α_C is non-zero only under time-reversal symmetry breaking⁴³, i.e., α_C is odd in H_j in the paramagnetic case, as observed previously⁸, or odd in M_j at H_j → 0, as found here for the ferromagnetic ordering. This means that for sufficiently small values of M_j, M_E should be linear in M_j, the relation fulfilled by our data displayed in Fig. 4a, b as a function of the magnetic field at 2 K and as a function of temperature at H = 0, respectively. Accordingly, as shown in Fig. 4b, the magneto-electric effect vanishes near T_C, about 7.5 K for Sample A. Quantitatively, in weak magnetic fields, M_E(H_j) = ±βM(H_j), where β = 2 × 10⁻³ for Sample A and V_AC = 1V.

As presented in Fig. 4c, we have also checked that the application of the DC bias of either polarization has no effect on M_E(H) generated by V_AC. This finding is consistent with the fact that the DC electric field is non-zero only across the HfO₂ film.

As far as the origin and magnitude of the parameters b and β are concerned, we note that, in general, crystalline magnetic anisotropy, which we perturb by the electric field via the inverse piezoelectric effect, originates from a single-ion contribution and anisotropy of spin-spin coupling. In magnetic compounds in question, the spin-orbit interaction within anions transmitting exchange between localized spins plays a dominant role^44,45, so that this anisotropy should be relatively small for the light nitrogen ions. This fact suggests that single-ion anisotropy, rather sizable for Mn³⁺ ions with the orbital momentum L = 2, dominates in (Ga,Mn)N. The signs of b and β are consistent with this expectation. However, though we know rather precisely the crystal field parameters of Mn³⁺ ions in GaN and their dependence on ξ (ref. ⁸), there are persistent difficulties in relating that information with parameters of the macrospin model⁴⁶, particularly in dilute magnetic systems⁴⁷ and in the percolating case. This fact makes a theoretical determination of b and β difficult within the macrospin model. We have, therefore, decided to adapt a stochastic Landau-Lifshitz-Gilbert (LLG) atomistic spin model^25,26,27 for the case under consideration.

Landau-Lifshitz-Gilbert atomistic spin model

The way that approach^26,27 is adopted to ferromagnetic (Ga,Mn)N is detailed in Methods. Our atomistic LLG model was also employed to describe magnetization of (Ga,Mn)N in the paramagnetic regime, where magnetic properties are determined by single Mn ions, nearest neighbour Mn pairs, and triangles, so that the outcome of the LLG model which treats the Mn spin classically, could be compared to the results of computations describing the spins quantummechanically⁴⁸.

The magnetic energy of Mn spins is described by a spin Hamiltonian ({{mathcal{H}}}) taking into account the Zeeman energy, Heisenberg spin-spin interactions between Mn spins located up to the 18th nearest neighbour positions, and the magnetocrystalline anisotropy energy, ({{mathcal{H}}}={{{mathcal{H}}}}_{{{rm{Z}}}}+{{{mathcal{H}}}}_{{{rm{Exch}}}}+{{{mathcal{H}}}}_{{{rm{Anizo}}}}), where ({{{mathcal{H}}}}_{{{rm{Anizo}}}}) includes the Jahn-Teller contribution and the trigonal anisotropy term taken as,

Here, K_Trig ∝ ξ is the anisotropy parameter describing the trigonal deformation that can be controlled by the electric field E, and S_i is a unit vector describing the local magnetic moment direction at the site i. The explicit forms of ({{{mathcal{H}}}}_{{{rm{JT}}}}) and ({{{mathcal{H}}}}_{{{rm{Exch}}}}) together with the values of relevant exchange integrals, are given in Methods.

The results on thermal equilibrium values of macroscopic magnetizations M_j(H_j) without and with an electric field E and for two directions j of the external magnetic field are obtained for a simulation box 21 × 21 × 18 nm³ containing 169,000 cation sites, 9464 (5.6%) of which encompass randomly distributed spin vectors. The site number is limited only by the available hardware employed to trace the time evolution towards thermal equilibrium values M_j(H_j) over an acceptable time frame. The simulations start in 50 kOe, in which Mn spins are aligned along the magnetic field, and next the field is decreased in steps specified in “Methods”.

Considering the lack of a clear correspondence between parameters of the quantum and classical spin Hamiltonians, our strategy involves two steps. First, by fitting the LLG model to the experimental data on M_j vs. H_j at E = 0 and for the two orientations of the external magnetic field, as presented in Fig. 1, we have evaluated K_JT = 0.85 meV and K_Trig = 0.1 meV appearing in system Hamiltonian. A positive sign of K_Trig is in agreement with the data on ξ obtained from x-ray diffraction (XRD) measurements, shown in Fig. 2a. However, it should be noted that the simulations point to larger remanence than observed experimentally. Moreover, even rescaling Curie temperature obtained from the numerical modeling for classical spins by a factor (S + 1)/S, we obtain ({T}_{{{rm{C}}}}^{{{rm{sim}}}}cong 4.5) K, which differs significantly from the experimental value T_C ≅ 7.5 K.

Second, we have performed computations for the determined value of K_JT and for K_Trig in the absence and presence of the electric field,

where ξ = 0.0029(3) for Sample A according to Fig. 2a. Here, in order to reduce uncertainties resulting from statistical fluctuations of computational results, we considered E corresponding to V_AC = 5 V, and reduce the numerical values of M_E(H_j) by a factor of ten to compare to experimental results depicted in Fig. 3 for V_AC = 0.5 V. As seen in Fig. 3, the LLG theory is successful only qualitatively, i.e., predicts properly the sign and shape of M_E(H_j). The obtained results bring us to the conclusion that the LLG approach, when applied to dilute ferromagnets, predicts too large magnitudes of remanence in the case of both M(H_j) and M_E(H_j). We assign this outcome to time scale of our experiment, which is sufficiently long to allow demagnetization of most Mn spins outside the percolation cluster. Such a demagnetization process does not occur in the necessary much shorter time frame of the LLG simulations, so that most of Mn spins contribute to the theoretical remanence, as shown by solid lines in Fig. 1.

Nonlinearities and irreversibilities

We now discuss surprising switching effects found in the hysteretic regime. The magnetoelectric signal M_E(H) in this range and at 2 K is plotted in Fig. 5a as a function of the in-plane magnetic field H⊥c for three different AC modulation voltages V_AC. As seen, the width of the coercive field H_c gets sizably reduced (from 120 Oe to about 15 Oe) by increasing V_AC to 5 V. The effect is better visualized in Fig. 5b, c, where we present changes of the coercive field with the applied voltage V_AC to the whole structure and with the corresponding electric field E_AC across the (Ga,Mn)N layer of the thickness d = 10 and 160 nm for Sample A and B, respectively. We see that the M_E appears highly nonlinear in E_AC in the hysteretic region. According to Fig. 5c, this nonlinearity is stronger in Sample B in which, judging by the magnitude of H_c, ferromagnetism is softer than in Sample A. Altogether, these results imply that taping of the system by sufficiently high V_AC brings the system towards thermal equilibrium.

a Areal density of the in-plane magnetic moment amplitude M_E/A generated by three different values of the AC voltage V_AC as a function of the in-plane magnetic field. b, c Corresponding coercivity field as a function of V_AC and the electric field E_AC across the (Ga,Mn)N layers of thicknesses 10 and 160 nm for Samples A (green circles) and B (red squares), respectively.

Full size image

While a time-independent DC bias has no effect in our structures (see Fig. 4c), we find that step changes of the DC voltage, if applied in the non-equilibrium range of the hysteresis, moves magnetization toward the direction of the external magnetic field. Our experiment begins by sweeping the external magnetic field to the non-equilibrium range of the hysteresis, and then we change V_DC in a step-like manner by ΔV_DC at constant H monitoring M_E at selected voltage values V_DC = 0, 3,…, 18 V or, for the negative sign of ΔV_DC, at V_DC = 0, −3, …, −18 V. There is at least 1 sec delay between consecutive voltage changes.

Results of such investigations are collected in Fig. 6 which shows that: (i) Magnetization changes generated by voltage alterations occur only in the non-equilibrium regime and those changes are irreversible, i.e. magnetization does not change under removal of voltage. (ii) According to results depicted in Fig. 6a, b, outside the metastable regime, the observed drift in M_E magnitudes for consecutive voltage changes is of the order of time-dependent fluctuations in the M_E values, shown in Fig. 6c. (iii) Voltage-induced shifts in M_E(V_DC) values are independent of the DC voltage polarization, and a sufficiently high magnitude of ∣V_DC∣ reverses the magnetization direction (Fig. 6a, b). (iv) As shown in Fig. 6d, the consecutive values of M_E(V_DC) are also virtually independent of the number of steps leading to a given V_DC, at least in the range of step heights 1 ≤ ΔV_DC ≤ 15 meV.

a, b At V_DC = 0 V, the magnetic field is swept from −450 to +50 Oe or to −50 Oe (blue circles, full sweep is shown). The V_DC is then changed in step-like manner by ΔV_DC = 10 mV or −10 meV, and the resulting magnetoelectric signal M_E is measured at successive points V_DC = 0, 3,…, 18 V (orange points or magenta points in a) or V_DC = 0, −3, …, −18 V (orange points or green points in b). Irreversible switching of magnetization is observed independently of the DC voltage polarization if for H < 0, M > 0 (i.e., M_E < 0 for H⊥c configuration), meaning that the system is initially out of thermal equilibrium. c Magnetoelectric response M_E as a function of time for ΔV_DC = 0. d, Magnetoelectric response M_E as a function of V_DC for various magnitudes of DC voltage steps ΔV_DC.

Full size image

It is remarkable that irreversible shifts of magnetization toward the thermal equilibrium value by voltage alterations, as we report here, were also detected at 0.25 K by local SQUID magnetometry in gated structures of another dilute ferromagnetic semiconductor (Bi,Sb,Cr)₂Te₃ (ref. ³¹). As in our case, the effect was independent of gate voltage polarization and the magnetization alteration scaled with the magnitudes of the gate voltage excursion. Also, time-independent gate voltage had little effect on magnetization and hystereses. Furthermore, spatially-resolved magnetization studies pointed to the presence of superparamagnetic regions with noticeable magnetization relaxation over the laboratory time scale.

To explain the origin of that universal behaviour, we first recall that despite magnetic dilution, a well-defined magnetic domain structure was found in (Ga,Mn)As, the fact explained quantitatively in terms of a long-range character of spin-spin interactions in that type of ferromagnets⁴⁹. In order to address this question in (Ga,Mn)N, we show in Fig. 4a experimental data on virgin magnetization obtained after cooling the sample down to 2 K under near-zero-field conditions (H ≈ 0.15 Oe). The measured form of m(H) corresponds neither to the domain-wall motion nor to the domain-wall pinning in a typical multidomain case for a ferromagnet. In the former case one would expect a rapid increase of m, nearly along the y-axis, in the latter the initial m(H) would be horizontal along the x-axis. Actually, the observed behaviour is consistent with or even supports the view that in the case of short-range interactions, driven here by superexchange, the virgin magnetic properties are dominated more by superparamagnetic regions present in the percolation picture than by a regular domain structure which is formed by dipole interactions in dilute ferromagnets with long-rage exchange couplings⁴⁹.

More specifically, we expect in both (Bi,Sb,Cr)₂Te₃ and (Ga,Mn)N, a co-existence, at T_C > T > 0, of an infinite ferromagnetic percolating clusters with superparamagnetic regions, characterized by a broad spectrum of blocking temperatures and relaxation times. With no doubts, the presence of long-range dipole interactions, responsible for the domain formation in materials like (Ga,Mn)As, contributes to the complexity of spin dynamics. Therefore, (Bi,Sb,Cr)₂Te₃ and (Ga,Mn)N in the metastable regime can be regarded as a glass that is far from thermal equilibrium. A question arises: how do such media respond to a time-dependent perturbation? It has been theoretically demonstrated that in driven glassy or granular systems, the degrees of freedom governing slow relaxation thermalize to an effective temperature, whose deviation from bath temperature scales with the drive strength^28,29,30. Considering a steady decay of the coercive force with temperature, we see that the data on dilute ferromagnetic semiconductors with a short range of spin-spin interactions provide the experimental corroboration of theoretical predictions concerning tapped glassy systems. However, neither in (Ga,Mn)N, nor in (Bi,Sb,Cr)₂Te₃ a full magnetization reversal has been attained by time-dependent voltage manipulations, at least in the explored time-domain much longer than the inverse precession frequency. This fact is in accord with our scenario that implies the presence of clusters in which a ferromagnetic spin arrangement is robust to weak perturbations.

To the best of our knowledge, the experiments have revealed novel aspects of driven dynamics, such as the fact that the trajectory toward the thermal equilibrium depends only on the accumulated magnitude of the drive changes but not on the number of those changes. It is worth noting that the standard Langevin dynamics does not capture the effects of time-dependent perturbations, so that the nonlinear and irreversible effects do not show up in our LLG simulations.

In our work, we have exploited the recent progress in the MBE growth of high-quality quantum structures incorporating ferromagnetic (Ga,Mn)N to examine much needed spintronic functionality, i.e., the manipulation of magnetization by an electric field. Our results demonstrate that (Ga,Mn)N with a sufficiently large Mn concentration is one of rather rare piezoelectric ferromagnetic homogeneous compounds, in which the magnetoelectric effect can be detected and, moreover, examined theoretically in a quantitative manner. In another perspective, the magnetization manipulation enlarges a pallet of wide-ranging nitride functionalities⁵⁰ that have made GaN and related systems the second most important semiconductor family next to Si. In particular, GaN:Mn serves already as a semiinsulating substrate⁵¹ in device epitaxy. The spectrum of (Ga,Mn)N functionalities has recently been enriched by the demonstration of efficient generation and detection of spin currents in Pt/(Ga,Mn)N structures up to 50 K (ref. ⁵²) using spin Hall magnetoresistance (SMR) effect⁵³. The same SMR could be then employed to probe the magnetic state of the insulating (Ga,Mn)N in the electrical way, paving the way towards more practical applications.

Our experimental results have revealed the presence of two distinct regimes showing a different pool of phenomena. The first of them concerns a wide range of magnetic and electric fields, in which the magnetoelectric response M_E is linear in the AC electric field E and reversible, i.e., the magnetization magnitude returns to its initial value after switching off the electric field. We have developed for that case (i) the macrospin approach and (ii) the atomistic theory with Landau-Lifshitz-Gilbert equations containing the external magnetic field, molecular fields generated by neighbouring Mn spins, and Gaussian magnetic noise—the Langevin dynamics. The macrospin model, containing one fitting parameter, describes M_E for both magnetic field configurations, H⊥c and H∥c over a wide field range, including the hysteretic regime. In the case of the LLG theory, isotropic ferromagnetic Mn-Mn exchange integrals come from the theory which describes T_C(x), whereas parameters of the single-ion anisotropy, from fitting M(H) for the two configurations. A general agreement between experimental and theoretical results on M_E(H) demonstrates that the magnetoelectric response results primarily from the influence of the inverse piezoelectric effect on the wurtzite c lattice parameter, which affects the trigonal uniaxial single-ion magnetocrystalline anisotropy. At the same time, we have identified a limited predictive power of the LLG simulations concerning the magnitudes of remanent magnetization and magnetoelectric effect in dilute ferromagnets at H → 0.

The other magnetoelectric phenomena discussed in our work appear in the metastable range of the low-field hystereses. Our results, and also those found for (Bi,Sb,Cr)₂Te₃³¹, let us to suggest that dilute ferromagnetic semiconductors with short-range spin-spin couplings show, in the metastable regime, characteristics similar to glasses far from thermal equilibrium. One of pertinent questions is the response of such systems to time-dependent perturbations specified by a time scale long compared to both the inverse Larmor frequency and the transverse spin relaxation time T₂, governed by non-scalar spin-spin interactions, but comparable to T₁ of superparamagnetic clusters incorporating most of localized spins at T > 0. To explain characteristics of irreversible magnetoelectric responses in (Ga,Mn)N and (Bi,Sb,Cr)₂Te₃, we refer to extensive numerical simulations of glasses^28,29,30, indicating that time-dependent perturbations of such systems lead to fast thermalization to an enhanced temperature magnitude facilitating a unidirectional drift of the magnetization value and direction toward the thermal equilibrium state. Conversely, magnetoelectric data for relatively simple and easy controllable systems, such as (Ga,Mn)N and (Bi,Sb,Cr)₂Te₃, can motivate further theoretical studies of non-equilibrium physics in disordered systems.

Material

The (Ga,Mn)N layers investigated here have been grown by plasma assisted molecular beam epitaxy (MBE) method in a Scienta-Omicron Pro-100 MBE at 620 ^∘C under near stoichiometric conditions following closely the previously elaborated protocol¹⁸. In order to reduce the number of threading dislocations and simultaneously provide the back electrical contact to the final capacitive device we use commercially available conductive low threading dislocation density (LTDD) n-type (0001) GaN:Si templates deposited on a c-plane sapphire, which contain, generally termed, an SiN_x mask. During the growth the surface quality of the layers has been investigated in situ by reflection high-energy electron diffraction (RHEED). The post growth [11(bar{2})0] RHEED pattern of Sample B collected after the growth (at a temperature below 200 ^oC) is shown in Fig. 7. The streaky and sharp RHEED pattern with clearly visible Kikuchi lines indicates that the (Ga,Mn)N surface is smooth and relatively free of oxides or other contaminants.

RHEED pattern along [11(bar{2})0] azimuth of Sample B observed after the growth and cooling the sample down to 200 ^oC.

Full size image

Two samples have been prepared for this study. Sample A of a low thickness of (Ga,Mn)N layer, d_(Ga,Mn)N = 10 nm and a thicker Sample B with d_(Ga,Mn)N = 160 nm. The precise values of d_(Ga,Mn)N are established upon high resolution x-ray characterization described below. Both are having a similar Mn concentration x, about 6%, established by direct magnetic measurements.

Structural characterization

The crystal structure, epitaxial relationship, in- and out-of-plane lattice parameters, and the strain condition of the samples are established by high resolution x-ray diffraction (HRXRD) using a Philips X’Pert MRD diffractometer operating on wavelength 1.5406 Å of Cu_Kα1 line. Its incident beam is formed by an x-ray mirror and a monochromator [four-bounce (220) Ge asymmetric cut]. Figure 8 collects the most relevant findings. Panel a demonstrates a 2θ/ω scan measured for Sample B over a wide range of detector angle changes, indicating that within a dynamic range of intensity exceeding four orders of magnitude no reflections specific to inclusions of other phases are seen. The reciprocal space map of the (bar{1}bar{1})24 reflection is given in panel b. Two main features can be recognized there: the high intensity peak located at Q_z = 4.827 Å⁻¹ – related to the GaN template, and a much weaker one at smaller reciprocal lattice units, at Q_z = 4.817 Å⁻¹, from the thinner (Ga,Mn)N layer. The map indicates a nearly perfect pseudomorphic growth, meaning that the (Ga,Mn)N layer adopts the same in-plane lattice parameter as that of the substrate (the GaN LTDD template).

a Diffraction pattern of Sample B in 2θ/ω mode. Three GaN reflections denoted with a star are a result of a multiple reflection and are seen here due to a specific orientation of the sample. The reflections indexed on the grey background are from the sapphire substrate. b The reciprocal space map of the asymmetric (bar{1}bar{1})24 reflection of Sample B. The high-intensity node is related to the GaN template layer. The low-intensity node reflects a practically fully strained (Ga,Mn)N layer, as can be inferred from the location of this node very close to the vertical line, which, together with the inclined dashed line define the triangle of relaxation. c, d HR diffraction pattern for the 0002 Bragg reflection for Sample A and B, respectively (the dark blue lines). The orange lines marks results of simulations upon which the values of the layers thicknesses are established (given in the legends).

Full size image

HRXRD patterns for the 0002 Bragg reflection are marked as the thicker blue line in Fig. 8c, d for Samples A and B, respectively. Clearly resolved x-ray interference fringes imply a high structural perfection of the layers and good quality of the interfaces. The peaks corresponding to the (Ga, Mn)N epitaxial layers are shifted to smaller angles with respect to that of the GaN buffer, a result of a larger perpendicular lattice parameter of (Ga,Mn)N. The orange lines mark results of simulations performed in the frame of the dynamical theory of x-ray diffraction using commercially available PANalytical EPITAXY software. In these computations we put the elastic stiffness constants in (Ga,Mn)N as those in bulk GaN and c = 5.1850(5) Å, and a = 3.1885(5) Å, as established in ref. ⁵⁴. The best fit to the experimental data is obtained for d_(Ga,Mn)N = 9.5(1) nm for Sample A and d_(Ga,Mn)N = 161(1) nm for Sample B, from which simplified values of 10 and 160 nm values are adopted to characterize both layers. We can parenthetically add that allowing for 5% of lattice relaxation (to account for a small horizontal shift of the (Ga,Mn)N node seen in Fig. 8b) does not call for a greater modification of d_(Ga,Mn)N than the modeling uncertainty specified above.

The XRD measurements are also used to establish the values of the trigonal distortion parameter (xi =c/a-sqrt{8/3}). Generally, the value of c is obtain from the symmetrical 0002 reflection, a parameter is calculated based on data from both 0002 one and asymmetrical −1–124 reflections. Bragg angle of (bar{1}bar{1})24 is usually determined from a relevant reciprocal space map, as the one presented in Fig. 8b. The bare results for ξ obtained from XRD in Samples A and B are: ξ_A = 0.0029(3) and ξ_B = −0.0005(3), respectively. However, these values quantify the distortion at room temperature, whereas PEME measurements are performed at 2 K. So one can expect an additional deformation in response to the specific thermal contractions of all the constituent layers embedded within the structure. In order to assess the magnitude of the temperature induced modification of the trigonal distortion of the (Ga,Mn)N unit cell upon cooling we transform ξ(300 K) into ξ(4 K) by taking into account^55,56 relevant elastic constants⁵⁷ and thermal expansion coefficients of a and c in sapphire⁵⁸ and GaN, ref. ⁵⁹. We find that for this particular material combination [the layer stack is presented in a later section (Fig. 9)]. The magnitude of the expected change of ξ upon cooling to liquid helium temperatures is marginally small. The calculated change does not exceed +0.0001, which is only a third of the typical uncertainty of establishing ξ from HRXRD. We report all these values in Table 1.

Gold pads W at the corners are the contacts to the back gate layer of n⁺-GaN:Si, the topmost layer of the substrates on which the 10 nm-thick (Ga,Mn)N layers have been deposited. Small contact pads C are deposited on 60 nm HfO₂-dielectric oxide deposited by atomic layer deposition onto the surface of (Ga,Mn)N to prevent electrically shorts via the threading dislocations. Each pad forms a parallel plate capacitor with the back gate conductive layer of GaN:Si. Upon further electrical testing only the non-leaking small capacitor (up to 12) are connected in parallel to form the final metal-oxide-semiconductor capacitor structure with capacitance of up to 1.2 nF.

Full size image

Importantly for this project, the atomic force microscopy AFM investigations yield a low surface roughness. As reported in ref. ⁶⁰ for the same Sample B and in ref. ¹⁸ for a range of sistering samples, the root mean square of the roughness of the surface does not exceed 2 nm on a 5 × 5 μm² field of view. The main reason of the surface roughening is the formation of spiral hillocks known to form at the locations where the dislocations piercing the (Ga, Mn)N layer reach the free surface, otherwise clear monoatomic steps are observed. The hillocks’ tops are these locations on the surface of GaN, which provide short-circuiting electrical passes for the vertical current across the layer⁶¹, so they have to be blocked for achieving effective gating of the material. It was previously estimated⁶⁰ that the surface density of hillocks in our samples, grown on LTTD substrates, does not exceed 10⁸ cm⁻², which represents an improvement exceeding two orders of magnitude compared to typical growth on plain GaN-templated sapphire substrates.

The growth has been performed on quarters of 2” wafers. In order to avoid a high lateral gradient of the Mn concentration, and so of the other magnetic characteristics caused by the inhomogeneous substrate temperature during the MBE growth¹⁸, for the current study the specimens are cut from the center of the quarter-wafers substrates, where the homogeneity of Mn concentration and related magnetic characteristics is the best^18,62.

Magnetic characterization

All the magnetic measurements are performed using a Quantum Design MPMS XL superconducting quantum interference device (SQUID) magnetometer following a procedure elaborated earlier^8,63,64 and adjusted for measurements of ultra-thin (Ga,Mn)N films⁶⁵ deposited on sapphire substrates⁶⁶. As shown previously for MBE-grown epilaers^17,18,20, and checked for the present samples, our high-resolution magnetic data do not reveal the presence of any high-T_C ferromagnetic precipitation.

SQUID magnetoelectric measurements

For the gate-dependent measurements Ga_1−xMn_xN layers need to be embedded in metal-insulator-semiconductor (MIS) structures. As the back gate electrode we use the conductive GaN:Si top buffer layer. To get an access to it we open by reactive ion etching (RIE) at least two 1 × 1 mm² windows in the corners of our 4 × 5 mm² specimens (marked “W” in Fig. 9). After placing photoresistive masking over these windows the whole surface is covered with 60 nm gate oxide (HfO₂) grown by atomic layer deposition. Next, after re-opening of the windows by a lift-off process, we lithographically pattern between 20 to 60 of separated gold pads over the surface covered by the gate dielectric. At the same process we deposit Au in the window areas to make the electrical contacts the GaN:Si layer—the back gate electrode(s). The gold pads are of an area about 0.2 × 0.2 mm² each and they form top capacitor plates to the back gate GaN:Si conductive layer. Figure 9 shows one of the investigated samples right after completing the patterning process. The capacitors are then electrically tested individually at room temperature. Only capacitors characterized by a high breakdown voltage V_G ≥ 20 V are connected in parallel using silver-filled electrically conductive epoxy. A resulting total capacitance of the investigated structures varies between 0.5 to 1.2 nF. A simplified general layout of the investigated structures is given in the central part of Fig. 10.

The sample is located at the center of the magnetometer’s pick up coils and the electrical gate induced changes of magnetic moment of (Ga,Mn)N are detected by phase sensitive fashion using a Lock-in amplifier.

Full size image

The completed capacitor-like structure is mounted on a modified sample holder provided by Quantum Design for performing electrical measurements in their MPMS magnetometers. In particular the modifications allows for both in plane and perpendicular orientation of the sample for electrical measurements. We perform the magnetoelectric measurements following the protocol described in ref. ⁸ with a number of small practical improvements.

The schematic diagram of the experimental set up is given in Fig. 10. Contrary to the standard method of SQUID magnetometry in which the sample is transported up and down across the SQUID pick-up coils, in our approach the sample is kept stationary at the central location of the SQUID’s 2nd order gradiometer. Having the sample stationary in SQUID a time-dependent response is provided by the application of an AC modulation voltage to the gate (capacitor) over the (Ga,Mn)N layer. The AC voltage V_AC of an acoustic frequency f (3 ≤ f ≤ 17 Hz) via the inverse piezoelectric effect modulates the magnetic properties of the (Ga,Mn)N layer, and these changes are picked up by SQUID coils. After amplification by the SQUID electronics, these signals are fed to the lock-in amplifier, where the amplitude of the changes is established in the phase sensitive fashion. We apply both the DC and AC external voltages: V_DC to change the magnetic anisotropy and, usually much smaller, V_AC to detect the resulting M_E. The method allows for detection changes down to 10⁻¹⁰ emu, however on an expense of progressively longer acquisition times. At f = 17 Hz of AC modulation it takes about 30 mins to complete each experimental point, including 5 min of initial delay to stabilize the field flux in the bore of the magnetometer’s superconducting magnet. Using about 0.1 s delay time during the data acquisition we collect about 12,000 single readouts. By plotting a histogram of the sample data, we obtain the magnetoelectric response, M_E, as the mean value of a fitted normal distribution. It should be emphasized that this technique does not sense time-independent magnetic properties of the substrate or other elements of the whole set-up. The in-phase AC magnetic response is generated only in (Ga,Mn)N located between the plates of the selected individual capacitors.

Density functional theory calculations

The first-principles spin-polarized calculations are performed using a plane-wave basis set and the Perdew-Burke-Ernzerhof exchange-correlation functional, as implemented in the Quantum-ESPRESSO package⁶⁷. The non-collinear magnetism is taken into account by performing calculations that include spin-orbit coupling and fully relativistic norm-conserving pseudopotentials. The plane-wave cutoff is set to be 50 Ry. Five different unit cells are used, namely 2 × 2 × 1, 2 × 2 × 2, 2 × 2 × 3, 3 × 3 × 2, and (3times 2sqrt{3}times 2) with 16, 32, 48, 72, and 96 atoms, respectively. In each case one Ga ion is replaced with Mn and the positions of all ions are fully relaxed. For each unit cell, the lattice parameter a is fixed to the theoretical value for GaN (i.e., to the substrate), and the c lattice parameter is fully optimized. The structural optimization is done using ultrasoft pseudopotentials and a 30 Ry kinetic energy cut-off (240 Ry for the electronic density). For the Brillouin zone sampling 8 × 8 × 8 and 3 × 3 × 3 k-point meshes are used for the smallest and largest unit cells, respectively. The calculated equilibrium values of lattice parameters for GaN are a = 0.3217 nm and c = 0.5241 nm, and agrees well with earlier theoretical reports⁶⁸.

Atomistic spin model simulations

The magnetic energy of localized spin vectors S_i is described by the Hamiltonian,

with the three terms relating to the Zeeman energy ({{{mathcal{H}}}}_{{{rm{Z}}}}=-{mu }_{S}{sum }_{i}{{{bf{S}}}}_{i}{{bf{H}}}), the exchange interaction energy in the Heisenberg form ({{{mathcal{H}}}}_{{{rm{Exch}}}}=-{sum }_{ine j}{J}_{i,j}{{{bf{S}}}}_{i}{{{bf{S}}}}_{j}), and the magnetocrystalline anisotropy energy (MAE)⁴⁸. Here, S_i is a unit vector denoting the local magnetic moment direction at the site i, μ_S = gμ_BS is an actual value of magnetic atomic moment with a spin S = 2, the Landé g-factor g = 2.0, and H is the external magnetic field. Following ref. ³⁸, we parameterize the exchange energies J_ijvs. Mn-Mn distances R_ij by the exponential function ({J}_{ij}={J}_{0}{e}^{-{R}_{ij}/b}), with b = 1.11 Å and J₀ ≅ 52.9 meV.

The anisotropy term ({{{mathcal{H}}}}_{{{rm{Aniso}}}}) contains both the trigonal (uniaxial along [001]∣∣c) ({{{mathcal{H}}}}_{{{rm{Krig}}}}) and triaxial ({{{mathcal{H}}}}_{{{rm{JT}}}}) components (refs. ^{8,40,48,69,70}). The uniaxial term includes the single ion magnetic anisotropy of the Mn³⁺ ions, driven by the trigonal distortion specific to the deviation from the perfect wurtzite atoms arrangement in the two-atom basis of GaN lattice. The triaxial term describes the influence of the Jahn-Teller deformation, which in wurtzite crystals grown along c-axis is directed along three inclined out-of-plane directions. These three Jahn-Teller centers are distorted along ({e}_{{{rm{JT}}}}^{A}=[sqrt{2/3},0,sqrt{1/3}]), ({e}_{{{rm{JT}}}}^{B}=[-sqrt{1/6},-sqrt{1/2},sqrt{1/3}]) and ({e}_{{{rm{JT}}}}^{C}=[-sqrt{1/6},sqrt{1/2},sqrt{1/3}]). It is important to remark here that it is precisely the presence of the three energy barriers that impedes in plane reversal of M and allows formation of hysteresis loops. The anisotropy terms take the form,

where K_Trig and K_JT are anisotropy parameters for uniaxial and Jahn-Teller deformations, respectively. The trigonal anisotropy energy can be expressed in more than one way due to the normalization condition ({S}_{x}^{2}+{S}_{y}^{2}+{S}_{z}^{2}=1). The lowest-order term is (-frac{1}{2}{K}_{{{rm{Trig}}}}{S}_{z}^{2}=-frac{1}{2}{K}_{{{rm{Trig}}}}(1-{S}_{x}^{2}-{S}_{y}^{2})). For K_Trig > 0 the z∣∣c axis is magnetic easy axis and the x, y plane is the hard one. Since constant term can be omitted, we have opted for the form presented in Eq. (6).

For the modelling of the magnetoelectric effect, the trigonal term is of a particular importance. It is exactly this term of ({{mathcal{H}}}) which is affected by the application of the electric field E and we compute the magnitude of gate driven effects through the relation K_Trig(E) ∝ ξ(E).

Simulations of the dynamical properties of the system are based on the stochastic Landau-Lifshitz-Gilbert (LLG) equation applied at the atomistic level^26,27,48, which is given by

where γ is the gyromagnetic ratio and α_G = 1 is the precession damping term. The total effective magnetic field acting on spin i is composed of two terms ({{{bf{H}}}}_{{{rm{eff}}}}^{i,tot}={{{bf{H}}}}_{{{rm{eff}}}}^{i}+{{{bf{H}}}}_{th}^{i}) (refs. ^25,26,27). The net effective field ({{{bf{H}}}}_{{{rm{eff}}}}^{i}=-(1/{mu }_{S})partial {{mathcal{H}}}/partial {{{bf{S}}}}_{i}) is obtained from the derivative of the total Hamiltonian with respect to the atomic moment μ_SS_i. The thermal field in each spatial dimension is represented by a normal distribution Γ(t) with a standard deviation of 1 and mean of zero, and is given by

where Δt = 5 ⋅ 10⁻⁶ ns is the integration time step and T is temperature.

In atomistic spin model simulations, the parallelization of the code is implemented. We divide the simulated region into sections, where each processor simulates specified part of the complete system. The computations are accelerated using graphics processing units (GPUs). We use Euler’s method as the integration scheme. In order to speed up the calculations, strict synchronization of the processors (cores) is applied every fifth step of the LLG integration. The scaling and parallelization of the code is found to be very good. We obtain practically the same results independent of number of nodes or cores used in simulations as well as independent of used computer architecture (based on GPU or CPU). As the velocity of magnetization precession depends on the strength of the magnetic field, we assume that the initialization and averaging steps of the LLG simulations depend on the field value, as given in Table 2.