# Influence of flexoelectricity on the spin cycloid in (110)-oriented BiFe O 3 films

D. Sando, F. Appert, S. Burns, Q. Zhang, Y. Gallais, A. Sacuto, M.

Cazayous, V. Garcia, S. Fusil, C. Carretero, et al.

► To cite this version:

D. Sando, F. Appert, S. Burns, Q. Zhang, Y. Gallais, et al., Influence of flexoelectricity on the spin cycloid in (110)-oriented BiFe O 3 films. Physical Review Materials, 2019, 3 (10), 10.1103/PhysRev-Materials.3.104404. [arXiv:2002.08230]

HAL Id: hal-02310230

https://normandie-univ.hal.science/hal-02310230v1

Submitted on 10 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Influence of flexoelectricity on the spin cycloid in (110)-oriented BiFeO₃ films

D. Sando,<sup>1,2,*</sup> F. Appert,<sup>3</sup> S. R. Burns,<sup>1</sup> Q. Zhang,<sup>1</sup> Y. Gallais,<sup>4</sup> A. Sacuto,<sup>4</sup> M. Cazayous,<sup>4</sup> V. Garcia,<sup>5</sup>
S. Fusil,<sup>5</sup> C. Carrétéro,<sup>5</sup> J. M. Le Breton,<sup>3</sup> A. Barthélémy,<sup>5</sup> M. Bibes,<sup>5</sup> J. Juraszek,<sup>3</sup> V. Nagarajan<sup>1</sup>

<sup>1</sup>School of Materials Science and Engineering, UNSW Sydney, High Street, Kensington, 2052, Australia

<sup>2</sup>Mark Wainwright Analytical Centre, UNSW Sydney, High Street, Kensington, 2052, Australia

<sup>3</sup>Normandie Univ., UNIROUEN, INSA Rouen, CNRS, GPM, 76000 Rouen, France

<sup>4</sup>Laboratoire Matériaux et Phénomènes Quantiques (UMR 7162 CNRS), Université Paris Diderot-Paris 7,
75205 Paris Cedex 13, France

<sup>5</sup>Unité Mixte de Physique, CNRS, Thales, Univ. Paris-Sud, Université Paris-Saclay, 91767 Palaiseau, France

daniel.sando@unsw.edu.au

**Abstract:**

The influence of film orientation, strain relaxation, and flexoelectric fields on the stability of the spin cycloid in (110)-oriented BiFeO₃ epitaxial films grown on LaAlO₃ substrates is investigated. By means of advanced x-ray diffraction techniques, we show that thinner films have very large strain gradients which give rise to high flexoelectric fields. Using low-energy Raman spectroscopy and conversion electron Mössbauer spectroscopy (CEMS) we show that films up to 53 nm thick possess collinear antiferromagnetic order, with no cycloidal modulation. This suppression of the cycloid is proposed to be from strain and strain-gradient-induced flexoelectric fields. On the other hand, films thicker than 90 nm show a complex spin texture consistent with two separate cycloids, likely with different propagation directions. Interestingly, CEMS analysis suggests that the two cycloids have the same spin rotation plane. The multiple cycloids are suggested to arise from different ferroelastic domains (in turn influenced by twinning in the substrate) with different strain relaxation behaviors. These results offer insight into the factors that influence cycloid stability in the less common (110) film orientation and have implications for future magnonic devices.
# 1. Introduction

BiFeO₃ (BFO) is a room temperature multiferroic [1]. Intense research interest into BFO, and other multiferroic materials, has been driven by the promise of devices in which energy consumption can be reduced by virtue of such materials to, in principle, allow the writing of a magnetic bit by electric field [2]. While progress is being made towards this goal, BFO has shown other qualities, such as conductive domain walls [3] (valuable for paradigm-changing data storage devices [4]), attractive optical [5–7] and photovoltaic [8] responses (good for optical and energy harvesting devices) as well as highly strain-tunable magnetic [9], ferroelectric (FE) [10–12] and multiferroic [13] characteristics.

From a magnetic point of view, bulk BFO has G-type antiferromagnetic (AFM) order [14] with a Néel temperature of ~640 K. The exact spin structure is more complex, however, with a long-period spin cycloid [15], typically confined to high-symmetry crystallographic directions [16], along with a weak ferromagnetic (FM) moment when the cycloid is suppressed (by high magnetic (H) field, strain, or otherwise). The spin cycloid in BFO gives rise to magnon (spin wave) modes [17] which can be probed by optical spectroscopy techniques [18,19]. Evidence for cross coupling between FE and magnetic orders has been evidenced by neutron diffraction [20] and low-energy Raman spectroscopy [21].

The influence of external perturbations (such as electric field, magnetic field, strain, etc.) on the cycloid in bulk BFO is well documented [20–25]. Furthermore, the modification of the unit cell by chemical strain (doping) can destabilize the cycloid [26], and hydrostatic pressure modulates the magnon modes through various structural phase transitions [27]. The influence of strain and external fields (electric *E*, magnetic *H*) on the cycloid is thus well understood, at least in the ‘standard’ (001) film orientation and under biaxial strain. On the other hand, and although the cycloid has been evidenced in thick (~1 µm) (110) films on SrTiO₃ substrates [28], the influence of uniaxial-like substrate-induced strain, and more importantly strain gradient (which can be accessed using strain relaxation in typically 20-100 nm thick partially-strained films) has received less attention.

Flexoelectricity refers to the generation of an electric field by a strain gradient. Although the magnitude of strain gradients in bulk materials is rather low (0.1 m⁻¹) [29], in nanoscale epitaxial films strain gradients can be huge (up to 10⁷ m⁻¹) [30] and thus generate enormous flexoelectric fields. In ferroelectric (FE) materials, these fields can influence self-polarization [31] or induce polarization rotation [32]. Using nanoscale mechanical pressure and the related localized strain gradient has been shown to switch FE polarization [33]. While the influence of flexoelectricity on FE
properties is thus well documented, there are relatively few reports discussing the influence of strain gradients on magnetic properties in *multiferroics*. Understanding the influence of strain, substrate orientation, and strain relaxation on the cycloid's existence has implications for the design of magnonic devices [34] in the future.

In this paper, we consider such issues by studying a set of BFO films grown epitaxially on (110)-oriented LaAlO3 (LAO) substrates. This substrate is chosen since the misfit strain of 4.5 % allows to induce large strain gradients through strain relaxation. We find that thinner films do not exhibit a cycloid, likely due to strain levels that are too high, while thicker films show a complex spin texture consistent with the existence of two cycloids, corresponding to different ferroelastic domains in the films. At intermediate thicknesses, although strain levels are low enough to allow a cycloid [9], the flexoelectric field destabilizes it.

## 2. Film Growth and Structural Characterization

Epitaxial films of BFO (19-144 nm in thickness) were grown by pulsed laser deposition (PLD) on single-crystal LAO (110) substrates. A frequency-tripled Nd:YAG laser ($\lambda$ = 355 nm, $f$ = 2.5 Hz, fluence ~1 J/cm$^2$) was focused on a $Bi_{1.1}^{57}FeO_3$ ceramic target (10 % excess of Bi to account for volatility of this element; ~100 % $^{57}$Fe substitution to facilitate CEMs measurements, as described in **Section 4**). The oxygen pressure was 0.01 mbar, and the substrate temperature was ~580 °C. Subsequent to film growth, the samples were cooled slowly to room temperature in an oxygen pressure of 300 mbar.

For films of thickness $t$ up to ~50 nm, the thickness was measured by x-ray reflectometry (XRR) (not shown). Using the calibration of the growth rate from thinner films, the thickness of the films > 60 nm was extrapolated from the growth time. As such, the film thicknesses reported here should be considered approximate.

X-ray diffraction (XRD) was performed to measure the structure of the films (Philips X'pert MRD diffractometer using Cu $K_{\alpha-1}$ radiation, equipped with a two-bounce Ge 220 monochromator). XRD $2\theta-\omega$ scans in the vicinity of the 110 and 220 LAO reflections are shown in **Fig. 1(a)**. In this manuscript we use monoclinic notation for BFO, for which the conversions to the LAO rhombohedral (≈ pseudocubic) frame are presented in **Table S1**, Supplemental Material<fnref n="1" />. Upon

<footnote n="1">¹ Unless denoted with 'mc', all crystallographic directions and reflections discussed in this paper are considered in the pseudocubic frame.</footnote>
increasing film thickness, the coexistence of several different peaks, with relative intensities varying with thickness, becomes apparent. We denote with a star (*), what appears to be for the 144 nm film, a tetragonal-like phase [13] BFO<sub>tetra</sub> with a characteristic elongated out-of-plane lattice parameter. The other three phases are, respectively, a ‘highly strained’ phase: BFO<sub>hstr</sub> (■); a phase with larger out-of-plane lattice parameter: BFO<sub>a</sub> (●); and phase with slightly lower out-of-plane lattice parameter: BFO<sub>b</sub> (▲). Such a phase mixture of various structures is typical in BFO films grown on LAO (110) substrates [35], and is likely due to the large (~4.5 %) misfit strain between substrate and film. Having seen the existence of multiple phases in thicker films we consider that the 19 nm film comprises a mixture of the BFO<sub>a</sub> and BFO<sub>b</sub> phases.

The out-of-plane lattice parameter $a_{BFO}$ for the three phases is plotted in **Fig. 1(b)**. From the XRD peak intensities (fit with Gaussian functions), we estimate the volume fraction of each phase, see **Fig. 1(c)**. In all the films the ‘tetragonal-like’ and ‘highly strained’ phases comprise an extremely low (0-2 %) volume fraction, therefore we neglect them for the remainder of the paper. For the other two phases, their volume fraction is complementary and BFO<sub>a</sub> becomes the dominant phase for thicknesses above ~90 nm.

<figure><img src="image_2.png" /></figure>

**Figure 1.** (Color online) (a) High angle 2θ-ω scans for BFO/LAO (110) films 19-144 nm in thickness. Multiple phases are identified (see text). Thickness dependence of (b) the out-of-plane lattice parameter, and (c) phase fraction of the different phases. (d) Average strain $ε_H$ for the BFO<sub>a</sub> phase.
More detailed characterization was carried out by performing XRD reciprocal space mapping (RSM) (Figs. S2 and S3, Supplemental Material). These analyses show that for the dominant phase BFOa there are two possible ferroelastic domains. Using the RSMs and 2θ-ω scans, we have extracted the full lattice parameters of BFOa. Here the average (homogeneous) strain εH (in percent) is calculated using the lattice parameter along the specific direction scaled down by the BFO unit cell volume. The strain values calculated using this formulation are presented in Fig. 1(d). This plot shows that the thinnest film (t = 19 nm) is under more than 2% compressive strain along the in-plane [001] direction and remarkably, experiences 1.2% tensile strain in the other in-plane [110] direction. Such an observation is rather counterintuitive since LAO nominally imposes 4.5% compressive strain. It is presumably the result of the strong uniaxial-type strain imposed by the LAO substrate along the [001] direction.

As the film thickness is increased, the magnitude of the in-plane strain decreases almost monotonically, as expected through strain relaxation [Fig. 1(d)]. It is interesting to note that the strain in the out-of-plane [110] direction (which is typically derived from the peak position in standard 2θ-ω scans) is close to constant. This observation indicates that using 2θ-ω scans alone to understand the relaxation behavior of (110)-oriented films is not sufficient, since strain relaxation occurs through the in-plane parameters, rather than manifesting as changes in the out-of-plane parameter.

To gain further insight into the structure, particularly in regard to strain relaxation, we have measured RSMs around a large number of peaks and performed a so-called Williamson-Hall (W-H) analysis [36]. Williamson Hall analysis allows to estimate the inhomogeneous strain εi (and therefore the strain gradient) in the out-of-plane and in-plane directions. The process described in Refs. [31,37,38] was followed, and a summary is outlined for a representative (t = 53 nm) sample in Fig. 2(a). The in-plane εi along [001] was estimated by fitting the film peaks for the 220, 221, and 222 reflections, while εi along [110] was estimated using the 110, 210, 220, and 320 reflections. Finally, the value of εi in the out-of-plane [110] direction was derived from peak fitting to the 110, 220, and 330 peaks from the 2θ-ω scan.

Using the relation

$$ \beta \cos \theta = \frac{K\lambda_w}{D} + 4\epsilon_i \sin \theta \quad (2) $$

where β = βfilm - βinstr is the linewidth of the diffraction peak [βinstr is the instrument resolution and is estimated from the line width of a nearby substrate peak], λw is the x-ray wavelength (1.5406 Å), D
is the coherence length along scattering vector, and K is a geometrical constant (taken here as 1), the value of $\epsilon_i$ is estimated separately for the three directions as described above.

Next, we use the value of $\epsilon_i$ in a formulation to extract the value of strain gradient in the films. To simplify the analysis, we assume an exponential strain profile, and that relaxation occurs predominantly through misfit dislocations. These assumptions allow us to define the film strain as a function of distance z from the substrate-film interface, and film thickness t, as [36]

$$ \epsilon(z, t) = \epsilon_0 \exp \left[ \cosh \frac{z}{\delta} - \tanh \frac{t}{\delta} \sinh \frac{z}{\delta} \right] \quad (3) $$

where $\epsilon_0$ is the strain at the substrate-film interface, and $\delta$ is the ‘strain penetration depth’; both of which are sample-specific. To incorporate the inhomogeneous strain, we use the treatment presented in Ref. [36], namely the relation

$$ \frac{t}{2\delta} \frac{\epsilon_H^2}{(\epsilon_i^2 + \epsilon_H^2)} = \tanh \left( \frac{t}{2\delta} \right) \quad (4) $$

where $\epsilon_H$ is the value of homogeneous (average) strain, as defined above. Numerically solving Eq. (4) allows determination of the value of $\delta$, and then substituting this into Eq. (3) yields the strain profile. Carrying out this analysis for the three directions ([001], [110], and [110]) allows to obtain the strain profile and thus the lattice constant as a function of z. The strain gradient is then calculated from the strain at the film surface and substrate-film interface (crudely as it approximates the strain profile as linear) by

$$ \sigma_{grad} = \frac{\epsilon(z=t,t) - \epsilon_0}{t} \quad (5) $$

The results of such analysis for the 53 nm film are shown in Fig. 2(b). Also plotted are the lattice parameters for the 19 nm thick BFO film, considered as an approximation to the expected lattice parameter near the interface for the 53 nm film. The agreement here is reasonable; however, the W-H analysis appears to overestimate (underestimate) the lattice parameter in the [001] direction ([110] and [110]) directions). This implies that the value of strain gradient for this film may in fact be larger than suggested by the present treatment.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Williamson-Hall plot analysis for strain gradients. (a) Williamson-Hall plot along the three cartesian directions. (b) predicted lattice constant (solid lines) according to Eq. (3). The average strain in each direction is shown as a dashed line. Data points are the average lattice parameter for the 19 nm thick film. (c) inhomogeneous strain in three directions as a function of thickness. (d) estimated flexoelectric field as a function of thickness (the dotted line is a guide to the eye).</figcaption></figure>

Applying the same treatment to the other films yields the value of inhomogeneous strain given in Fig. 2(c). Generally speaking, the inhomogeneous strain $ε_i$ decreases with thickness in the out-of-plane [110] and in-plane [1\overline{1}0] directions, while the $ε_i$ for the other in-plane [001] direction generally tends to increase upon increasing thickness. The 93 nm film, however, does not follow the trend as it has a very low value of $ε_i$ along the [001] direction. The origin for this is unclear; it may be related to some extrinsic effect related to film growth or substrate quality. Further, the 19 nm film has values of $ε_i$ that do not follow the trend; i.e., they seem to be too low. One can understand this when considering the lattice parameters and homogeneous strain $ε_h$ [Fig. 1(d)]. This film is under considerably higher *average strain* than the thicker films, thus implying that less strain relaxation has occurred; in turn meaning that the inhomogeneous strain is expected to be lower in magnitude.

Finally, the magnitude of the flexoelectric field estimated [31] from all the components of inhomogeneous strain is plotted in Fig. 2(d). The value of flexoelectric field estimated to be as high as 11 MV/m for the 19 nm film, and upon increasing thickness its value reaches a minimum for the 93 nm film, and then slightly increases to ~2 MV/m for the thicker films [Fig. 2(d)]. These very high
values of flexoelectric field almost certainly play a role in the self-polarization of the film (as discussed in Ref. [31]), and they become important later in **Section 5** when we consider the magnetic structure of the films.

### 3. Low-Energy Raman Spectroscopy

To understand how the previously described strain relaxation, strain gradients (and related flexoelectric fields), and ferroelectric/ferroelastic domain structure influence the cycloidal magnetic order in the films, we used low-energy Raman spectroscopy. This technique probes the spin wave modes corresponding to in-cycloidal-plane ('cyclon') and out-of-cycloidal-plane ('extracyclon') excitations. These different modes are selected by Raman spectroscopy using crossed and parallel polarizers respectively. Remember that Raman scattering probes excitations with a total wave-vector close to zero. For a first-order Raman process with only one excitation involved, the Raman spectrum represents the excitations of the compound at the center of the Brillouin zone (BZ). In a simple antiferromagnet, there are two magnon branches and two magnons are expected at the center of the BZ, the phase and the amplitude modes. The phason mode is gapped due to the magnetic anisotropy in the plane. The anisotropy out of the cycloidal plane induces 'folding' of the two branches at the center of the BZ and the higher energy modes (cyclon and extra-cyclon modes) observed in the Raman spectra. Several peaks are even split by higher harmonics of the spin cycloid. As first predicted by de Sousa and Moore [17] and then observed by Cazayous et al. [18], in the bulk, the cyclon mode energy follows a geometric progression with mode number with a slope given by the cyclon energy. On the other hand, the extracyclon modes follow a slightly different progression, due to the opening of a gap attributed to the pinning of the cycloid plane to the ferroelectric polarization. The spin excitations of BFO are fully described by a microscopic model including easy-axis anisotropy, and Dzyaloshinskii-Moriya interactions [40–42].

Low-energy Raman spectra (collected at room temperature using the 532 nm laser line [9,43]) are presented in **Fig. 3**. Note that the 19 nm film was not measured as the average strain was considered to be too large to sustain the cycloid; an assumption that is justified as we will see later when CEMs data are presented. The Raman spectrum for the 53 nm film [Fig. 3(a)] shows a single peak for both configurations. Such an observation is consistent with the presence of standard G-type antiferromagnetic order [44]; i.e. no cycloid. Moving to the 93 nm film [Fig. 3(b)], a series of peaks for each configuration is observed. Using the scaling laws from Refs. [17,18] we estimate the period of this cycloid to be about $\lambda = 72$ nm. This is larger than the bulk value of 62 nm, but is
comparable to previously-reported values for films [24,43,45]. The lengthening of the cycloid period is consistent with a strain-induced anisotropy [46].

<figure><img src="image_4.png" /><figcaption>**Figure 3.** (Color online) Low energy Raman spectroscopy for BFO//LAO (110) films of various thicknesses. (a) The 53 nm film shows a single peak for each configuration, consistent with G-type AFM order (i.e. no cycloid); (b) 93 nm film shows a cycloid with period ~72 nm; (c-d) the 126 nm and 144 nm films show a complex spectrum suggesting the coexistence of two cycloids.</figcaption></figure>

For the thicker films ($t \ge 120$ nm), the picture becomes rather intriguing. The Raman spectra for the 126 nm film [Fig. 3(c)] and the 144 nm film show a larger number of peaks. In each spectrum, one can distinguish two sets of peaks. This could reflect the presence of two separate cycloids in these samples, each set of peaks being associated to a cycloid. The period of both cycloids, estimated from the peak separation for each series, is close to 76 nm. In the bulk, the position of the peaks is a function of easy axis anisotropy *K* along the polarization axis. The shift between the two set of peaks ($\approx 2$ cm⁻¹) might be induced by the fact that the both cycloids do not experience the same magnetic anisotropy due to different propagation directions. In BFO bulk, this constant has been evaluated to 0.0035 meV (Ref. [41]). The shift of 2 cm⁻¹ is roughly proportional to an additional variation of *K* of 0.0005 meV.

Regardless of the presence of one ($t = 93$ nm) or two ($t = 126$ nm and $144$ nm) cycloids, the Raman data clearly show that the structural changes induced by increasing thickness influences the magnetic structure in these films. To try to understand if the two cycloids observed in the Raman spectra for the thicker films arise from two cycloids with different propagation directions (and thus two different cycloidal planes), we next used conversion electron Mössbauer spectroscopy (CEMS). An advantage of this technique is that it is sensitive to the *average* spin direction (i.e. the plane in
which the spins rotate), and supplementary measurements with the sample tilted in various orientations allows unambiguous determination of the cycloid plane.

## 4. Conversion Electron Mössbauer Spectroscopy (CEMS)

The ~100 % ⁵⁷Fe enrichment in the films allows conversion electron Mössbauer spectroscopy measurements to be performed even on the thinnest film of 19 nm with good signal-to-noise ratio. In the case of a magnetic cycloidal modulation of spins, the Fe nuclear energy levels are perturbed by combined electric quadrupolar and anisotropic magnetic hyperfine interactions, resulting in line broadening and asymmetry of the Mössbauer spectrum of BFO [47]. Therefore, the CEMS technique can be used to probe the existence or absence of a magnetic cycloid in BFO thin films [9,48,49].

We begin with CEMS data taken with the gamma ray incident normal to the film surface, as shown in Fig. 4(a). The CEMS spectrum taken at 295 K in this geometry for the 144 nm film is shown in Fig. 4(b). The spectrum exhibits a six-line Zeeman-split magnetic hyperfine pattern with isomer shift (δ = 0.41 mm.s⁻¹) and hyperfine field values (Bhf ≈ 48.8 T) typical for octahedrally-coordinated Fe³⁺. The first salient feature is that the spectrum shows a clear asymmetry, as evidenced by the difference in height of peaks 1 and 6 (and peaks 2 and 5). Such asymmetry of the experimental spectral data can, in principle, be reproduced in the calculated spectrum by a distribution of both orientations and hyperfine field values, correlated to quadrupole shift values [9].

Neutron diffraction experiments have shown that in (110)-oriented BFO films, the cycloid propagation direction is modified from the ‘bulk-like’ <110> family of directions, to a unique [112] direction [28,50]. In such a cycloid, the spins rotate in, and are confined to, the (110) plane. Based on these observations, we have attempted to fit [9] the CEMS spectrum for the 144 nm film using a cycloidal modulation of Fe spins inside such a (110) plane; however, as shown in Fig. 4(b), the fit is not satisfactory. Notably, the intensity of peaks 2 and 5 is too low. If the direction of the Fe spins makes an angle β with the direction of the incoming gamma ray, the line intensity ratio R₂₃ between lines 2 and 3 (or lines 5 and 4) of a magnetic hyperfine sextet can be expressed as

$$R_{23} = \frac{4 \sin^2 \beta}{(1+\cos^2 \beta)} \quad (6)$$

In normal incidence geometry, this expression implies that if all the spins lie in the plane of the film, β = 90° and thus R₂₃ = 4. On the other hand, if all the spins are oriented out of the film plane, then β → 0° and R₂₃ → 0. The fact that in our fit of the spectrum shown in Fig. 4(b) gives a ratio R₂₃ that
is too low indicates that the magnetic structure of the film is more complex than a simple harmonic cycloid, and that it comprises, according to Eq. (6), a greater proportion of spins aligned perpendicular to the direction of the incident gamma ray beam, i.e. in the film plane. This could be due either to i) a mixed spin arrangement comprised of a cycloidal order and collinear antiferromagnetic (AFM) order (with the spins confined to the film plane), or ii) a strong anharmonicity [51] of the cycloid, preferentially aligning the spins in a direction close to $β = 90°$, i.e. along the [001] direction of the film plane.

To test this first hypothesis, we next fit the CEMS spectrum using a model based on a mixed cycloidal / collinear antiferromagnetic (AFM) order. In this model, the cycloid is assumed to have a [112&#x0304;] propagation direction and (1&#x0304;10) plane, while the collinear AFM phase is confined to be in the film plane, as illustrated schematically in Fig. 4(c). Using this model, the fit is substantially improved as shown in Fig. 4(d). For the 144 nm film, the fit yields approximately 68 % cycloid and 32 % collinear AFM order. This demonstrates that about 1/3 of the Fe spins lie in the film plane and that a perfectly harmonic cycloid alone cannot reproduce the data.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (Color online) Mössbauer spectroscopy in normal incidence on BFO//LAO (110) films. (a) Geometry of the measurement with the gamma ray incident along $\vec{k}1$, that is, along the [11&#x0304;10] direction. (b) Normal incidence spectrum fit with a cycloid with (1&#x0304;10) cycloid plane. (c) Adding a collinear component improves the fit considerably. (d) Spectra and fits for the full set of samples.</figcaption></figure>
Applying the same analysis to the CEMS data recorded for the other films yields the fits presented in **Fig. 4(e)**. For all films the fits are very good, attesting to the appropriateness of the analysis presented here. Note that the line intensities of peaks 2 and 5 progressively increase upon decreasing film thickness, and that the asymmetry in the spectra disappears for the thinnest films. This trend reflects a progressive change in the spin structure with the film thickness. In the thickest film ($t = 144$ nm), the cycloid is most prevalent and the $R_{23}$ value is close to 2 (as for a random orientation of Fe spins), but for the 52 and 19 nm films, the line intensity ratio $R_{23}$ is close to 4 indicating that strain- and/or shape-induced anisotropy appears to force the spins in the plane of the film, and the cycloid is no longer present.

Further insight on the true spin direction and cycloid plane was gained by conducting CEMS measurements with the sample tilted relative to the incoming γ-rays. In the first tilted experiment [the geometry for which is shown in **Fig. 5(a)**], the γ-ray is incident along the $[0\bar{1}0]$ direction, which is tilted 45° from the $[1\bar{1}0]$ direction and 45° from the $(1\bar{1}0)$ plane. The second tilted experiment [**Fig. 5(b)**] has the γ-ray incident 45° from the $[1\bar{1}0]$ direction, in the $(1\bar{1}0)$ plane. This choice of tilted geometries allows better characterization of the Fe spin orientation, and simultaneous fitting of the spectra for the three geometries allows to considerably reduce the uncertainty.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (Color online) Mössbauer spectroscopy in inclined geometry. (a) geometry of the 1st tilted experiment where the γ-ray is incident along the [01̄0] direction, i.e. tilted 45° to the (11̄0) plane. (b) geometry of the 2nd tilted experiment where the γ-ray is incident at an angle 45° from the [110] direction in the (11̄0) plane. (c) Mössbauer spectra for the 1st tilted experiment for the five films. Note that for the thickest film, the spectrum resembles the standard type-1 cycloid, while for the thinnest film, the spectrum is consistent with spins orthogonal to the γ-ray direction (see text). (d) the spectra for 2nd tilted geometry show minimal dependence on thickness. The thickest film shows a spectrum very similar to that observed in normal incidence, while the thinnest film shows a symmetric spectrum, consistent with spins tilted ~45° from the γ-ray direction (see text). (e) geometry of the tilted Mössbauer experiments and the relevant angles. (f) angles β for the collinear AFM phase, obtained from the fits of the data shown in (c-d). For a pseudo-collinear phase with spins aligned collinear along [001], the expected values of these angles are shown as dashed lines.</figcaption></figure>

The CEMS spectra for the 1st tilted experiment and 2nd tilted experiment are presented in **Fig. 5(c-d)** respectively. Upon inspection of these spectra it is apparent that the thickness-dependent changes in the intensity of peaks 5 and 6 is much smaller than was observed for normal incidence, as shown in **Fig. 5**. Analysis of these spectra (once again considering a mixed cycloidal and collinear phase) indicates that the spins are indeed confined to the (11̄0) plane (for the thicker films) and that for thinner films the spins are directed almost completely along the [001] in-plane direction. Indeed, for the three used experimental geometries (normal incidence, and the two tilted experiments), the
orientation of the hyperfine field $B_{hf}$ of the collinear phase with respect to the incident gamma ray (as shown in **Fig. 5(e)**) was deduced from the fit of the Mössbauer spectra, with corresponding spin angles $\beta_1$, $\beta_2$, and $\beta_3$ of the collinear phase presented in **Fig. 5(f)**. Here for all films, we find that the $\beta_1$ and $\beta_2$ fitted values are close to 90° and that $\beta_3$ is close to 45°, evidencing that the spins of the collinear phase are, approximately, confined to the in-plane [001] direction. The latter result is consistent with the fact that the unit cell of BFO is in an $M_B$ monoclinic phase under about 1 % biaxial strain in the *ab* plane, thus causing the spins to point close to the [001] direction [48].

Using the analysis of mixed collinear / purely harmonic cycloid, fitting for all films gives a fraction of collinear phase that increases upon decreasing thickness. For the 19 nm and 53 nm films, the collinear phase constitutes 100 % of the CEMS spectrum fit, thus implying that there is *no* cycloid in these samples. For the 93 nm, 126 nm, and 144 nm films the collinear phase fraction decreases upon increasing film thickness.

Such an observation can be understood by considering the phase coexistence as shown in **Fig. 1(c)**. Recall that the $BFO_b$ phase constitutes close to 50 % of the film volume for thinner films ($t = 19$ and 53 nm), while as thickness is increased the fraction of $BFO_b$ decreases. Also recall that the structure of this phase has a strain of ~0.5 %, but, more importantly, it has a strongly distorted unit cell which prohibits the presence of the cycloid [52]. Such a strain value implies that the average direction of the Fe spins should be along [001], that is, in the plane of the film. Therefore, the observation of increased collinear AFM fraction with decreasing thickness is consistent with the increasing fraction of the $BFO_b$ phase.

**5. Discussion**

The Raman and Mössbauer data show that the cycloid is destroyed in films with thickness $t \le 53$ nm, while for films with $t \ge 93$ nm, the cycloid survives. Here we summarize our data and suggest, based on the detailed structural and strain gradient analysis of **Section 2**, the origin for the observed cycloidal order as a function of film thickness. A summary of the observations presented in this paper is given in **Fig. 6**.

For the thinnest film ($t = 19$ nm), the average biaxial strain in the *ab* monoclinic plane is above +1 % [Fig. 6(a)(i)], and the flexoelectric field is very large at ~11 MV/m. CEMS evidences antiferromagnetic order with no cycloidal modulation, consistent with the large strain value,
implying that the film is under too much strain and/or strain gradient to allow a cycloid. The strain-
induced anisotropy forces the spins to lie in the in-plane [001] direction.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** (Color online) (a) Dependence on thickness of (i) homogeneous strain $\epsilon_H$, (ii) flexoelectric field $E_{flexo}$, (iii) fraction of BFO$_b$ phase, and (iv) fraction of pseudo-collinear phase. Regions where the cycloid is stable are shown in orange, while blue denotes non-cycloidal (pseudo-collinear) states. (b) Dependence of collinear phase fraction on the phase fraction of BFO$_b$ phase, showing a monotonic trend. Note that extrapolation to 0% of BFO$_b$ suggests that the cycloidal phase has a prominent anharmonicity (see text).</figcaption></figure>

For the 53 nm film, although the average strain is +0.65 %, i.e. low enough to allow the cycloid [9], no evidence for this is shown in Raman and CEMS measurements. In addition, the CEMS spectra from measurements with the sample tilted suggest that the spins, just as for the 19 nm film, are oriented close to the [001] in-plane direction. Recall that the W-H analysis presented in **Section 2** suggests that the flexoelectric field in this film is quite high at ~4 MV/m (see Fig. 6(a)(ii)). Recently, it was predicted that external electric field of 6 MV/m applied to bulk BFO is sufficient to switch the cycloid rotation plane, while higher fields can create a 'plane cycloid' [25]. The influence of electric field on the cycloid has indeed been demonstrated experimentally, where the domain populations and cycloid period were changed with applied electric field in films of BFO (Ref. [53]). In the context of these predictions and experimental observations, we argue that in the 53 nm film, the relatively large flexoelectric field destabilizes the cycloid. Although the field of 6 MV/m appears to be rather modest to be able to exert such a strong influence (when compared to the theoretical predictions), one should keep in mind that generally theoretical approaches overestimate coercive fields in ferroelectrics and multiferroics by orders of magnitude [54,55].
For the 93 nm film, Raman data show the existence of a cycloid with a period $\lambda = 72$ nm while CEMS analysis suggests about 24 % cycloid and 76 % collinear order. Once again, the average strain and strain gradient are both low enough to expect the cycloid to be present.

For the thickest films ($t = 126$ nm and $144$ nm), the existence of a cycloidal modulation from Raman data, with a decrease of the cycloid proportion from CEMS, is consistent with the reduced average strain, and reduced $E_{flexo}$. The detection of two separate cycloids by Raman spectroscopy can be explained by existence of two ferroelastic domains as found using XRD.

This discussion is summarized in **Fig. 6(a)**, where the average strain $E_H$, flexoelectric field $E_{flexo}$, phase fraction of BFO$_b$, and the fraction of collinear phase (from Mössbauer analysis), are plotted as a function of thickness. In these plots, blue zones correspond to condition where the cycloid is destabilized, while orange regions indicate phase space where the type-2 cycloid [9,28] is stable.

Finally, to understand the decreasing fraction of collinear order upon increasing thickness, we have plotted the collinear phase fraction (deduced from Mössbauer) against the fraction of BFO$_b$, see **Fig. 6(b)**. A linear trend is observed, and, notably, extrapolation to the y-axis (i.e. if there was 0 % BFO$_b$ phase), suggests that the Mössbauer could indicate a small fraction of the collinear phase. Note that if error bars from the regression analysis are considered, the intercept could lie between the values of 2.5 and 15 %, the lower limit thereof being within the experimental uncertainty of our CEMS technique. If we consider that the collinear phase is attributed to the BFO$_b$ phase, this could suggest that the cycloid in these films has a prominent anharmonicity, i.e., in a phase pure film, a fit with mixed cycloid / AFM character could also be interpreted as cycloid with strong anharmonicity.

## 6. Conclusion

In summary, we have presented a comprehensive study of the structure, strain relaxation, and magnetic order of (110)-oriented BFO thin films on LAO substrates. Despite the large (~4.5 %) misfit strain, the films crystallize in the R-like phase, under levels of strain that are dependent on the thickness. For the thinnest films of 19 nm and 53 nm, Raman and CEMS data clearly show the absence of the cycloidal modulation, which we explain by the presence of large strain and high flexoelectric fields in these samples. For the thicker films of 93, 126, 144 nm, CEMS data show a cycloid fraction that increases with thickness, related to two factors: the first is that the average strain
is reducing, and secondly the fraction of BFO<sub>b</sub> phase (which is too distorted to allow the cycloid) decreases upon increasing thickness. On the other hand, Raman spectroscopy measurements reveal the existence of a single cycloid for the 93 nm film, while for the 126 nm and 144 nm films two separate cycloids, as shown by two sets of magnon peaks, are observed. These results show that strain gradients and flexoelectric fields in multiferroic films can be used to engineer magnetic order. The work presented here has implications for the design of magnonic devices, and bring insight to the rich interplay between strain, strain gradients, ferroelectric domain structure, and the magnetic order in multiferroic thin films.

## Acknowledgements

This work was supported by an Australian Nanotechnology Network (ANN) travel fellowship, and by the Australian Research Council through Discovery Grants. S.R.B would like to thank AINSE Limited for providing financial assistance (Award - PGRA) to enable this work. M.C. would like to acknowledge support from the French National Research Agency (ANR) through the DYMMOS project. F.A., J.M.L.B., and J.J. would like to acknowledge support from Region of Normandy and the European Regional Development Fund of Normandy (ERDF) in the frame of the MAGMA project.
# Supplemental Material

## Conversion between monoclinic $M_B$ and pseudocubic frames

Since we use the monoclinic $M_B$ unit cell when considering the BFO structure [56], we have provided a table converting between the monoclinic and pseudocubic axes in **Table S1**.

<table><caption>**Table ST1.** Conversion of the axes for monoclinic $M_B$ BFO and the rhombohedral LAO lattice (≈ pseudocubic lattice), and reflections used for the RSMs.</caption><thead><tr><th></th><th>Monoclinic $M_B$<br>BFO</th><th>LAO rhombohedral<br>≈ pseudocubic</th></tr></thead><tbody><tr><td rowspan="3">Axes;<br>directions</td><td>$a_{BFO} = [100]$</td><td>$[110]$</td></tr><tr><td>$b_{BFO} = [010]$</td><td>$[\bar{1}10]$</td></tr><tr><td>$c_{BFO} = [001]$</td><td>$[001]$</td></tr><tr><td rowspan="7">Reflections</td><td>200</td><td>110</td></tr><tr><td>400</td><td>220</td></tr><tr><td>600</td><td>330</td></tr><tr><td>401</td><td>221</td></tr><tr><td>402</td><td>222</td></tr><tr><td>310</td><td>210</td></tr><tr><td>420</td><td>310</td></tr><tr><td></td><td>510</td><td>320</td></tr></tbody></table>

## Characterization of strain state and ferroelastic domain structure

Information about the in-plane lattice parameters of the films was obtained using XRD reciprocal space mapping (RSM). Symmetric and asymmetric reflections were measured to see the full picture of the film structure. In **Fig. S1**, RSMs of reflections for which the x-ray plane of diffraction is parallel with the $[\bar{1}10]$ direction are presented. With this azimuthal ($\phi$) angle, we access the symmetric 220 LAO reflection, as well as the asymmetric 310 LAO reflection<fnref n="2" />. The symmetric RSMs [Fig. S1(a-e)] reveal the existence of $BFO_a$ and $BFO_b$ for all film thicknesses, consistent with the results of standard 2θ-ω scans of **Fig. 1(a)**. Note that we do not observe peaks corresponding to $BFO_{hstr}$ phases in these RSMs. The maps in the vicinity of the asymmetric 310 LAO (420 BFO) reflection reveal two peaks, corresponding to $BFO_a$ and $BFO_b$; however, from these maps, neither of

<footnote n="2">² These reflections are thus called because for the symmetric reflections the XRD diffractometer angle $\omega \approx 2\theta/2$, while for the asymmetric reflections the scattering vector is required to be tilted such that $\omega \neq 2\theta/2$.</footnote>
these phases appears to be ‘fully-strained’, since the in-plane reciprocal lattice component for film and substrate are not equal. This reflects strain relaxation of the film layer. The lattice parameter *b*<sub>BFO</sub> is obtained for the two phases using these maps.

<figure><img src="image_8.png" /><figcaption>**Figure S1.** Reciprocal space mapping with φ = 0, scattering plane parallel with the [110] direction. (a-e) RSMs around the symmetric (220) reflection reveal the existence of a strained (BFO<sub>a</sub>) and relaxed (BFO<sub>b</sub>) phase for all films; (f-j) around the asymmetric (310) reflection show that all films are relaxed along the [110] in plane direction.</figcaption></figure>

Changing the azimuthal angle to φ = 90° in the XRD experiment gives us access once again to the symmetric 220 LAO reflection, and to the 221 LAO reflection. The symmetric reflection [Fig. S2(a-e)] once again reveals the presence of the BFO<sub>a</sub> and BFO<sub>b</sub> phases for all thicknesses, but also shows the ‘highly strained’ BFO<sub>hstr</sub> phases in the 93, 126 and 144 nm films. A salient feature of the symmetric maps is the splitting of the BFO<sub>a</sub> peak for the 126 nm film. This splitting implies that there are ‘tilted domains’ in this BFO phase. Such a behavior has been observed for thick films in the (110) orientation, and was attributed to strain relaxation [56].

The asymmetric maps [Fig. S2(f-j)] reveal the presence of two monoclinic ferroelastic domains, evidenced by the characteristic splitting [39] of the BFO<sub>a</sub> peak. One peak corresponds to the 401<sub>mc</sub> reflection of BFO<sub>a</sub>, while the other corresponds to the 401<sub>mc</sub> reflection. Interestingly, the ‘relaxed’ phase BFO<sub>b</sub> is not visible in these maps: this may be due to the peak being obscured by the peaks of the monoclinically-split BFO<sub>a</sub>. The domain population for BFO<sub>a</sub> appears not to follow a
particular trend; that is, for the 53 nm [Fig. S2(g)], 126 nm [Fig. S2(i)] and 144 nm [Fig. S2(j)] films, the lower peak (401<sub>mc</sub>) is stronger, while for the 93 nm [Fig. S2(h)] film it is the upper peak (401<sub>mc</sub>) that is more intense. The origin for this phenomenon is not entirely clear, it may be related to slight differences in growth conditions, differences in the miscut angle of the substrates, or variations in the proportion of the ferroelastic twins in the LAO substrate.

<figure><img src="image_9.png" /><figcaption>**Figure S2. Reciprocal space mapping with φ = 90°, scattering plane parallel with the (001) direction.** (a-e) RSMs around the symmetric (220) reflection reveal the existence of BFO<sub>a</sub> and BFO<sub>b</sub> phases for all films, and some tilted domains for the 126 nm thick film; (f-j) Maps around the asymmetric (221) reflection show that all films comprise two ferroelastic domains with proportions that vary with thickness.</figcaption></figure>

It is convenient to define an average homogeneous strain in the *ab* plane ε<sub>H ab plane</sub> (that is, the *ab* plane of the monoclinic BFO unit cell), calculated using the *a*<sub>BFO</sub> and *b*<sub>BFO</sub> parameters using

$$ \epsilon_H \text{ ab plane} = \left[ \frac{\sqrt{\frac{a_{BFO}b_{BFO}}{2} - \frac{3}{\sqrt{V_{pc}}}}}{\frac{3}{\sqrt{V_{pc}}}} \right] \times 100 \%, \quad (1) $$

where the lattice parameters for BFO are defined in the monoclinic M<sub>B</sub> frame (Table ST1), and V<sub>pc</sub> is the volume of the pseudocubic unit cell given by V<sub>pc</sub> = (a<sub>BFO</sub>×b<sub>BFO</sub>×c<sub>BFO</sub>) / 2. Using this strain value allows to compare the strain experienced by the present films with the more conventional (001) orientation. In Fig. 5(a) are plotted the values of ε<sub>H ab plane</sub> as a function of film thickness. The
thinnest film is under an average of +1 % strain, while in the thickest film the average strain reduces to +0.4 %. Note that this is tensile strain in the monoclinic $M_B$), and in lattice parameter, as (001) films grown for example on $SmScO_3$ or $NdScO_3$ (Refs. [5,9]).

Experimental details for conversion electron Mössbauer spectrometry (CEMS)

The CEMS experiments were performed at room temperature using a ~50 mCi source of $^{57}$Co in Rh matrix in constant acceleration mode. The conversion electrons were collected with a home-made gas flow proportional counter operating with a (He-5% $CH_4$) gas mixture [57]. A collimator consisting of a lead screen with a 3 mm diameter hole was placed in front of the detector window. The data were fitted using the MOSFIT code using the histogram method [58]. Isomer shift values are given with respect to α-Fe at 300 K.
References

[1] G. Catalan and J. F. Scott, Adv. Mater. **21**, 2463 (2009).

[2] M. Bibes and A. Barthélémy, Nat. Mater. **7**, 425 (2008).

[3] J. Seidel, L. W. Martin, Q. He, Q. Zhan, Y.-H. Chu, A. Rother, M. E. Hawkridge, P. Maksymovych, P. Yu, M. Gajek, N. Balke, S. V Kalinin, S. Gemming, F. Wang, G. Catalan, J. F. Scott, N. A. Spaldin, J. Orenstein, and R. Ramesh, Nat. Mater. **8**, 229 (2009).

[4] P. Sharma, Q. Zhang, D. Sando, C. H. Lei, Y. Liu, J. Li, V. Nagarajan, and J. Seidel, Sci. Adv. **3**, (2017).

[5] D. Sando, Y. Yang, E. Bousquet, C. Carrétéro, V. Garcia, S. Fusil, D. Dolfi, A. Barthélémy, P. Ghosez, L. Bellaiche, and M. Bibes, Nat. Commun. **7**, 10718 (2016).

[6] D. Sando, P. Hermet, J. Allibe, J. Bourderionnet, S. Fusil, C. Carrétéro, E. Jacquet, J.-C. Mage, D. Dolfi, A. Barthélémy, P. Ghosez, and M. Bibes, Phys. Rev. B **89**, 195106 (2014).

[7] D. Sando, Y. Yang, C. Paillard, B. Dkhil, L. Bellaiche, and V. Nagarajan, Appl. Phys. Rev. **5**, 041108 (2018).

[8] T. Choi, S. Lee, Y. J. Choi, V. Kiryukhin, and S.-W. Cheong, Science. **324**, 63 (2009).

[9] D. Sando, A. Agbelele, D. Rahmedov, J. Liu, P. Rovillain, C. Toulouse, I. C. Infante, A. P. Pyatakov, S. Fusil, E. Jacquet, C. Carrétéro, C. Deranlot, S. Lisenkov, D. Wang, J.-M. M. Le Breton, M. Cazayous, A. Sacuto, J. Juraszek, A. K. Zvezdin, L. Bellaiche, B. Dkhil, A. Barthélémy, and M. Bibes, Nat. Mater. **12**, 641 (2013).

[10] D. Sando, A. Agbelele, C. Daumont, D. Rahmedov, W. Ren, I. C. Infante, S. Lisenkov, S. Prosandeev, S. Fusil, E. Jacquet, C. Carretero, S. Petit, M. Cazayous, J. Juraszek, J.-M. Le Breton, L. Bellaiche, B. Dkhil, A. Barthelemy, and M. Bibes, Philos. Trans. R. Soc. A **372**, 20120438 (2014).

[11] I. C. Infante, S. Lisenkov, B. Dupé, M. Bibes, S. Fusil, E. Jacquet, G. Geneste, S. Petit, A. Courtial, J. Juraszek, L. Bellaiche, A. Barthélémy, and B. Dkhil, Phys. Rev. Lett. **105**, 057601 (2010).

[12] C. Daumont, W. Ren, I. C. Infante, S. Lisenkov, J. Allibe, C. Carrétéro, S. Fusil, E. Jacquet, T. Bouvet, F. Bouamrane, S. Prosandeev, G. Geneste, B. Dkhil, L. Bellaiche, A. Barthélémy, and M. Bibes, J. Phys. Condens. Matter **24**, 162202 (2012).
[13] D. Sando, B. Xu, L. Bellaiche, and V. Nagarajan, Appl. Phys. Rev. **3**, 011106 (2016).

[14] S. V. Kiselev, R. P. Ozerov, and G. S. Zhdanov, Sov. Phys. Dokl. **7**, 742 (1963).

[15] I. Sosnowska, T. Peterlin-Neumaier, and E. Steichele, J. Phys. C Solid State Phys. **15**, 4835 (1982).

[16] J.-G. Park, M. D. Le, J. Jeong, and S. Lee, J. Phys. Condens. Matter **26**, 433202 (2014).

[17] R. de Sousa and J. E. Moore, Phys. Rev. B **77**, 012406 (2008).

[18] M. Cazayous, Y. Gallais, A. Sacuto, R. de Sousa, D. Lebeugle, and D. Colson, Phys. Rev. Lett. **101**, 037601 (2008).

[19] U. Nagel, R. S. Fishman, T. Katuwal, H. Engelkamp, D. Talbayev, H. T. Yi, S.-W. Cheong, and T. Rõõm, Phys. Rev. Lett. **110**, 257201 (2013).

[20] D. Lebeugle, D. Colson, A. Forget, M. Viret, A. M. Bataille, and A. Gukasov, Phys. Rev. Lett. **100**, 227602 (2008).

[21] P. Rovillain, R. De Sousa, Y. Gallais, A. Sacuto, M. Méasson, D. Colson, A. Forget, M. Bibes, A. Barthélémy, and M. Cazayous, Nat. Mater. **9**, 975 (2010).

[22] Y. F. Popov, A. K. Zvezdin, G. P. Vorob’ev, A. M. Kadomtseva, V. A. Murashov, and D. N. Rakov, JETP Lett. **57**, 69 (1993).

[23] M. Tokunaga, M. Azuma, and Y. Shimakawa, J. Phys. Soc. Japan **79**, 064713 (2010).

[24] I. Gross, W. Akhtar, V. Garcia, L. J. Martínez, S. Chouaieb, K. Garcia, C. Carrétéro, A. Barthélémy, P. Appel, P. Maletinsky, J.-V. Kim, J. Y. Chauleau, N. Jaouen, M. Viret, M. Bibes, S. Fusil, and V. Jacques, Nature **549**, 252 (2017).

[25] A. F. Popkov, N. E. Kulagin, S. V. Soloviov, K. S. Sukmanova, Z. V. Gareeva, and A. K. Zvezdin, Phys. Rev. B **92**, 140414(R) (2015).

[26] I. Sosnowska, W. Schäfer, W. Kockelmann, K. H. Andersen, and I. O. Troyanchuk, Appl. Phys. A **74**, S1040 (2002).

[27] J. Buhot, C. Toulouse, Y. Gallais, A. Sacuto, R. de Sousa, D. Wang, L. Bellaiche, M. Bibes, A. Barthélémy, A. Forget, D. Colson, M. Cazayous, and M.-A. Measson, Phys. Rev. Lett. **115**, 267204 (2015).

[28] W. Ratcliff, D. Kan, W. Chen, S. Watson, S. Chi, R. Erwin, G. J. McIntyre, S. C. Capelli, and
I. Takeuchi, Adv. Funct. Mater. **21**, 1567 (2011).

[29] P. Zubko, G. Catalan, and A. K. Tagantsev, Annu. Rev. Mater. Res. **43**, 387 (2013).

[30] D. Lee and T. W. Noh, Philos Trans A Math Phys Eng Sci **370**, 4944 (2012).

[31] B. C. Jeon, D. Lee, M. H. Lee, S. M. Yang, S. C. Chae, T. K. Song, S. D. Bu, J.-S. Chung, J.-G. Yoon, and T. W. Noh, Adv. Mater. **25**, 5643 (2013).

[32] G. Catalan, A. Lubk, A. H. G. Vlooswijk, E. Snoeck, C. Magen, A. Janssens, G. Rispens, G. Rijnders, D. H. A. Blank, and B. Noheda, Nat. Mater. **10**, 963 (2011).

[33] H. Lu, C.-W. Bark, D. Esque de los Ojos, J. Alcala, C. B. Eom, G. Catalan, and A. Gruverman, Science. **336**, 59 (2012).

[34] V. V Kruglyak, S. O. Demokritov, and D. Grundler, J. Phys. D. Appl. Phys. **43**, 264001 (2010).

[35] H. Liu, P. Yang, L. You, Y. Zhou, Z. Fan, H. R. Tan, J. Wang, J. Wang, and K. Yao, J. Appl. Phys. **118**, 104103 (2015).

[36] G. Catalan, B. Noheda, J. McAneney, L. J. Sinnamon, and J. M. Gregg, Phys. Rev. B **72**, 020102 (2005).

[37] D. Lee, A. Yoon, S. Y. Jang, J. G. Yoon, J. S. Chung, M. Kim, J. F. Scott, and T. W. Noh, Phys. Rev. Lett. **107**, 057602 (2011).

[38] D. Sando, C. Carrétéro, M. N. Grisolia, A. Barthélémy, V. Nagarajan, and M. Bibes, Adv. Opt. Mater. **6**, 1700836 (2017).

[39] G. Xu, J. Li, and D. Viehland, Appl. Phys. Lett. **89**, 222901 (2006).

[40] R. S. Fishman, N. Furukawa, J. T. Haraldsen, M. Matsuda, and S. Miyahara, Phys. Rev. B **86**, 220402 (2012).

[41] R. S. Fishman, J. T. Haraldsen, N. Furukawa, and S. Miyahara, Phys. Rev. B **87**, 134416 (2013).

[42] R. S. Fishman, J. H. Lee, S. Bordács, I. Kézsmárki, U. Nagel, and T. Rõõm, Phys. Rev. B **92**, 094422 (2015).

[43] A. Agbelele, D. Sando, C. Toulouse, C. Paillard, R. D. Johnson, R. Rüffer, A. F. Popkov, C. Carrétéro, P. Rovillain, J. M. Le Breton, B. Dkhil, M. Cazayous, Y. Gallais, M. A. Méasson,
A. Sacuto, P. Manuel, A. K. Zvezdin, A. Barthélémy, J. Juraszek, and M. Bibes, Adv. Mater. **29**, 1602327 (2017).

[44] D. Wang, J. Weerasinghe, and L. Bellaiche, Phys. Rev. Lett. **109**, 067203 (2012).

[45] W. Saenrang, B. A. Davidson, F. Maccherozzi, J. P. Podkaminer, J. Irwin, R. D. Johnson, J. W. Freeland, J. Íñiguez, J. L. Schad, K. Reierson, J. C. Frederick, C. A. F. Vaz, L. Howald, T. H. Kim, S. Ryu, M. Veenendaal, P. G. Radaelli, S. S. Dhesi, M. S. Rzchowski, and C. B. Eom, Nat. Commun. **8**, 1583 (2017).

[46] Z. V. Gareeva, A. F. Popkov, S. V. Soloviov, and A. K. Zvezdin, Phys. Rev. B **87**, 214413 (2013).

[47] D. Lebeugle, D. Colson, A. Forget, M. Viret, P. Bonville, J. F. Marucco, and S. Fusil, Phys. Rev. B **76**, 024116 (2007).

[48] A. Agbelele, D. Sando, I. C. Infante, C. Carrétéro, S. Jouen, J. M. Le Breton, A. Barthélémy, B. Dkhil, M. Bibes, and J. Juraszek, Appl. Phys. Lett. **109**, 042902 (2016).

[49] H. Hojo, R. Kawabe, K. Shimizu, H. Yamamoto, K. Mibu, K. Samanta, T. Saha-Dasgupta, and M. Azuma, Adv. Mater. **29**, 1603131 (2017).

[50] S. R. Burns, D. Sando, B. Xu, B. Dupé, L. Russell, G. Deng, R. Clements, O. H. C. Paull, J. Seidel, L. Bellaiche, N. Valanoor, and C. Ulrich, NPJ Quantum Mater. **4**, 18 (2019).

[51] M. Ramazanoglu, W. Ratcliff II, Y. J. Choi, S. Lee, S.-W. Cheong, and V. Kiryukhin, Phys. Rev. B **83**, 174434 (2011).

[52] D. Sando, F. Appert, B. Xu, O. H. C. Paull, S. R. Burns, C. Carrétéro, B. Dupé, V. Garcia, Y. Gallais, A. Sacuto, M. Cazayous, B. Dkhil, J.-M. Le Breton, A. Barthélémy, M. Bibes, L. Bellaiche, V. Nagarajan, and J. Juraszek, Submitt. to Appl. Phys. Rev. (2019).

[53] W. Ratcliff, Z. Yamani, V. Anbusathaiah, T. R. Gao, P. A. Kienzle, H. Cao, and I. Takeuchi, Phys. Rev. B **87**, 140405 (2013).

[54] S. Lisenkov, D. Rahmedov, and L. Bellaiche, Phys. Rev. Lett. **103**, 047204 (2009).

[55] R. Landauer, J. Appl. Phys. **28**, 227 (1957).

[56] D. Kan and I. Takeuchi, J. Appl. Phys. **108**, 014104 (2010).

[57] J. Juraszek, O. Zivotsky, H. Chiron, C. Vaudolon, and J. Teillet, Rev. Sci. Instrum. **80**, 043905 (2009).
(2009).

[58] J. Teillet, F. Varret, and J. Juraszek, MOSFIT Program, Unpubl. (n.d.).