Simulations of magnetization reversal in FM/AFM bilayers with THz frequency pulses

It is widely known that antiferromagnets (AFMs) display a high frequency response in the terahertz (THz) range, which opens up the possibility for ultrafast control of their magnetization for next generation data storage and processing applications. However, because the magnetization of the different sublattices cancel, their state is notoriously difficult to read. One way to overcome this is to couple AFMs to ferromagnets—whose state is trivially read via magneto-resistance sensors. Here we present conditions, using theoretical modelling, that it is possible to switch the magnetization of an AFM/FM bilayer using THz frequency pulses with moderate field amplitude and short durations, achievable in experiments. Consistent switching is observed in the phase diagrams for an order of magnitude increase in the interface coupling and a tripling in the thickness of the FM layer. We demonstrate a range of reversal paths that arise due to the combination of precession in the materials and the THz-induced fields. Our analysis demonstrates that the AFM drives the switching and results in a much higher frequency dynamics in the FM due to the exchange coupling at the interface. The switching is shown to be robust over a broad range of temperatures relevant for device applications.


Introduction
Antiferromagnets (AFM) are materials with alternating magnetic moments at the atomic level resulting in no net magnetization.Their inherently fast THz magnetization dynamics [1,2] and absence of stray fields potentially means faster magnetization reversal and higher information density which makes them strong candidates for next generation spintronic and magnetic recording media.Due to the antiferromagnetic order, however, controlling them remains a challenge as there is no net Zeeman energy when placed in a magnetic field, which complicates device development, although this is starting to be overcome by taking advantage of spin-orbit effects, e.g.spin-orbit torques.
The use of current-induced Néel spin-orbit torques (NSOT) has been shown to be able to induce magnetization dynamics in several AFMs.Two materials, CuMnAs [3][4][5][6] and Mn 2 Au [7][8][9][10], have been identified with a strong spin-orbit coupling that creates sufficient NSOTs for manipulation of the Néel vector.Experiments so far have not taken advantage of the intrinsically fast dynamics of AFMs, with switching via electrical control being limited by the low frequency mode of the antiferromagnet [11].Computationally, It has been shown that square staggered fields can lead to reversal in Mn 2 Au [12][13][14] with pulse durations of the order of a few picoseconds.While a square field is a good first approximation, experimental realisation of such a profile is challenging.However, advances in the experimental generation of ultrashort THz pulses ´ µ ´ µ 0 1 2 ´ µ ´ µ ´ [ [15][16][17][18] has opened the possibility for excitation and reversal of the Néel vector using sub-picosecond single and multi-cycle fields.
Even though the manipulation of the order parameter can be carried out electrically via NSOTs, the opposing net zero magnetization in AFMs means the readout of the order parameter, otherwise known as the Néel vector, still remains a challenge.The opposing sublattices mean that detection of the magnetic order using conventional methods such as the anisotropic magnetoresistance mechanism (AMR) often yields weak signal outputs compared to typical ferromagnetic materials [11,19].However, a recent demonstration of 100% magnetoresistance in an antiferromagnetic stack of MnPt and Mn 3 Pt [20] which could pave the way for future antiferromagnetic spintronic devices.Another approach of inferring the Néel vector is to couple AFMs directly to FM materials and measure the orientation of the AFM through measurements of the magnetization state of the FM.Recently, it has been shown Mn 2 Au can be coupled ferromagnetically at the interface to Permalloy with one-to-one imprinting of the AFM order on the FM domains [10], though little is known of the dynamics of switching and whether one can simultaneously take advantage of the higher frequency dynamics of the AFM whilst being coupled to the sluggish FM.
In this work, we conduct atomistic spin dynamics (ASD) simulations of pure Mn 2 Au as well as an FM/AFM bilayer consisting of a Mn 2 Au coupled to ferromagnetic Permalloy (Py).Sec. 2 introduces the atomistic model, the applied magnetic field profile and material properties used for Mn 2 Au.ASD simulations of reversal in Mn 2 Au as a result of a THz frequency field at 0 K and 300 K are then presented in Sec. 3 where we show that, due to the tetragonal anisotropy in the plane, the magnetization can be reoriented along the [110], [ 110], [ 1 10] and [ 11 0] directions.We show that staggered linearly polarized fields of the order of ∼ 1 Tesla are required for deterministic reversal.In Sec. 4 we determine phase diagrams for magnetization reversal in Mn 2 Au/Py bilayers with varying interface exchange and Py film thicknesses.We compare cases of excitation of both the Mn 2 Au and Py sublattices to excitation of Mn 2 Au alone to probe the relative roles of each layer in the switching.Finally, we simulate switching at 300 K and 600 K to determine whether the temperature scaling of the anisotropy affects the field strengths and pulse durations required for switching.

Methods
The spin dynamics of the magnetic moments is described by the well known stochastic Landau-Lifshitz-Gilbert (LLG) equation: where S i is a normalised unit vector of the spin at site i, λ i is the effective damping parameter, γ = 1.76 × 10 −11 is the gyromagnetic ratio, µ i is the atomic magnetic moment and H i is the effective field acting on the spin at site i.The effective field is given by the equation: where ζ i (t) describes the coupling to the thermal bath and B(t) is the applied external field.For atomistic simulations of pure Mn 2 Au, the following Heisenberg Hamiltonian is used: with anisotropy constants taken from Ref. [12] d z = −1.19meV, d zz = −0.015meV and d xy = 0.04 meV following the original calculations in Ref. [21].We use a magnetic moment of µ Mn = 3.87µ B .The exchange interactions are taken from Ref. [22].We use a value of λ = 0.01 in Mn 2 Au as seen in previous atomistic modelling [12,13,22].The Mn sublattices are initialised along the [1,-1,0] and [-1,1,0] directions respectively.The LLG equation is integrated using the Heun method [23].
The staggered field, B(t), is applied along the perpendicular direction in the easy-plane to maximise the torque.Such staggered fields can be generated in metallic antiferromagnets such as Mn 2 Au as a result of electrical currents [8,24,25] .The field is modelled as a Gaussian envelope with fixed frequency of 0.67 THz, which is close to the in-plane resonant resonant mode.The field equation is given by: Where H is the magnitude of the Gaussian envelope, σ is the standard deviation, and f is the frequency of the pulse.The exchange for Mn 2 Au consists of the 13 interactions for each site of which 9 are AFM coupling and 4 correspond to coupling in the FM planes.This set of exchange parameters is the same as used in Ref. [22] and yields a Néel temperature of ≈ 1350 K. Experimental work puts the Néel temperature between 1300 and 1600 K [26].

Reversal in Mn 2 Au
We firstly consider a system of pure Mn 2 Au subject to a staggered field at 0 and 300 K. Initial calculations (not presented here) show that the switching at T=0 K is precessional and therefore we can simulate a single unit cell.However, in the room temperature simulations we simulate 70 × 70 × 70 unit cells of Mn 2 Au, each consisting of 4 atoms yielding a total system size of 1,372,000 Mn atoms with periodic boundary conditions at every surface.We simulate 9 values of H in Eq. ( 4) varying from 80 to 720 mT in steps of 80 mT.It is worth mentioning that H does not correspond to the peak field amplitude, B max .A value of H = 720 mT for σ = 0.3 ps yields a maximum field amplitude of B max ≈ 442 mT.The standard deviation and frequency remain fixed in this instance at σ = 0.3 ps and f = 0.67 THz respectively.Fig. 1 shows the Néel vector reorientation in the easy-plane following the THz pulse pulse at 0 K (left pane) and 300 K Fig. 3: A snapshot of the system close to the interface.The Py atoms are shown in green, the Au in black and the Mn sublattices in blue and red respectively.The Py and Mn is coupled ferromagnetically at the interface.Each axis has been scaled differently to aid in the visualisation.(right pane).It is worth noting that the out-of-plane angle remains close to zero through the entire process due to the large uniaxial easy-plane anisotropy constant d z in Eq. ( 3).The onset of 90°switching occurs B max = 294 mT with 180°being observed at B max = 442 mT for T = 0 K.For 300 K, the simulations were repeated 8 times for all field values to account for the stochastic thermal effects.Dotted lines show the results from individual simulations, and shaded area represents the upper and lower bound for each field value.Upon comparing both panes in Fig. 1, the first instance of reversal has reduced from B max = 294 mT for 0 K, to B max = 245 mT for 300 K due to the reduction in anisotropy field.In the room temperature simulations, the switching was deterministic; however, the path and final reorientation angle were not -as seen for a field amplitude of B max = 392 mT, where both 90°and 180°switching occurred.In the study by Selzer et al. [12], much smaller switching fields of ≈ 76 mT were required for deterministic switching close to room temperature although it is worth noting that this is as a result of a 'step-like' staggered field compared the cyclic field profiles used in Fig. 1.Previous attempts to characterise the staggered spin orbit field in terms of an electric field give varying strengths.In Ref. [12], Electric fields of 10 7 V/m yield staggered fields of about 76 mT.Much lower magnitudes of the staggered field were calculated in Refs.[8,13] but orbital contributions were not taken into account.For the THz staggered switching fields of around ≈ 300 mT discussed in the previous paragraph, the electric fields required to generate such a field should be achievable in experimental facilities [27][28][29].
For precessional switching, it is possible to describe the easy-plane magnetization dynamics of Mn 2 Au using a single equation of motion [30] 1 where ϕ L is the angle of rotation in the easy plane, B is the staggered field given in equation ( 4), α G is the Gilbert damping constant, H an and H ex are the anisotropy and exchange fields respectively, σ c is the conductivity and λ NSOT is the torquance from the staggered field.By fitting the 0 K magnetization dynamics shown in Fig 1 to Eq. ( 5), a value for λ NSOT can be extracted and compared directly to experimental results irrespective of the field profile.Fitting to the 0K dynamics, using a conductivity of 1.5 A selection of the fits used to extract λ NSOT are shown by the empty circles in the left hand pane of Fig. 1.A more detailed discussion of the fitting process can be found in Supplementary Information S1.
For 0 K, it is possible to sample a large phase space of B max and σ with little computational cost.In Fig. 2 we present a switching phase diagram at T = 0K for a range of applied field amplitudes and pulse durations for a fixed frequency of 0.67 THz.The color shows whether the system undergoes 90 or 180 degree switching (relative to the initial state) following a THz staggered pulse.The field and pulse duration required for Néel vector reversal increases for a multi-cycle signal because of the changing sign of the external field.For σ 0.4 ps, the switching window appears less structured, which is likely due to the transition from a single-cycle to a multi-cycle field with increasing σ.What is surprising is that pulses as short as 0.2 ps with sub-Tesla staggered fields can drive antiferromagnetic switching.The phase diagram is more complex than that seen previously for antiferromagnetic NiO [31] due to the use of an un-staggered square field profile with a single uniaxial anisotropy constant in that study.
4 Reversal in Mn 2 Au/Py Bilayers Now that we have established the feasibility of ultrafast THz switching in pure Mn 2 Au, in this section we demonstrate magnetization reversal of ferromagnetically coupled Mn 2 Au/Py bilayer structures.The introduction of Py requires further material parameters, increasing the complexity of the parameter space.
Here we focus on: (i) thickness of the Py (ii) thickness of Mn 2 Au and (iii) the coupling at the interface.In this work, we keep the thickness of Mn 2 Au fixed at 25 units cells along the length of the system equating to roughly 21 nanometers (nm).We chose to vary the thickness of Py using 6, 8, 10 and 20 unit cells.This corresponds to thin film thicknesses of approximately 2.13 nm, 2.84 nm, 3.55 nm and 7.1 nm respectively -similar thicknesses of Py used in the Mn 2 Au/Py bilayer experiments of Bommanaboyena et al. [10].We include open boundary conditions at either end of the chain and periodic boundaries on all others.The Py is stacked on top of the Mn2Au in the (001) direction such that the the same AFM sublattice couples to the FM at the interface.For Mn 2 Au we use lattice constants of a Mn 2 Au = 3.33 Å and c Mn 2 Au = 8.537 Å.
We also assume the Py matches a Mn 2 Au given a comparable value of a Py = 3.55 Å for bulk Permalloy.A schematic of the chain close to the interface can be found in Fig. 3.
For the Py, we use a simplistic model with nearest-neighbour (n.n.) totalling 12 interactions per Py site.At the interface, exchange is treated as nearest neighbour at the two neighbouring atoms at the interface.In the simulations, the exchange coupling at the interface is varied from 0.5 × 10 −21 J (≈ 15% of the n.n J ij used for Py) to 5.0 × 10 −21 J.The interface coupling is chosen as a variable quantity in our simulations as it can be manipulated experimentally by the insertion of a non-magnetic spacer between the FM and AFM [32] or doping at the FM/AFM interface.For the simulations of Mn 2 Au and Permalloy bilayers, there are two additional Hamiltonian terms: The first term in the above can be found in Eq. ( 3).For the Permalloy, we use an average moment of 0.95µ B with nearest neighbour exchange and exclude any anisotropy.For the 12 neighbouring interactions in Py, we use a value of J Py ij = 3.01 × 10 −21 J which yields a Curie temperature ≈ 720 K.
The applied field is once again staggered for each Mn sublattice.An identical field to the sublattice containing the Mn atom at the interface is applied to the Py sublattice.Firstly, we consider the case of Py at 6 unit cells, and vary the strength of the coupling across the interface.Fig. 4 shows four different coupling strengths, namely 0.5, 1.0, 2.0 and 5.0 × 10 −21 J.For low σ, the phase diagram varies minimally in structure across an order of magnitude increase in the interface exchange.While this region remains similar, there is a significant change in the phase diagram for σ 0.2 which transitions from small and scattered switching windows, to larger more continuous regions of magnetization reversal.Simulations are also conducted for constant exchange coupling with varied Py thickness, we chose the largest value of J Inter ij = 5.0 × 10 −21 J and vary the thickness from 6 to 20 unit cells, similar in range to the bilayer experiments in Ref. [10].As was seen in Fig. 4, the phase diagram remains similar in shape for σ 0.2, but there is a notable increase in the fields required for reversal in this region with almost no 180°reversal observed for the case of 20 Py unit cells.For σ 0.2, the change in the phase diagram as we transition from 6 to 20 unit cells is more disordered.In contrast to switching in pure Mn 2 Au, there exist large bands in Figs 4 and 5 where no magnetization reversal is observed, which suggests the FM is hindering switching.
To further reiterate this point, simulations were conducted where only the Mn sublattices were subject to an applied staggered field.Fig. 6 shows the change in phase diagram as a result of an excitation of the entire system compared to the case of excitation only in the AFM.For σ 0. Changes are only observed at the boundary between two reorientation angles by, in most cases, a single shift in the field strength or pulse width increment.The AFM is therefore almost almost entirely responsible for the switching in this region.For larger pulse widths, the FM is given more time to respond to the applied field and therefore the role of the field in the FM becomes more important as shown by the large changes in reorientation for σ 0.3.
It's well understood that elevated temperatures lead to a reduction in anisotropy fields [33][34][35].Such heating is the key principle behind Heat Assisted Magnetic Recording (HAMR) [36][37][38][39].Here, we conduct finite temperature simulations for the bilayer containing 6 unit cells of Py along the chain with an interface exchange of 0.5 × 10 −21 J -which shows the most favorourable switching characteristics.Due to the increased computational cost of finite temperature ASD simulations a smaller phase space is sampled.We simulate in the ranges 0.1 ≤ σ ≤ 0.2 ps and 0 ≤ H ≤ 6 T (recall H does not correspond to the maximum amplitude, B max ).Simulations are repeated 8 times for each field strength and pulse width.Fig. 7 shows the switching probability for the two aforementioned temperatures where there is a lowering and broadening to the switching band transitioning from 300 to 600 K.At 300 K, the equilibrium magnetization is roughly 0.84 M/M S and 0.92 M/M S for the Py and Mn 2 Au sublattices respectively.At 600 K, it is approximately 0.42 M/M S and 0.81 M/M S .The area inside the red box shows the finite temperature probability.The area outside the red box shows switching at 0 Kelvin.The switching is always deterministic at 0K as there is no stochastic noise processes.There is no differentiating between 90°and 180°reversal in Fig. 7; what is shown is the the probability of either event occurring.
The reversal path for the Mn sublattices remains circular in the xy-plane across all temperatures because of the large negative uniaxial anisotropy constant d z .The Py reversal path is also circular with no significant reduction in the magnetization following the application of the field.However, interestingly, the lack of anisotropy allows for the precession of the Py sublattice out of the xy-plane.Fig. 8 shows the magnetization dynamics for the two most common reversal paths for an interface exchange of 0.5 × 10 −21 J with a Py thickness of 6 unit cells.Because of the highly circular precession in the x-y plane resulting from the strong uniaxial anisotropy d z , the reversal pathways are derived from the dynamics of the Mn sublattices, making it simple to characterise the path in terms of a rotation about the z axis.Fig. 8a shows the dynamics corresponding to a 270°clockwise rotation of the Mn sublattices and 8b shows a reorientation after a clockwise rotation of 180°.In total, these two paths account for 51% of all the switching events at 300 K with 8a and 8b accounting for 29% and 22% of the switching events, respectively.If we characterise the reversal time as the time taken for the Mn sublattices to reach the final steady state, the average reversal time is ∼ 5 ps and ∼ 9 ps for Fig. 8a and 8b respectively.There exist 8 other paths with occurrences of between 1% and 9% that account for the the remaining 49% of switching events.While the switching events occur in These two paths account for 52% of all switching processes at 300 K.
only a few picoseconds, the system continues to oscillate at its resonant frequency for an extended period of time.The frequency of this resonance is characterised by both the thickness of the Py layer and the strength of the interface coupling.Increasing the thickness of Py reduces the resonant frequency, whereas increasing the coupling strength, J Inter ij , at the interface causes an increase in the resonant frequency.Both of these quantities could act as a way to tune to desired frequencies for experimental applications.Further discussion of the resonant frequency can be found in Supplementary Information section S2.

Discussion
Here we have presented results for magnetization reversal in pure Mn 2 Au as well as an FM/AFM bilayer consisting of an average moment model of Permalloy and Mn 2 Au.We have shown that reversal across the bilayer is possible as a result of sub-picosecond THz fields applied in an in-plane direction that maximises the torque, with amplitudes no larger than a few Tesla.The magnetization dynamics have been used to accurately parameterise the torquance, λ NSOT , which can be directly compared to experiment irrespective of the field profile.The importance of these findings is based on the ability to generate staggered fields in Mn 2 Au via NSOTs and, more practically, create FM/AFM bilayers which have one-to-one mapping of the Néel vector on the FM domain structure as seen in Ref. [10].In this work we have focused on several important physical parameters including the field amplitude, field standard deviation, the interface exchange and the thickness of the FM layer.It is also worth noting that the FM used in the bilayer system should be treated as generic and could easily be interchanged with another weakly anisotropic material.Nonetheless, we have shown that reversal in FM/AFM Bilayer systems are possible using sub-picosecond fields with amplitudes of ∼ 2 Tesla.Regions of switching have been identified where staggered fields in the AFM alone are responsible for magnetization reversal.Simulations at finite temperature show a small reduction in the fields required for switching.A choice of FM with a similar critical temperature to Mn 2 Au would allow switching at further elevated temperatures, likely resulting in further reductions of the required field amplitudes.

Fig. 1 :
Fig. 1: Magnetization reversal in Mn 2 Au for varying maximum field amplitudes, B max , from approximately 50 to 450 mT at (a) T = 0 and (b) T = 300 Kelvin for a constant pulse width and frequency of σ = 0.3 ps and f = 0.67 THz.Inset in (a) shows the field profile between 0 and 1 ps or a σ = 0.3 ps.The dotted lines in (b) represent sublattice dynamics from individual simulations, the shaded area corresponds to the upper and lower bound for each field strength.

Fig. 2 :
Fig. 2: Switching phase diagram for pure Mn 2 Au following a THz pulse.Turquoise and yellow regions show regions of 90°and 180°switching respectively.
2 ps, there is little observed change between Fig 4a and Fig 6a, as well as 4d and Fig 6b.

Fig. 6 :Fig. 7 :
Fig. 6: Changes in phase diagram from an excitation of both Mn and Py sublattices to just Mn sublattices for coupling strengths of (a) 0.5 × 10 −21 J and (b) 5.0 × 10 −21 J. Dark blue areas represent areas of no change in switching between the two cases.Other possible scenarios can be found in the colorbar.

Fig. 8 :
Fig. 8: Magnetization dynamics at 300 K for the two most common reversal paths.(a) an anti-clockwise rotation of 270°in the xy-plane for the Mn sublattices and a clockwise rotation of 90°for the Py.(b) an anti-clockwise reorientation of 180°in the xy plane for Mn sublattices and a rotation of 180°out of the xy plane such that the Py sublattice is parallel to the z axis roughly halfway through the switching process.These two paths account for 52% of all switching processes at 300 K.