Switching the spin cycloid in BiFeO3 with an electric field

Bismuth ferrite (BiFeO3) is a multiferroic material that exhibits both ferroelectricity and canted antiferromagnetism at room temperature, making it a unique candidate in the development of electric-field controllable magnetic devices. The magnetic moments in BiFeO3 are arranged into a spin cycloid, resulting in unique magnetic properties which are tied to the ferroelectric order. Previous understanding of this coupling has relied on average, mesoscale measurements. Using nitrogen vacancy-based diamond magnetometry, we observe the magnetic spin cycloid structure of BiFeO3 in real space. This structure is magnetoelectrically coupled through symmetry to the ferroelectric polarization and this relationship is maintained through electric field switching. Through a combination of in-plane and out-of-plane electrical switching, coupled with ab initio studies, we have discovered that the epitaxy from the substrate imposes a magnetoelastic anisotropy on the spin cycloid, which establishes preferred cycloid propagation directions. The energy landscape of the cycloid is shaped by both the ferroelectric degree of freedom and strain-induced anisotropy, restricting the spin spiral propagation vector to changes to specific switching events.


Introduction
In the push to realize spintronic devices for computation, controlling magnetism using electric fields at the device scale remains an elusive challenge 1,2 .The multiferroic BiFeO3 (BFO) exhibits both room-temperature large ferroelectric polarization and canted antiferromagnetism with deterministic coupling between these order parameters through the crystal symmetry [1][2][3][4][5][6] .Thus, the magnetic nature of BFO can be controlled with electric fields, potentially paving the way for ultra-energy-efficient magnetic and spintronic devices with more favorable scaling [6][7][8][9][10] .Indeed, electric-field control of magnetism in coupled ferromagnetic layers using BFO been the subject of considerable attention 4,[11][12][13] , and while this functionality has long excited researchers, understanding of the process has generally relied on mesoscale imaging and transport measurements to infer the structure and behavior of the BFO.This has left gaps in our understanding of the microscopic mechanism and details of this process.
Only recently has the opportunity to directly image and understand complex magnetic structures at their native scales become possible using techniques such as nitrogen-vacancy (NV) diamondbased scanning probe magnetometry (henceforth NV microscopy) 14,15 .While recent works have shed new light onto the nanoscopic magnetic structure of BFO 16,17 , a significant fundamental question remains as to how the complex magnetism interacts with the ferroelectric order and electric fields.In BFO, ferroelectric switching pathways are complex 11,18,19 , proceeding through coupled, multi-step rotations of the polarization.The corresponding magnetic-switching phenomena can be tied to the ferroelectric polarization evolution 4,12,20 , but observation of these coupled phenomena in real space on the nanometer scale has not been possible.Here, using NV microscopy, we directly image the stray magnetic field at the surface resulting from the spin cycloid as it couples to ferroelectric domains and complex (71°, 109°, and 180°) ferroelastic and ferroelectric switching events.
In the bulk, BFO exhibits an antiferromagnetic (AFM) spin cycloid that exists in the plane defined by  and the propagation direction, , which points along a ⟨110⟩ that is orthogonal to  [21][22][23][24][25] .This vector connects second-nearest neighbor iron sites which, in an unperturbed G-type AFM, would be ferromagnetically coupled within a {111} [25][26][27] .The spin cycloid itself has been modeled as a Néel-type, rotating uncompensated magnetization, (), that exists in the plane defined by  and  (i.e., the (112 ̅ ), where  is along the [111] and  along the [1 ̅ 10], unless otherwise noted) with a period of ~65 nm (Figure 1a).This has been described as where  is the volume-averaged magnetization, || = 2/,  is a coordinate in real space, and   and   are the unit vectors in the directions of  and , respectively 14 .A second-order canting also exists due to the Dzyaloshinskii-Moriya interaction (DMI) arising from the antiferrodistortive octahedral rotations 26,[28][29][30] , causing the magnetization to buckle slightly out of the - plane (Figure 1b).This second-order spin-density wave, noted here as   , can be described 16 by ( Example solutions to eqs. 1 and 2 are provided in Supp. Figure 1.The direction of , and the magnetization are thus intimately tied to the ferroelectric polarization of BFO, consistent with the fact that the spontaneous polarization is the primary order parameter of multiferroicity in this system.In NV microscopy, the stray magnetic fields (()) from the sample surface perturb the energy states of a single NV center implanted into a diamond scanning probe tip.By measuring the optically detected magnetic resonance spectra of the NV center, one can detect small magnetic fields (<2  √ −1 ) with spatial resolution down to at least 10 nm.Combined with piezoresponse force microscopy (PFM), it is possible to locally map the vectorial polarization distribution, (), in the ferroelectric domains and thereby correlate the relationship between the magnetic and ferroelectric structure in BFO.While the interaction of the cycloid with ferroelectric domains has been suggested in previous work, where poling areas of the film results in a locally uniform  (the cycloid propagation direction), the interaction between ferroelectric domain walls and the cycloid has not been directly demonstrated 14 .More importantly, it is not well understood if and how the cycloid propagation direction changes during ferroelectric switching, a question especially relevant to electric-field manipulation of magnon transport 20,31 and exchange coupling across heterointerfaces 4,11 .Here, it is shown that direct mapping of the canted antiferromagnetic texture to the ferroelectric domains can be achieved.Of greater importance 20,31 , it is shown that upon applying an electric field to switch the ferroelectric polarization, the relationship between  and  is conserved, with  showing a strong anisotropy along the in-plane [110] and [1 ̅ 10] which are both perpendicular to .It is also observed that this anisotropy persists agnostic to the direction of electric field, through both in-plane and out-of-plane ferroelectric switching events.

Anisotropy of the cycloid propagation
The model heterostructures studied here are ~100 nm thick, (001)-oriented BFO thin films deposited on DyScO3 (DSO) (110) substrates using pulsed-laser deposition, both with and without metallic SrRuO3 bottom electrodes (Methods).In both configurations (with and without bottom electrodes), X-ray diffraction confirms that heterostructures are constrained by the substrate (Supp.Figure 2).Subsequent PFM imaging reveals 71° stripe-like ferroelectric domains (Figure 1c); a commonly observed characteristic 32 of BFO.In turn, NV microscopy of the same shows the ~65 nm period sinusoidal modulation of magnetization due to the spin cycloid 14,17 .By carrying out PFM and NV microscopy at the same location, we observe that the 71° ferroelectric domain walls match exactly with changes in the cycloid propagation, switching from  1d).In bulk BFO, all three of these variants are allowed 16 .In the thin film heterostructures studied here, however,  appears to select only the purely in-plane directions, [110] or [1 ̅ 10], switching every ~200 nm commensurate with the ferroelectric domain walls (Figure 1e,f).
In the as-grown configuration,  is observed to be along both the [110] and [1 ̅ 10], in each case selecting the allowed direction without an out-of-plane component dependent on the ferroelectric polarization.We hypothesize that the preference for the in-plane  vectors arises from epitaxial constraint imposed on the BFO film from the DSO substrates (e.g.differences in symmetry, lattice constant, etc.).It has been previously observed that, in thin-film BFO, the relationship between  and  can be controlled through the choice of substrate and, at small compressive strains, the cycloid follows similar behavior to bulk 14,16,24 , where  is bound to a ⟨110⟩ that is orthogonal to  (so-called type-I).At higher epitaxial strains, however, a type-II cycloid is observed, which propagates along a ⟨112 ̅ ⟩ perpendicular to  16,17,33 .It appears, then, that while constraints imposed by the DSO substrate allow for the bulk-like, type-I cycloid, it creates a strongly anisotropic landscape which inhibits out-of-plane projections of .This effect on the anisotropy of , however, has not been determined.To better understand the impact of this anisotropy on the spin texture, density functional theory (DFT) calculations have been performed to explore the canted magnetic moment in BFO.
While DFT simulations can provide powerful insight into the local atomic structure, predicting the ground-state magnetic order of BFO is computationally intractable using plane-wave DFT due to the long period of the cycloid (~65 nm) and the corresponding large system size required to fully simulate it 34 .A practical alternative is to use the generalized Bloch condition for q-spirals which can account for rotations in the cycloidal plane, but this cannot predict moments that are canted out of this plane due to boundary conditions.To then help understand the ground-state magnetization of the system, here we discretize the magnetic structure into subsections manageable by first principles.Starting from a 2x2x2 unit cell with G-type antiferromagnetism, where anti-aligned iron moments point along the vector ±, we systematically rotate the initial  within a fixed plane (given by the angle , Figure 2a) to simulate the cycloid and resolve the canted magnetism.In this case,  is in the (112 ̅ ), defined by  and , which hosts the cycloid.
Additionally, calculations are performed applying on-site Hubbard U and Hund J corrections to the O-p manifold, in addition to the Fe-d and Bi-p, with further detail in Methods.We first demonstrate the validity of our approximations and methodology by reproducing the established spin texture, , as reported in the literature.Specifically, we rotate  within an orthogonal plane, (111) which leads to an approximately two-times greater energy cost than  in (112 ̅ ), relative to the reference minimum energy spin quantization axis (Figure 2b). and is reported as the vector sum of the iron spins.From our simulations,   follows the same period as the cycloid and reaches a maximum value when  is parallel to [1 ̅ 10].This is consistent with the expectation from symmetry that   emerges due to the DMI from octahedral rotations with their axis along the polarization direction 26,29,30 , where   • (  ×   ) is maximized when   and  , are orthogonal, in this case  , || .The value of   reaches a maximum of 0.02   , which is consistent with previous predictions 34,35 and is of the order of previous experimental results 14,16,22,29 .Discretizing the magnetic cycloid in this way then produces results that agree exactly with previous experimental and theoretical interpretations of the canted   component of the cycloid [28][29][30] .Additionally, reducing the ferroelectric polarization or rotating the cycloid in another plane removes this phenomenon (Supp.Figure 3), showcasing its physical origin.
With the magnetic structure reproduced using DFT, we next consider how the anisotropy introduced into the system through epitaxial constraints.Fixing the unit cell to the in-plane lattice constants of DSO, a similar calculation is performed to that above where the iron moments are rotated in the planes defined by  and the three possible  directions (i.e., [1 ̅ 10], [1 ̅ 01], and ).From these data, there is a clear anisotropy favoring the  ∥ [1 ̅ 10] due to epitaxy, with a mean energy of approximately two-times lower than the other symmetry allowed , in line with the experimental results.This is in contrast to the bulk, zero-strain state, where all possible directions of  are symmetry and energetically equivalent (Supp.Figure S5).Though, as mentioned above, we cannot accurately simulate the interaction between neighboring iron moments along  that give an additional potential energy gain due to the cycloid 26 , it is noteworthy that in this colinear limit of the antiferromagnetic structure, we would predict a ground state with  pointing along [1 ̅ 10] with  along [111]; the same as suggested previously 11 .Understanding, then, that substrate constraints intrinsically break the degeneracy of the allowed  directions, one can then ask whether the state of the cycloid can be deterministically changed with an electric field and whether this anisotropy is strong enough to persist in the switched state.

In-plane electric field switching of the cycloid
Through application of an in-plane electric field, perpendicular to the striped 71° domain walls, we  As in the case of the as-deposited films,  remains parallel to the surface of the film.If the polarization of a domain reorients from [111] to [1 ̅ 11],  selects the [110] out of the three symmetry allowed axes.We believe that this may be due to the biaxial strain state imposed by the DSO substrate, as this strong in-plane anisotropy is not observed in BiFeO3 deposited on other REScO3 substrates 16,17 , such as GdScO3, where the  || [011] propagation appears to be favored.With this assumption, that  is confined to the in-plane directions and we will observe a 90° change in  under 71° ferroelectric switching, we measure the change in  in situ at a single location under electric field.
With an in-plane electric field, individual ferroelastic domain walls will tend to remain stationary and the polarization of individual domains reorients by 71° (Figure 3e) 4,11,20 .The stray magnetic field measured at a single location in situ and the corresponding ferroelectric domains are shown (Figure 3f-k) where, e.g., in the center domain, the directionality of  changes from [110] to [1 ̅ 10]   and back to [110] under successive switching events.Again, after ferroelectric switching, the relationship between  and  is conserved, with both rotating 90° about the [001]-axis.In the case of BFO, this observation is not necessarily surprising.In previous works, electrical switching of the polarization in BFO has been proposed to follow successive rotations of  (i.e., ferroelastic switching) instead of through an intermediate state which is a nonpolar configuration (i.e., ferroelectric switching) 11,18 .From this framework of the deterministic rotation of  and the FeO6 octahedra, the case of an in-plane field becomes easy to understand: if  rotates 90° about the [001] (as in a 71° switch), and  ⊥ , we would expect that  will also rotate correspondingly by 90° about the [001].

Out-of-plane electric field switching of the cycloid
In the case of an electric field applied along the [001], however, this switching pathway becomes more complex.Using a BFO heterostructure deposited with a SrRuO3 back electrode, we can explore how electric fields in the out-of-plane direction, and more complex ferroelectric switching events, influence the reorientation of the spin cycloid.An out-of-plane ([001]-oriented) electric field was applied locally to the sample using the PFM tip as the top electrode while grounding the bottom SrRuO3 layer (Methods).A voltage of +8 V was applied to the heterostructure within a ~5 m area, and -8 V was applied to a ~2.5 m area within the previously switched region to return it to the same   as the as-grown state (a so-called box-in-a-box structure).PFM is then used to measure the switched area at multiple scan angles to vectorize the polarization (Methods and Figure 4a,b).Using the polarization vectors and calculating the angle between them before and after switching allows us to construct a map of ferroelectric switching events (Figure 4c).This data can then be compared to NV magnetometry measurements to study the dependence of the cycloid propagation direction on local switching events.4f,i) that there are 71°, 109°, and 180° ferroelectric switching events during the poling process.As in the as-grown material,  is confined to the [110] and [1 ̅ 10] within the (001) and is oriented perpendicular to the final state of , regardless of the local switching event being a 71°, 109°, or 180° rotation of the polarization.This is true for both the singly and doubly switched areas.The interpretation is then that  is dominated by the magnetoelastic energy from the substrate constraints and the final polarization vector direction.
Previous reports on BFO have demonstrated that, when subjected to an out-of-plane oriented electric field,  will not rotate continuously through 180° or through a nonpolar intermediate state, but switching will be accomplished by successive 71° and 109° switching events 11,18 .The data here mirror these observations, but the anisotropy of  is showcased by the fact that the highest symmetry operations to rotate  do not preserve the experimentally observed directions of .
These potential pathways and the operations to preserve the directionality of  are shown (Supp.

Figure S8
).In other words, if  were not anisotropic in the (HK0), it would be expected to change with the highest symmetry operation of .This observation mirrors previous results on the deterministic nature of the switching in BFO 11 , where switching only proceeds along experimentally observed pathways if the antiferromagnetic order parameter, , is fixed to the ⟨110⟩ in the (001), the same as  observed here.
Here we show that the antiferromagnetic spin cycloid in BFO is intimately tied to the direction of polarization, even after ferroelectric switching events.In BFO thin films deposited on DSO, the propagation vector of the cycloid, , is confined to directions orthogonal to  with no out-of-plane component, [110] and [1 ̅ 10].Discontinuities in the cycloid propagation map directly to 71° ferroelectric domain walls, where  must reorient to maintain this relationship.As the magnetic structure can be not only tuned in the ground state, but the functionality can be constrained inplane using the substrate (e.g., symmetry, lattice mismatch), this may provide a method for designing particular material architectures to optimize spin and magnon transport 20,31 .The fixed in-plane anisotropy of the antiferromagnetic cycloid is especially important for the implementation of vertically oriented spintronic devices 9,10 .Additionally, the ability to potentially engineer local frustration in the magnetic structure may open the door for electrical creation of magnetic topological defects at particular domain wall interactions [36][37][38][39] , as confinement-induced frustration is demonstrated to be a powerful method for engineering complex textures in ferroics [40][41][42] .

Sample Fabrication
BFO thin films were deposited on DSO substrates at 700 °C using pulsed laser deposition with a laser energy density of ~1.5 J cm -2 in a background oxygen pressure of 90 mTorr.X-ray diffraction was measured on a Panalytical diffractometer with CuK source.Test structures were patterned using standard lithography techniques.
Ferroelectric domains were mapped using piezoresponse force microscopy in an asylum MFP-3D, both in-plane and out-of-plane, at two orthogonal sample rotations to vectorize the polarization.

NV microscopy
BFO samples were measured at room temperature using a commercial scanning NV microscope The magnetic field is quantitatively determined by measuring the optically detected magnetic resonance spectrum of the NV center at each point in space, yielding the projection of the magnetic field along the NV center axis.Large area images were collected in the qualitative "dual iso-B" mode, where the response of the NV center to two different microwave frequencies is used to track the magnetic field, as detailed in Ref 44 .This dramatically reduces data acquisition time with respect to collecting the full optically detected magnetic resonance spectrum.

Electronic Structure Calculation Details
DFT calculations were carried out in VASP.Hubbard U and Hund J corrections for the valence states of all species were calculated using collinear DFT using a linear response (LR) workflow developed in the atomate code framework 45 .Linear response analysis was performed without the inclusion of SOC to reduce computational cost, and because SOC has been found in other systems to have a relatively small effect on the on-site corrections from LR 45 .These values were calculated to be U, J = 5.2, 0.4 eV for Fe-d, U, J = 0.8, 0.8 eV for Bi-p, and U, J = 9.7, 1.9 eV for O-p.Because MAE is on the order of µeV, geometry relaxation, spin-orbit coupling (SOC), and on-site Hubbard U and Hund J corrections are all accounted for in a holistic manner to address their interdependence.In all calculations, Hubbard and Hund corrections are applied to all outer shell manifolds (Fe-d, Bi-p, and O-p).These structures were calculated for the structure endpoint reported in Ref. 18 .Using these calculated on-site Hubbard corrections, full geometry relaxation of the structure was performed with SOC included until self-consistency was reached for an electronic energy tolerance of 10 -6 eV.All calculations were performed for the 2x2x2 supercell (eight times the formula unit) to accommodate the G-type antiferromagnetic structure.These computational subtleties are addressed in greater depth in the Supp. Figure S4 and Note S1.

Figure 1 |
Figure 1 | Spin cycloid in BiFeO3.a Iron moments in BFO are antiferromagnetically aligned along the [111], modulated by the cycloid propagation along , [1 ̅ 10].Other allowed directions of  also lie within this (111).The canting of the AFM alignment gives rise to an uncompensated magnetization, (), which rotates primarily in the - plane with the same period as the antiferromagnetic moments, ~65 nm.b  is [110] to [1 ̅ 10] (horizontal and vertical in Figure 1c-f) with changes in .This behavior is due to the magnetoelectric coupling between  and  and demonstrates a nondegeneracy in the cycloid landscape.Generally, because  must remain orthogonal to , when  is along [1 ̅ 11],  can propagate along [110], but when  changes across a domain wall to [111],  must reorient to an allowed direction orthogonal to the new , either [1 ̅ 10], [1 ̅ 01], or [01 ̅ 1] (Figure

Figure 2 |
Figure 2 | Resolution of MSDW from first principles.a Schematic showing the initialization angle  within the cycloid plane for the antiferromagnetically aligned Fe spins.b Comparison of the magnetocrystalline anisotropy when the Fe spins are rotated in the (111) and (112 ̅ ).This agrees with the expectation that the cycloid rotates within (112 ̅ ), as the mean value of the energy is 2x higher when moving the rotation to the (111) plane.c,d Illustrations of the rotation planes in b. e Graphical representation of the magnetic cycloid directly from DFT, showing the relaxed Fe spins (red, blue), as well as the resulting canted moment (green),   , that spontaneously arises.f Components of the relaxed canted moment, showing the periodicity of MSDW which follows the cycloid period.The components sum to a vector of varying magnitude along [112 ̅ ] and   is largest when the Fe spins are pointed along [1 ̅ 10].g Relative energy along the three possible  directions when the unit cell is epitaxially strained to DSO.The mean energy is 2x lower when the cycloid can reorient the in-plane component of polarization in BFO to produce a net polarization (  ) along the [100] or [1 ̅ 00], composed of ([111 ̅ ] and [11 ̅ 1 ̅ ]) or ([1 ̅ 11 ̅ ] and [1 ̅ 1 ̅ 1 ̅ ]) polarized domains, respectively.Test structures (Methods) used to apply this in-plane field and the corresponding ferroelectric switching behavior are shown (Figure 3a,b).The ferroelectric domains measured by PFM are shown for devices, respectively poled into the   ∥ [010] (Figure 3c) and   ∥ [01 ̅ 0] (Figure 3d) configurations.Mapping the spin cycloid, measured through NV microscopy, to these domain images, the changes in  map directly to the ferroelectric domain walls and the sense of the cycloid such that the relationship  perpendicular to  is preserved.Here, for example, [111 ̅ ] and [1 ̅ 1 ̅ 1 ̅ ] polarized domains show equivalent  axes along [1 ̅ 10].

Figure 3 |
Figure 3 | In-plane switching of ferroelectric domains.a Schematic of the test structure used to apply an in-plane electric field along the [010] direction, showing the orientation of the ferroelectric domain walls perpendicular to the applied field.b Ferroelectric hysteresis loop of the device in a, where ferroelectric switching happens at -30/+40 V. c,d Ferroelectric domains overlayed with NV magnetometry data for devices poled both with net polarization left and right.Ferroelectric domain walls are shown as black lines

Figure 4 |c
Figure 4 | Out-of-plane switching of the ferroelectric domains.a,b Vector map constructed from the PFM of BiFeO3 samples in the (a) as-deposited and (b) out-of-plane poled configurations.The sample is poled using the PFM tip as the top electrode and +/-8V, corresponding to an electric field of ~2 MV cm -1 .c Map of the reorientation of ferroelectric domains, calculated from the difference of the polarization vectors in a and b. d Illustration showing the different ferroelectric switching events under the application of an out-of-plane electric field.e,h zoomed-in areas of the ferroelectric polarization, shown by the blue and purple boxes labeled 1 and 2 in b.The in-plane projection of the ferroelectric polarization is shown by the blue arrows.f,i Calculated switching events from c at the same locations as e,h.g,j NV microscopy images showing the spin cycloid at the same locations in e,h. and  are shown by blue and yellow arrows.

(
Qnami ProteusQ) which combines a confocal optical microscope with a tuning-fork based atomic force microscope.Diamond tips with a parabolic taper containing single NV centers were used to increase photon collection efficiency (Quantilever MX+).The orientation of the NV center in the )/MgO/CoFeB(0.9nm)/Ta(5nm)sample with perpendicular magnetic anisotropy.