Expansion of the spin cycloid in multiferroic BiFeO3 thin films

Understanding and manipulating complex spin texture in multiferroics can offer new perspectives for electric field-controlled spin manipulation. In BiFeO3, a well-known room temperature multiferroic, the competition between various exchange interactions manifests itself as non-collinear spin order, i.e., an incommensurate spin cycloid with period 64 nm. We report on the stability and systematic expansion of the length of the spin cycloid in (110)-oriented epitaxial Co-doped BiFeO3 thin films. Neutron diffraction shows (i) this cycloid, despite its partly out-of-plane canted propagation vector, can be stabilized in thinnest films; (ii) the cycloid length expands significantly with decreasing film thickness; (iii) theory confirms a unique [112¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$11\bar 2$$\end{document}] cycloid propagation direction; and (iv) in the temperature dependence the cycloid length expands significantly close to TN. These observations are supported by Monte Carlo simulations based on a first-principles effective Hamiltonian method. Our results therefore offer new opportunities for nanoscale magnonic devices based on complex spin textures.


INTRODUCTION
Understanding spin-lattice interactions in magnetic oxide materials is critical to achieve emergent phenomena in such systems. For example magneto-transport, a property crucial for advancing spinbased technologies, can be affected by magneto-striction, which relies on the intimate coupling between lattice and spin. [1][2][3] When confined to the nanoscale, these interactions are further enhanced by the imposed size constraint and inherent symmetry breaking at interfaces. 4,5 In this context, the last few years have seen intense efforts to harness the spin degree of freedom in bismuth ferrite BiFeO 3 (BFO), a room temperature multiferroic, to yield exotic magnetic arrangements at the nanoscale. [6][7][8][9][10][11] In the bulk, BFO crystallizes in the R3c space group and possesses G-type antiferromagnetic order with a magnetic phase transition at T N = 643 K 12-16 and a ferroelectric transition at T C = 1123 K. [14][15][16][17] Additional Dzyaloshinskii-Moriya interactions (DMI) establish both spin canting and a long-period spin cycloid of length~64 nm. [14][15][16] The propagation direction of the spin cycloid can be modified by switching the ferroelectric polarization 18 thus demonstrating prominent magnetoelectric coupling. Upon destruction of the cycloid-achievable by magnetic [19][20][21][22] or electric 23 field, strain, 24 doping, 25 or hydrostatic pressure 26 -a weak canted ferromagnetic moment is released. 27,28 To date, the majority of understanding has focused on the manipulation and engineering of this weak canted moment in thin films. 4 On the other hand, the presence of a spin cycloid gives rise to specific magnon modes. [29][30][31][32] Magnons are collective spin-wave excitations and show a sensitivity to applied electric field 33 or hydrostatic pressure, 26 permitting their exploitation for use in magnonic devices in which spin waves are used to carry and process information. 34,35 In this regard, epitaxial thin films present a model platform to understand how the spin cycloid-and its related magnonic response-reacts to controlled boundary conditions such as strain or electric field. Further, by virtue of their quasi-two-dimensional morphology, epitaxial thin films lend themselves to a systematic understanding of finite-size effects, including their influence on the exchange interactions responsible for the cycloid, particularly near critical points.
Here, we investigate the stability and systematic variation of the spin cycloid length in the temperature and film-thickness parameter space in a series of 2% cobalt substituted BFO thin films using neutron diffraction and first-principles-based calculations. Co-doping in addition to strain engineering along with careful choice of substrate orientation, allows us to stabilize a cycloidal order (with a partly out-of-plane canted propagation vector) in films as thin as 25 nm, i.e., smaller than length of the spin cycloid. Theory confirms a new [112] cycloid propagation direction. Crucially, we find that the cycloid period increases upon decreasing film thickness, and also as the sample is heated to the magnetic phase transition.

Growth and characterization of thin films
Epitaxial films of 2% Co-doped BiFeO 3 , 25-192 nm in thickness, were grown by pulsed laser deposition (PLD) on (110)-oriented SrTiO 3 (STO) substrates (pseudocubic notation is used throughout). The (110) orientation was chosen as it enables one to impose a moderate uniaxial compressive strain which does not destroy the cycloid in the thinnest BFO films. 10,24 High-angle X-raydiffraction (XRD) θ − 2θ scans (see Supplementary Fig. 1) indicate phase-pure (110)-oriented growth, while XRD reciprocal-space mapping (RSM) shows that the films are strained along the inplane [001] direction ( Supplementary Fig. 1). The 50, 142, and 192 nm thick films comprise a single ferroelastic domain, while the 25, 65, and 90 nm films are found to consist of two domains, one of which constitutes approximately 3:1 to 4:1 volume fraction (see Supplementary Information). X-ray reflection measurements indicated that the surface and interface roughnesses are less than 2 nm for all films. Furthermore, there is no evidence for an interfacial or surface layer with different composition (see Supplementary Fig. 2). This further demonstrates the excellent quality of the grown films.

Spin cycloid in BFO thin films
The precise structure of the spin cycloid was determined by neutron diffraction experiments using a triple-axis spectrometer (see Methods). Figure 1 presents neutron diffraction and XRD data for the 50 and 192 nm films. Neutron RSMs were measured around the ð 1 2 1 2 1 2 Þ and the inverse ð 1 2 1 2 1 2 Þ G-type antiferromagnetic Bragg peak positions (see Fig. 1a, e). The sample alignment was performed with respect to the lattice parameters of the STO substrate. Due to the epitaxial nature of the PLD grown films, the position of the magnetic signal of the BFO film is matched in the in-plane Q L direction and shifted in the out-of-plane Q HK direction. This is in agreement with the X-ray RSM as shown in Fig. 1d, h, respectively. The X-ray maps further confirm that both thin films comprise a single monoclinic domain. In addition to the magnetic signal in the neutron RSMs arising from the BFO film, a peak at ð 1 2 1 2 1 2 Þ is visible. This peak is caused by a second-order contamination originating from the strong structural (111) reflection of the STO substrate. The elliptical shape of the Bragg peaks is due to the finite instrumental resolution, which is different for different directions in Q-space.
For both films the neutron diffraction data show a peak splitting around one of the two equivalent G-type antiferromagnetic Bragg peak positions, i.e., around ð 1 2 1 2 1 2 Þ. This is confirmed in the neutron diffraction data along the ½112 direction across the magnetic Bragg peaks, as shown in Fig. 1c (50 nm) and 1g (192 nm). In contrast, the complementary reflections ð 1 2 1 2 1 2 Þ, displayed in the insets of Fig. 1a, e, show a single peak, consistent with cycloidallysplit peaks that are superimposed on each other in this Q direction due to instrumental resolution (see Supplementary Fig. 3). The incommensurate splitting of the peaks can be attributed to a cycloidal spin arrangement. This was demonstrated previously by polarized neutron diffraction experiments on thin-film BFO 2 Þ data a weak second-order contamination (termed λ/2) originating from the substrate, is also present. All neutron diffraction data were taken at a temperature of 150 K S.R. Burns et al.
samples. 10,36 Following the procedure outlined in refs. 6,10 , we determine that the splitting of the peaks occurs along the ½112 direction and not ½110 as in bulk BFO. The magnitude of the incommensurable splitting is δ = 0.0043 r.l.u. for the 192 nm sample, allowing to deduce a cycloid period of λ 0 = 65.1(2.5) nm consistent with previous reports on bulk BFO 12,[14][15][16]18 or 100-200 nm thin BFO films. 6,10,37 On the other hand, for the 50 nm film, the splitting is less pronounced. Possible explanations for such weaker splitting include (i) a longer cycloid period, yielding peaks that are closer to each other in reciprocal space; or (ii) the coexistence of different domains, one of which has G-type antiferromagnetic order (one peak) and the other which has cycloidal modulation (two peaks), resulting in a superposition of three peaks. However, since the film is shown to have a single ferroelastic domain, we can rule out case (ii). Furthermore, Fig. 1c indicates that the neutron scattering data can be fitted with two well-separated magnetic Bragg peaks with the linewidth corresponding to the experimental resolution where no additional peak in-between is required. The occurrence of the peak splitting around ð 1 2 1 2 1 2 Þ and the lack thereof in the complementary Q L direction confirms a single cycloid with a propagation vector along ½112 for both films.
To verify that the reduced peak splitting of the 50 nm film vis-àvis the 192 nm film is possible from a change in the cycloid period (case (i) above), we have simulated the neutron diffraction pattern using the cycloid period found in bulk (64 nm) and a longer period of 84 nm (for details see Supplementary Information). For the 192 nm film, the simulated pattern for a cycloid period corresponding to the bulk agrees well with the measured neutron data (see Fig. 1f). On the other hand, for the 50 nm film, a period of 84 nm shows the most convincing agreement to the data (Fig. 1b). Line profiles along the ½112 direction for each of the reflections around ð 1  12,[14][15][16]18 As will be shown later, there is a systematic increase of the period length with decreasing film thickness. Since the strain state of the films is virtually thickness independent in our films (see Supplementary Fig. 1), in contrast with previous studies on (001) oriented films, 37 the data suggest that the spin cycloid is sensitive to finite size or surface and interfacial effects. Figure 1 thus presents an exciting finding: although it is possible, in the (001) film orientation, to stabilize a cycloid with propagation vector in the film plane for thicknesses less than the period, 38 here we demonstrate that this phenomenon is indeed possible for cycloid propagation vectors with a component out of the film plane. To rationalize our interpretation of the magnetic structure found by neutron diffraction, XRD RSMs of the (221) reflection were acquired (Fig. 1d, h). These reflections were chosen since they are in the same scattering geometry as the neutron RSMs, and crucially, a direct comparison between neutron-and Xray-diffraction data can be made. For both samples, a single XRD film peak is observed, reflecting a single ferroelastic domain state in these films. This implies unambiguously that the magnetic splitting observed in Fig. 1a, e arises from cycloidal order with a single propagation direction, namely ½112.
Expansion of the spin cycloid in thinnest films To investigate the length variation of the spin cycloid period in more detail, we have measured 2% Co-doped BFO films with thicknesses from 25 to 192 nm. Figure 2 presents neutron diffraction data for films of thickness 65, 90, and 142 nm. The analysis in case of the 25, 65, and 90 nm films is more complex due to the multiple domains present in these films (see Supplementary Fig. 1). The neutron diffraction RSM for the 65 nm film (Fig. 2a) shows a strong peak at ð 1 2 1 2 1 2 Þ arising from second-order contaminations of the substrate peak, as well as a weakly split, rather broad film peak. Simulation of the neutron data (Fig. 2b, c) based on the XRD RSM ferroelastic variant densities (details are given in Supplemental Information), shows that the magnitude of the splitting corresponds rather to a cycloid period of 64 nm (as seen for the 192 nm film) than to a longer cycloid with 84 nm period (as observed for the 50 nm film). The 90 nm film behaves similarly (see Fig. 2d). In contrast, the singledomain 142 nm film (Fig. 2d) shows distinct peak splitting with magnitude similar to that observed for the 192 nm film. Nonetheless, the propagation direction remains the same for all samples, namely ½112.
The cycloid period for all films, derived from extracted line profiles along the ½112 direction, is shown in Fig. 2f (for more details see Supplementary Information). A distinct increase of the length of the spin cycloid in BFO with film thickness becomes apparent. For films thicker than~60 nm the period maintains a value comparable to bulk BFO as indicated by the horizontal line; while for thicknesses below 60 nm it asymptotically increases with decreasing thickness. Recall that the strain state of this series of films is not thickness dependent. This suggests that the increase in period must arise from factors related to finite-size and/or surface/ interface-driven effects. Gareeva et al. have theoretically shown that the period of the cycloid in BFO can increase upon modification of specific anisotropy terms. 39 In the present case, such a change in anisotropy could arise from local deviations in the DMI and other exchange interactions, which are sensitive to surface and interface effects. Note that the ½112 propagation direction of the spin cycloid corresponds to a direction which is 36°canted out of the film plane, and, as such, the endpoints of the spin chain are exposed to proximity effects (spin pinning) at the film surface and the film-substrate interface. This would be expected to have a more pronounced effect in films of thickness comparable to the cycloidal period. Real-space imaging techniques would give further information on the precise driving factors for the period lengthening of the spin cycloid. 38 Monte Carlo simulations To gain an understanding why the observed cycloid propagation vector is different from the ½110 direction reported for bulk BFO, we have investigated the magnetic structure of strained BFO theoretically, using Monte Carlo (MC) simulations based on a firstprinciples effective Hamiltonian (H eff ) method. This approach has been used to reproduce the structural and magnetic ground states of the R3c phase of bulk BFO, in particular the magnetic cycloid with propagation along the ½110 direction. 40,41 Here, we freeze in some components of the strain tensor according to the experimentally found lattice parameters (further detail is given in the Methods section).
The critical feature of our model is the incorporation of an energetic term analogous to the converse spin-current model, 42,43 in which the parameter C is a coefficient that captures the spincurrent (see the full Hamiltonian given in Methods). To characterize the magnetic structure found by MC simulations, we compute the discrete Fourier transform of the local magnetic moments in an 18 × 18 × 18 supercell. For a specific range of the spin-current parameter C, we find that the dominant k-points in reciprocal space (weighing more than 99% accounting of all possible k-points) correspond to a real-space propagation vector along the ½112 direction, as found experimentally in Figs. 1 and 2. The computed spin configuration for a single cycloid in the supercell is presented in Fig. 3.
The computations find (i) the angle between two neighbor magnetic moments along the propagation direction is 120°, which implies that such angle should be 2π 6 N radians in (larger and more realistic) N × N × N supercells; and (ii) all magnetic moments are found to lie mostly in the ð110Þ plane, with a weak but finite outof-plane component. More precisely, the out-of-plane component is found to be at least one order of magnitude smaller than the ð110Þ in-plane components, consistent with the experimentally observed magnetic structure.
Elongation of the spin cycloid with temperature In addition to finite-size effects, temperature also plays a dominant role on exchange interactions. Therefore we now address the influence of temperature on such interactions. Neutron diffraction data for the 142 nm thick 2% Co-doped BFO film acquired at different temperatures are presented in Fig. 4. Line profiles across the incommensurate split magnetic structure at ð 1 2 1 2 1 2 Þ are shown in Fig. 4a for temperatures from 175 to 575 K and indicate a decrease in peak separation upon increasing temperature. Fitting of two Gaussian functions to the data allows to quantify the reduction of the peak splitting and intensity with increased temperatures. A power-law fit to the temperature-dependent peak intensities (Fig. 4b) allowed us to determine the Néel temperature, T N . The value of 585(10) K is consistent with the hypothesis that Co-doping, like other B-site doping 11 reduces the magnetic transition temperature. In addition, theoretical calculations have demonstrated that the paraelectric-to-ferroelectric phase transition temperature decreases with film thickness in ultrathin ferroelectric films. 44 Figure 4c shows the extracted relationship between the cycloid periodicity and temperature, with a sharp increase in the spin cycloid period as the magnetic phase transition is approached. A similar behavior has been observed for the length of the spin helix in bulk SrFeO 3 . 45 The abrupt increase in cycloid period as T N is approached is significantly stronger than the change in single-crystal BFO (see solid line in Fig. 4c, from ref. 46 ) which may reflect a significant difference in the relative exchange interaction energies in thin films, where the inverse DMI responsible for the cycloid may be more sensitive than it is in bulk BFO. Therefore, thermal fluctuations may result in a pronounced weakening of the spin cycloid in thin films as compared to bulk BFO. Indeed, the critical temperature at which the cycloid period diverges appears slightly below T N , while in previous work on bulk BFO the cycloid appears to be in fact more stable as T N is approached. 46 Our findings confirm that finite size effects and/or surface and interface effects play a critical role in the stability of incommensurate spin structures. The systematic elongation of the cycloid period with thickness and as a function of temperature suggests The cycloidal length of bulk BFO of 64 nm is indicated as horizontal line. [14][15][16]18 The error bars correspond to 1 s.d. of Gaussian fits to the diffraction peaks. The inset shows the experimental scattering geometry of a (110)-oriented thin film with two ferroelastic domains with ferroelectric polarization vectors P 1 and P 2 S.R. Burns et al.
that the DMI in thin films is sensitive to both, size effects and thermal fluctuations. This demonstration of the modulation of the spin cycloid as a function of dimension and temperature allows us to gain a better understanding of the energy scales that govern spin order parameter and spin-lattice interactions in order to stabilize non-commensurate spin structures. These results thus offer important insight into the correlations between size and spin-lattice coupling in nanoscale multiferroic thin films, particularly in regimes where order-parameter interplay is predicted to yield unique multiferroic responses.

DISCUSSION
In summary, we have explored the systematic length variation of the spin cycloid in (110)-oriented Co-doped BiFeO 3 thin films. We find a cycloidal order with a propagation vector tilted out-of-plane in films as thin as 25 nm. In all investigated films, the propagation direction of the cycloid is found to be ½112 and effective Hamiltonian calculations based on first-principles allow us to explain why a uniaxially strained state gives rise to such a propagation direction. Remarkably, the length of the spin cycloid is rather robust with only a sharp increase for the thinnest film, probably due to finite size effects in thin films or spin-pinning effects at the surface and interface. Temperature-dependent neutron diffraction measurements indicate that the magnetic phase transition temperature is reduced by Co-doping or finite size effects and that the cycloidal length shows a sharp increase on approaching the antiferromagnetic phase transition. The knowledge of length variation of the spin cycloid in BiFeO 3 thin films will form important guidelines for design of future magnonic devices, since the magnon mode energies are dependent on the cycloid period. Our results offer therefore new perspectives for innovative magnonic and spintronic devices underpinned by controllable cycloid dimensions and properties.

Film fabrication and structural characterization
Epitaxial 2% Co-doped BFO films were grown on bare (110)-oriented SrTiO 3 substrates (size 8 × 8 mm 2 ) by PLD. The substrate temperature was maintained at 590°C, and the oxygen pressure at 100 mTorr. A ceramic Bi 1.1 Fe 0.98 Co 0.02 O 3 target was ablated using a KrF excimer laser (wavelength 248 nm; fluence~1 J/cm 2 ), at a repetition rate of 10 Hz. Such conditions provided growth rates of~0.04 Å/pulse. In contrast to our previous investigation, 10 no conductive SrRuO 3 intermediate layer was used. Standard X-ray diffraction characterization was performed using a Philips X-pert Pro MRD diffractometer, using K α−1 radiation. X-ray reciprocal space maps (RSM) were acquired using a 1D detector (Pixcel). It was found that the films were uniaxially strained along the [001] in-plane direction, while they were partially relaxed in the other directions. Full characterization of the crystallographic structure and surface quality of the obtained films is presented in Supplementary Figs. S1 and S2, respectively.

Neutron diffraction
Neutron diffraction measurements were performed using the triple-axis spectrometer TAIPAN, located at the Opal research reactor at ANSTO in Sydney, Australia. TAIPAN provides both the strong signal-to-noise ratio and low background necessary for measuring nanometer-thick films with neutron diffraction. 10,47 The instrument uses a vertically focussing highlyordered pyrolytic-graphite monochromator and analyzer which were both operated with a fixed wavelength of λ = 2.346 Å set for elastic neutron scattering. Two 40″ collimators were placed in the neutron beam before and after the sample to achieve the required resolution. In order to suppress contaminations arising from second-order diffraction emerging from the significantly thicker substrate, two pyrolytic graphite filters of 60 mm total thickness were placed in the neutron beam. The BFO films were oriented within the [110] × [001] scattering plane in order to access the ð 1

Effective Hamiltonian calculations
Following the work of refs. 40,41 , the total energy of the effective Hamiltonian for BFO is written as a summation of two main contributions: E tot ¼ E FEÀAFD ðfu i g; fη H g; fη i g; fω i gÞ þ E Mag ðfm i g; fu i g; fη i g; fω i gÞ where E FE-AFD contains energy terms arising from non-magnetic variables, 48 including the local mode (u i ) being proportional to the electric dipole of unit cell i, the homogeneous (η H ) and inhomogeneous strain in this unit cell (η i ), and the antiferrodistortive (AFD) mode (ω i ) that is associated with the tilting of oxygen octahedra in unit cell i. The second term E Mag is related to magnetism. It includes the mutual interaction between magnetic moments of Fe ions (m i ) with a fixed magnitude of 4 μ B (as consistent with first-principles and measurements 49,50 ) as well as couplings between magnetic moments and the other degrees of freedom, viz., local modes, AFD motions, and strains. Technically, the analytical form of E Mag can be expressed as: 40,41 E Mag ðfm i g; fu i g; fη i g; The first term in the above expression represents the magnetic dipolar interaction where i and j run over all the sites. The second term is the magnetic exchange coupling between magnetic moments up to thirdnearest neighbor (NN). The third, fourth, and fifth terms characterize the change in magnetic exchange interaction induced by the local modes, AFD motions, and strains. The index j for the second, third, fourth, and fifth term runs over first, second, and third NN of site i. Note that the first five energies can only lead to collinear magnetism. On the other hand, the sixth term, which involves the rotations of the oxygen octahedral and in which the j index runs over the six first NN of site i, is responsible for the spin-canted weak magnetization of BFO. [51][52][53][54] The last term characterizes the converse spin-current model, 42,43 with e ij being the unit vector along the direction joining site i to site j. Interestingly, these sites i and j were assumed to be second NN of each other in refs. 40,41,55 , which allowed to reproduce the ½110 propagation direction of the cycloidal magnetic structure of bulk BFO. 40,41,55 However, we numerically and presently find that first NNs are dominant over the second NNs to obtain the experimentally found cycloids. As a result and consistent with a recent work, 56 we now choose the j sites being (the six) first NN of site i, which implies that e ij is a unit vector lying along the different 〈100〉 pseudo-cubic directions. Such choice is found to lead successfully to the generation of cycloids propagating along ½110 in bulk BFO as well as along ½112 in (110) BFO films. Note also that the prefactor C is a material-dependent coefficient, and characterizes the coupling strength of the spin-current interaction. In our simulations, C is allowed to vary to accommodate the size of the supercell and periodic boundary conditions. Other parameters in the effective Hamiltonian are obtained from first-principles densityfunctional calculations. We carry out MC simulations using the total energy E tot with 12 × 12 × 12 (8640 atoms) and 18 × 18 × 18 (29,160 atoms) supercells, in terms of the 5atom pseudo-cubic perovskite unit cells. The geometry of the lattice vectors of the supercell is constrained according to the strain found in the experimentally used (110) epitaxial BFO films, i.e., +1.15% along the pseudo-cubic [110] direction, +0.25% along the ½110 direction, and −1.40% along the [001] direction, while the atomic internal coordinates are allowed to vary. Such a strain condition generates in our simulations the so-called M B monoclinic phase 57 for which the component of the electrical polarization along [001] is slightly smaller than the (identical) components of such polarization along [100] and [010], being in agreement with measurements. To find the stable magnetic structure, we first equilibrate the system up to 200,000 MC sweeps at 600 K, starting with the cycloidal structure of bulk BFO in the R3c phase, followed by MC simulations at 30 K for up to 800,000 sweeps. Finally, the system is further relaxed at 1 K for 200,000 sweeps for the analysis of the magnetic structure. Note that this low temperature is chosen for better statistics.
Moreover, the in-plane components of the magnetic structure in our 12 × 12 × 12 and 18 × 18 × 18 supercells were found to be rather well described by the following equation: where m i is the ð110Þ in-plane component of the magnetic moment vector at site i. m R and m I are two vectors that can be easily determined by two calculated magnetic moments at two non-equivalent real-space positions r i . k m is the k-point with the largest Fourier transform. A spin cycloid with the measured period can be constructed using this formula, as illustrated in Fig. 3, thus implying that such a cycloid is harmonic in nature.

DATA AVAILABILITY
Data supporting the findings in this study are available within the article or from the corresponding authors upon request.  2 Þ antiferromagnetic Bragg peak at selected temperatures. The lines are the results of Gaussian peak fits to the data. b Temperature dependence of the peak intensities. The analysis with a power law behavior indicates a magnetic phase transition temperature of T N = 585(10) K. c The length of the spin cycloid as extracted from the splitting of the two magnetic Bragg peaks, expands close to the magnetic phase transition. In addition, the change in the cycloidal length for bulk BFO is shown as gray line (extracted from Ramazanoglu et al. 46 ). The error bars in (a) correspond to 1 s.d. of the count rate while in (b, c) they correspond to 1 s.d. of Gaussian fits to the data