Topologically stable helices in exchange coupled rare-earth/rare-earth multilayer with superspin-glass like ordering

Existence of 2 π -planar domain walls (DWs) are often reported for transition metal-rare-earth (TM/RE) layered systems. The magnetization process of such two-dimensional randomly anisotropical system in the form of 2 π -DWs is directly correlated with topologically stable helices. Here, instead of TM/RE, we have investigated [Dy/Tb] 10 multilayers involving two different anisotropic layers of rare-earth/rare-earth (RE/RE). Using magnetization and susceptibility as function of temperature along with thermo-remanent magnetization measurements we have con ﬁ rmed superspin-glass type of behavior within this RE/RE system. Additionally, an exchange bias ﬁ eld up to – 0.88 kOe ( – 88 mT) was also revealed for such rare-earth based multilayers. Interestingly, using detailed analysis of the polarized neutron re ﬂ ectometry pro ﬁ les, we ﬁ nd evidences of superimposed helical magnetic con ﬁ gurations within both materials of Dy and Tb associated with spin-frustrated interfaces. Furthermore, magnetization ﬂ uctuations around the mean magnetization from vertically uncorrelated domains were observed with polarized off-specular neutron scattering. We believe that coexistence of helical ground states with superspin-glass-like ordering are fundamentally instrumental for topologically stability in RE/RE systems, which in principle, can be exploited in all-spin-based technology.

T hin films composed of rare-earth (RE) metals are known to form a helical antiferromagnetic (AF) structure due to competing magnetoelastic, crystalline anisotropy, and exchange interactions [1][2][3] . These interactions can lead to topologically stable spin configurations, which are virtually indispensable in all-spin-based technology 3 . Numerically, Vedmedenko et al. have shown that magnetic states corresponding to modulated helices with integer number of twists, commensurate with the chain's length, are topologically stable 4 . Such stability comes without the presence of chiral Dzyaloshinskii-Moriya interaction relevant in skyrmions, where its shape protects it from trivial unwinding. The quantized energy spectrum within finite magnetic chains results in stabilizing topologically protected configurations of helices. Such helices, initially twisted in an external magnetic field, stay stable even without the presence of a field as they are created by exchange interaction and does not allow a trivial unwinding 5 .
Using polarized neutron reflectivity (PNR) measurements, it was shown earlier that AF-coupling at a pure transition metalrare earth (TM/RE) (e.g., Fe/Tb) interface felicitates the formation of π or 2π planar domain walls (DWs) 6,7 . Within a multilayer, formation of such 2π-DWs resemble a helical form leading to double hysteresis loop (DHL). DHLs are associated with exchange-bias (EB)-like shifts along and opposite to the field cooling axis below the ordering temperature of RE. The energy it takes to form a DW within the soft layer, which remains frozen upon cooling, results in the EB shift. Subsequently, exchange coupling and helical phase characteristics were also reported within both layers of Fe and Dy, comprising the multilayer 8 . At low fields, the TM/RE layers exhibit ferrimagnetic spin alignment, which form twisted magnetic helices in the form of planar 2π-DWs. However, at higher fields, a more complicated but stable and continuous helical arrangement was observed. In another case, evidence of such stable helices within single layer Er films was shown using PNR 2 .
Even though TM/RE multilayers are commonly studied, RE/RE multilayers remain scantily probed. Some combinations of epitaxial RE/Y multilayers showing modified helical moment configurations had been reported in the past using neutron diffraction 9 . Epitaxial erbium thin films and Er-Y superlattices were reported to show the effects of basal-plane strain on the modulated spin structure by Borchers et al. 10 . However, all strain related phenomena are expected in epitaxial films and are insignificant in polycrystalline ones. Moreover, Y is a transition metal but is often classed as RE due to its similar electronic behavior. Note that exchange coupling has also been observed in several TM/RE systems [11][12][13] , but not in RE/RE systems. In contrast to the direct exchange in the transition elements, RE magnetism finds its explanation in an indirect exchange interaction known as the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction. Experimentally, it is found that in addition to normal ferromagnetism, some of the rare earths exhibit a helical AF state over a limited temperature range. On application of an external magnetic field, the helical state is radically distorted and eventually collapses into a ferromagnetic configuration. Presence of large orbital momentum in Tb, Dy, and Ho leads to strong spin-orbit coupling and larger magnetic anisotropy 14 . Moreover, a large difference in spin-orbit coupling for different RE elements, would consequently have a significant influence on the overall demagnetization processes. In a recent example of RE/RE system, Er-Tb superlattice has been shown to posses long range magnetic order despite their competing anisotropy directions. Probing the vertically correlated magnetic structures by off-specular polarized neutron scattering the existence of magnetic vortex-like domains associated with magnetic helical ordering within the Er layers was evidenced 3 .
Here, we investigate polycrystalline [Dy10/Tb10] 10 : [Dy(10.0 nm)/Tb(10.0 nm)] ×N=10 multilayers on Si substrates where N is the number of bilayer repetitions. The temperature ranges where Tb and Dy show AF helical magnetization, are distinctly different in bulk. For Tb, this range hovers around 220 K, whereas for Dy it ranges approximately between 180-80 K, below which both are ferromagnetic 6,8 . We found that the multilayer is magnetically very hard, with a saturation magnetization beyond 50 kOe (5000 mT). Most interestingly, we observed an EB coupling up to -0.88 kOe (-88 mT) due to strong pinning at spin-frustrated interfaces. Direct current (DC) magnetization as a function of temperature and alternating current (AC) susceptibility measurements for various frequencies and temperature confirmed such superspinglass (SSG) type of behavior for the multilayer. Additionally, thermo-remanent magnetization (TRM) confirms the dynamics in the SSG phase and test its scaling behavior. Using PNR and polarized off-specular neutron scattering, we determined the magnetic profile of the multilayer indicating spin configurations like helices and associated magnetic fluctuations within the stack. Double spin spirals within Dy and Tb leads to spin-frustrated SSG-like ordering at the interfaces. Stability of the 2π-DWs at a higher field ensures their topological stability. Our study is expected to open up a path for a broader base of researches as it describes novel phases of long range helical ordering in RE/RE polycrystalline system across spin-frustrated interfaces with innumerable possible combinations to explore in the future.

Results
X-ray characterization. For details on the initial characterization such as X-ray reflectivity (XRR) and X-ray diffraction (XRD), see the discussion in Supplementary Note 1 and Supplementary  Fig. 1.
Magnetization measurements. In order to characterize the magnetic properties of the multilayer, the magnetization M(T) was measured as a function of temperature for [Dy10/Tb10] 10 . We applied different magnetic field strengths H a = 0.25, 0.5, 1.0, 2.5, 5.0, 10.0, and 50.0 kOe (25, 50, 100, 250, 500, 1000, and 5000 mT) during the measurements. All samples were either initially cooled down to 5 K in presence of a field H FC = 70 kOe (7000 mT) to obtain the field cooled (FC) curves or with no magnetic field to obtain the zero field cooled (ZFC) curves. The M(T) curves are shown in Fig. 1a, b. The ZFC curves show a broad peak which can be referred to as blocking/freezing temperature (T F ) whereas the FC curves do not show a temperature independent platform-like behavior. A well-defined irreversibility temperature (T irr ), which is indicated by the temperature where FC and ZFC curves diverge, could be identified. The gradual convergence of ZFC and FC curves with increasing field, indicates the attainment of similar type of magnetic configuration near equilibrium as exemplified for 50 kOe (5000 mT). As both 5/5 ([Dy5/Tb5] 10 10 sample only. Note that for neutron measurements, sub-loops were measured for a limited FC condition of 10.0 kOe (1000 mT) as shown in Supplementary Fig. 4.
The transition from ferromagnetism to superparamagnetism (SPM) or super spin glass (SSG) behavior is generally expected for discrete nano-clusters (collection of nanoparticles) where individual magnetic moments within clusters are thermally unstable. The SSG state is believed to result from frustration generated by dipole-dipole interactions among superspins (magnetic moments of nano-clusters) and from disorders in the system (e.g., the random distributions of clusters, positions, sizes and anisotropyaxis orientations) 15 . The irreversibility temperature (T irr ) and the blocking or freezing temperature T F as a function of field is shown in Fig. 1c for the [Dy10/Tb10] 10 multilayer. T F is slightly below T irr , which in some respects resembles cluster-glassy behavior 16 . They both show a gradual decrease which suggests that the frozen state is suppressed by the field (a similar behavior is seen for [Dy5/Tb5] 10 : see Supplementary Fig. 2).
The shift of T irr to lower temperatures can follow the Almeida-Thouless (AT) line indicating a SSG-like behavior 17,18 . The AT line is expressed as where T irr (0) is the zero field freezing temperature and ΔJ is the width of the distribution of exchange interactions. The plot of the measuring field H 2=3 a as a function of T irr has been shown in Fig. 1d, which fits to a straight line following the AT equation. The line cuts the x axis at around T irr (0) = 205 K for [Dy10/ Tb10] 10 (T irr (0) = 210 K for [Dy5/Tb5] 10 : see Supplementary  Fig. 2) which is similar to T F . AT line separates a nonergodic (SSG) phase from a ergodic (SPM) one. In the inset of Fig. 1d we additionally plot H a 2/3 as a function of 1 − T irr (H a )/T irr (0), showing the transition line separating the SPM and SSG phases 19 .
In-plane magnetic field hysteresis loops were measured at seven different temperatures after field cooling in 70 kOe (7000 mT) from 300 K and are shown in Fig. 2a Fig. 2b for the [Dy10/Tb10] 10 multilayer. One can see a gradual decrease of the coercivity H c with increasing temperature. The H EB value decreases gradually with increase in temperature and vanishes at around 200 K. Note that the temperature dependence of the remanent magnetizations M R = (|M| (H) + |M| (−H) )/2 (red circles in Fig. 2b) goes to zero also at 200 K (a similar behavior was seen for the [Dy5/ Tb5] 10 multilayer as well).
Dy exhibits a rich magnetic phase diagram, including a few modulated magnetic phases and also posses one of the highest intrinsic magnetic moments (bulk value of 10.6 μ B /atom) 20,21 . The magnetic modulations propagate coherently over a long range, even with intervening nonmagnetic layers, via RKKY interaction. The magnetic moment of Tb (with bulk value of 9.0 μ B /atom and 6.4 μ B /atom at 10 K in epitaxial thin film 22 ) is slightly lower than Dy. The range of stability of the helical state in Tb is quite restricted compared to Dy. An applied field of less than 1 kOe (100 mT) will produce a transition to the ferromagnetic state in Tb, while in Dy the reminiscent of the helical state can be found even at lower temperatures 8 . Moreover, due to basal-plane anisotropies existing both in Dy (a axis) and Tb (b axis), a combination of these two RE elements may lead to an exchange coupling as well. With an increase in spin-frustration in the pinning layers at the interfaces, a transition from a helical AF domain state to a spin-glass type state might be expected which can eventually increase the H EB field.
If we assume that randomly oriented nano-clusters are following the usually observed T    10 ). One may note that the peaks in the respective ZFC curves, corresponding to the average SPM/SSG freezing temperature T F , is similar to T F (0). However, the same plot can also be subjected to an exponential fit which leads to an ambiguity. In such case of low interacting system, it often becomes difficult to distinguish between the SSG and SPM type of phases from the analysis of the hysteresis curves alone 19 . The rate of change of d 2 H EB /dT 2 with respect to temperature does not exhibit any peak as was reported earlier for SSG nanoparticles 24 . Thus, the appearance of H EB and the onset of spin blocking/freezing (T F (0)) may not be interdependent. The physical origin of exchange-coupling is rather due to the two hard magnetic layers with two different anisotropies 5 .
We measured the temperature dependence of the AC magnetic susceptibility over a frequency range of 10 to 10,000 Hz in the presence of a small AC field of 0.01 kOe (1 mT) to differentiate between the SPM an SSG characteristic. Figure 3a shows the frequency dependence of the real part of the susceptibility (χ′(T)) versus temperature for [Dy10/ Tb10] 10 . The curves show peaks centered at around T AC max ' 206 K. The paramagnetic Curie temperature at around T C = 230 K can also be estimated from a 1/χ′(T) versus temperature plot. The T AC max peaks decrease in magnitude and shift to higher temperature with increasing driving frequency. The shift gives the activation energy and is characteristic of the SPM-type or SSG-type behavior. Non-interacting SPM clusters are expected to show larger frequency dependence than SSG clusters, because the distribution of relaxation times is characteristic for the SSG phase. To qualitatively analyze the dynamical behavior, three common models are employed in discerning the supermagnetic dynamics. The models are embodied in the (i) Néel-Arrhenius law 25,26 , (ii) Vogel-Fulcher law 27,28 and (iii) power law 29 in addition to the (iv) empirical equation 7,30 .
where τ 0 = (2πf 0 ) −1 is the relaxation time for the attempt frequency f 0 and τ m = (2πf m ) −1 is the measuring time (~10 2 s for the DC magnetization measurement) for the measuring frequency f m . Here, ΔE (=K A V) is the activation energy, K A is the anisotropy constant, k B (=1.38 × 10 −16 erg/K) is the Boltzmann constant and V is the average cluster volume critical for SPM or SSG state. The value of τ 0 typically ranges from~10 −9 -10 −13 s for a SPM behavior 31 . In Fig. 3b, the dependence of T AC max À Á À1 on the natural logarithm of the measurement frequency is plotted. A linear fit to the data yields a slope of ΔE/k B = 190,344 K, which implies an unreasonably small value of τ 0 . Failure of the Néel-Arrhenius expression in describing the dynamics suggests the presence of non-ignorable magnetic dipole-dipole interactions. Weak interactions among spin nano-clusters are taken into account within the Vogel-Fulcher law 27,28 . Phenomenologically, it describes the frequency response of the relaxation time and is expressed by where T 0 is the characteristic temperature that accounts for the static interaction field of the surrounding clusters. In Fig. 3c, a plot of ln(τ m ) versus T AC max is shown. The values obtained from a fit to the equation 3 are ΔE/k B~1 12 K, T 0~2 03 K and τ 0~1 0 −13 s. From the linear fits to the plots of T AC max versus 1/ln(τ 0 /τ m ), there can be a few possible combinations where τ 0 can be within a reasonable range yielding different ΔE/k B values 32 . For example, τ 0 = 10 −17 s corresponds to ΔE/k B = 214 ± 12 K and T 0 = 201 ± 0.4 K while τ 0 = 10 −9 s corresponds to ΔE/k B = 35 ± 2 K and T 0 = 205 ± 0.2 K (see the discussion under AC susceptibilities within Supplementary Note 2 and subsequent Supplementary Fig. 5a).  These values are reasonably comparable with the conventional spin glass behavior for which the typical value of τ 0~1 0 −10 -10 −13 s or even lower 33 . A longer τ 0 yields dynamics on a shorter accessible time scale and are probed within the experimental time window available for a SSG type of behavior. The maxima observed in χ′(T), would therefore correspond to the freezing of weakly interacting moments of the nano-clusters.
The power law relies on the relation of critical slowing down of the relaxation time near the transition temperature. It is also known as the scaling hypothesis, which predicts the existence of an equilibrium phase transition 29 . The relaxation behavior corresponding to a characteristic temperature T SSG for superspin-glass behavior is expressed by where τ * is a relaxation time for each nano-cluster, zν is the dynamical scaling critical exponent constant 34 related to the correlation length ξ, which is a measure of the size of the lateral coarsening (ν describes the divergence while z is involved in the dynamical scaling hypothesis τ m~ξ z ). Tholence criterion requires a very small value for η ¼ ðT AC max À T SSG Þ=T SSG in order to comply with the SSG type of behavior. In Fig. 3d we plot the variation of τ m (in log scale) versus T AC max . A fit to the Eq. (4) yields zν~5.0 (typical range for SSG is between 4 to 12 ) , τ *~3 × 10 −15 s, and T SSG~2 06.8 K. Considering this T SSG value, the maximum value for η~0.01. Alternatively, linear fits to the plots of ln(τ m ) versus ln½ðT AC max À T SSG Þ=T SSG yield a few possible values, for e.g., zν = 4.0 ± 0 . 4 and τ * = 3 × 10 −13±1 s for T SSG = 207.0 K and zν = 6.0 ± 0 . 4 and τ * = 5 × 10 −17±1 s for T SSG = 206.6 K (see Supplementary  Fig. 5b). These values are fairly comparable with the typical values reported for a SSG type of behavior 35 . One may note that the values of T F (0), T irr (0) and T SSG all coincide roughly with the helical ordering temperature in bulk Tb.
Lastly, another useful and sensitive criterion to distinguish between the freezing and the blocking processes is to determine the relative shift of the χ′(T) peak with frequency using the empirical equation where T AC max is the mean value of the frequency dependent maximum in χ′(T), whereas ΔT AC max is the difference in T AC max over the frequency interval Δlog 10 (f m ) 35,36 . Typically, for SSG phase, the parameter p assumes values of 0.0045-0.06 whereas for SPM phase, it is between 0.1-0.5 30 . The parameter p for the SPM phase is typically an order or two higher than that for the SSG phase. Here, we obtain p = 0.0028, which again points towards an interacting SSG-type of behavior 7 . The obtained lower value of p indicates that the accessible time scale for SSG dynamics is indeed shorter.
For TRM measurements 37 , the sample is heated to a temperature above the superspin-glass transition temperature (T SSG = 205 K) and subsequently cooled to a measuring temperature (T m = 143 K) in the presence of a small excitation field, H TRM = 0.005 kOe (0.5 mT). After a waiting time t w , the field is switched off and the relaxation of the magnetization is measured over t = 30,000 s. A waiting time dependence was examined to evaluate the aging character of the dynamics. The field should remain small enough so that the response of the sample remains in the linear regime (non-perturbative probe) to avoid an influence of the field on aging. Figure 4a shows the relaxation of magnetization M TRM (t) for [Dy10/Tb10] 10 , normalized to the field cooled magnetization (M FC ) for values of waiting times t w varying from 30 s, 300 s, 1,000 s, 3,000 s to 30,000 s. The plot for t w = 0 s (t 0 ) is used for reference. It can be seen that the relaxation depends on the value of t w : this illustrates the aging character of the dynamics. Longer the t w , slower is the relaxation, indicating a stiffening of the sample response during the waiting time.  is slightly different from the commonly defined derivative term. Nevertheless, the quantity S is equivalent to the relaxation time distribution, a wide breadth of which is the reason behind the slow and non-exponential relaxation of the response function in a SSG state. The inflection point t infl in the curves can be found from the S(t) versus log(t) plot in Fig. 4b. Two features distinguish a SSG behavior from an atomic spin glass behavior. Firstly, a departure from approximate scaling in the TRM curves versus the scaling variable λ/t w plot in Fig. 4c where the effective time λ = t w [(1 + t/t w ) 1−μ − 1]/[1 − μ] with μ = 0.5, which is much smaller than is normally found for atomic spin glasses (0.7-0.9) and is similar to that of SSG (0.4) 38 . Values of μ ≠ 1 indicate by how much the effective age of a SSG deviates from the nominal age; i.e., the experimental waiting time, t w 15 . Secondly, log(t infl ) ≠ log(t w ), which is unlikely for an atomic spin glass where an approximate scaling is usually observed (Fig. 4d). Thus, the SSG behavior is well established for the multilayer under investigation.
Polarized neutron scattering measurements. For the investigation of the magnetic structure of the sample, PNR measurements were carried out. While specular scattering is used for depth information, off-specular scattering can deliver information on lateral structure. A sketch of the neutron scattering geometry is shown in Fig. 5a. Figure 5b shows the sketch of lateral magnetic domains with respect to the lateral projection of the neutron coherence length l || , which is deterministic in observing specular and off-specular scattering. For details see the "Methods" section. PNR measurements were performed for different magnetic histories of [Dy10/Tb10] 10 . Figure shows the specular polarized neutron intensity profiles along Q z (see the "Methods" section for details) and their fits. The measurements were done after cooling the sample to 10 K in zero field (ZFC) (Fig. 6a) and in +10.0 kOe or +1000 mT (FC) (Fig. 6b). Both measurements took place at an external field H a = -0.01 kOe (-1 mT). The corresponding hysteresis sub-loops measured at 10 K are shown for the FC condition in Supplementary Fig. 4. The fits were initialized using simple models of block-potentials. The parameters used for fitting were the individual layer thicknesses along with the nuclear and magnetic scattering length densities (SLDs) of individual layers. The errors in the thickness of the layers are ±0.2 nm, while that for the nuclear and magnetic SLDs ρ n and ρ m are ±0.1 × 10 −6 Å −2 and ±0.05 × 10 −6 Å −2 , respectively. As the top layer is susceptible to oxidation, we assumed a different ρ n and ρ m values when compared with the stack. The interface roughness is 2:0ð± 0:5Þ nm. Fitted parameters were obtained using the minimization of the reduced chi-squared (χ 2 or the goodness of fit) value. The insets in the figure for PNR show the fitted SLDs values ρ n and ρ m versus the thickness of the multilayer.
Due to the limitation in providing the maximum applied field at the experimental station, the saturation moments could not be determined. The ZFC data (Fig. 6a) at H a = -0.01 kOe (-1 mT) (virgin state) was found similar to the room temperature data with no net magnetic moment. Absence of any net magnetization is due to the formation of domain states in accordance with the ZFC data in Fig. 1a. The peaks at Q z = 0.038 Å −1 and Q z = 0.109 Å −1 are the 1st order and 3rd Bragg peaks due to the structural periodicity. The 2nd order structural Bragg peak is expected to be significantly reduced since the layer thicknesses are almost equal. For the structural SLDs, the fits yielded the following values: ρ Tb n ¼ 2:5 10 À6 Å −2 , ρ Dy n ¼ 4:8 10 À6 Å −2 . We also found an interdiffused layer of around 0.5 nm and 1.0 nm at the Dy/Tb and Tb/Dy interfaces, respectively. This diffusion layer was also evident in the XRR measurements (estimated from the fits to the XRR data in Supplementary Fig. 1). In the case of a random distribution of domain magnetization directions, the dispersion is rotation this dispersion is essentially zero. One should not expect here any SF specular signal due to the formation of domain states and no net magnetism, therefore the SF specular intensity visible is attributed purely to the inefficiencies (2-5%) of the optical elements. Note that the so-called leaked SF specular signal is at least two orders of magnitude smaller than the corresponding NSF specular signal. The plot in the inset show the depth profiles in terms of ρ n values for the two layers while the ρ m value is zero for both.
The nuclear SLDs, thus obtained, were used to fit the FC data ( Fig. 6b) at H a = -0.01 kOe (-1 mT) and obtain the magnetic SLD parameters close to the reduced-remanence field of the sub-loop hysteresis. Here again, the fit was started with a simple model of block-like-potentials with average magnetic SLD values to fit the profiles. At this field, the magnetic SLD values obtained from the fits are: ρ Tb m ¼ À1:4 10 À6 Å −2 , ρ Dy m ¼ À1:6 10 À6 Å −2 . The corresponding magnetic moments are Tb: 1.7 μ B /atom and Dy: 1.9 μ B /atom. These values are found to be much lower than those known for bulk or values reported earlier (polycrystalline Tb: 2.4 μ B /atom 6 and polycrystalline Dy: 7.1 μ B /atom 8 ) for thin films. This lowering of moment is caused by the fact that the FC process was done not in a saturation field but in an intermediate field.
The top layer was found to posses a similar average moment and was associated with a small rotation, oriented 10(±5) away from the applied field direction. This disorientation gave rise to a small SF signal. Such an increase in the SF signal also was reported earlier for another multilayer system 39 . The difference between the two SF curves is basically due to a difference in the instrumental (polarizer and analyzer) inefficiencies. The plot in the inset show the ρ n and ρ m values for the two layers as function of thickness.
Similarly, fits to the FC profiles (Fig. 7) at H a = -10.0 kOe (-1000 mT) were used to obtain the magnetic SLD parameters near the coercive field. Like that for the ZFC measurement at H a = -0.01 kOe (-1 mT), no true SF signal could be observed at this field, which indicates existence of either domain states or compensated moments 40 Figs. 6 and 7). Out of these, two of the most probable scenarios are discussed below.
In Furthermore, the magnetization vector of each layer to rotate in-plane (x − y) by an angle ϕ A(z) within the multilayer stack is considered, which mimics the existence of 2π-planar DWs within both layers of Dy and Tb 6 . The scenario is shown in Fig. 7b. The plot in the inset show the ρ n and ρ m values obtained for the two layers as a function of depth. The larger Dy and Tb moments are confined within a range of a few nanometers from the core of the layers. The planar DWs, which are compressed with an increase in field, are spanned on both sides of the remaining layer while they penetrate on either side of the adjacent layer. Extension of continuous helices, from one layer to the other, is owed to ferromagnetic coupling of the core layers and formation of DWs due to the pinning at interfaces. Within each layer, integration of magnetic moments show little or no variation around the core, whereas around the interfaces the DWs being laid spanned results in a lowering of magnetization. The improved quality of the fit with planar DWs indicates a possibility of our helical modulation model. Since no additional parameters were included other than those mentioned we believe that this is the most logical consequence. Earlier, helical magnetic structures or 2π-DWs have been reported for similar systems 2,3,8 .
Unlike TM/RE systems, here the layers are neither ferrimagnetically or antiferromagnetically aligned but are ferromagnetically aligned to each other. Signatures of magnetic helices in Dy has have been reported earlier in presence of low fields (0.03-0.1 kOe or 3-10 mT) even at 5 K 20,41 . The magnetic ordering of helical modulation is sensitive to the cooling history and interface clamping. Thus, one can expect a similar helical phase within Dy and Tb, which are merely reminiscent of the helical ordering realized at higher temperatures. The magnetic stability of the confined helices, demonstrated at -10.0 kOe (-1000 mT), can be attributed not only to the internal field that is created by exchange interaction which stores the magnetic energy but also to the strong anisotropy of the Dy layer (axial: 10 9 -10 8 erg/cm 3 ; basalplane: 10 7 -10 6 erg/cm 3 ). Note the difference in the basal-plane anisotropy constant for Tb (2.4 × 10 6 erg/cm 3 ) and Dy (7.5 × 10 6 erg/cm 3 ) at 4 K 42 . This difference in anisotropy leads to planar DW formation whose energy is modified by exchange interactions. The anisotropy in both Tb and Dy keeps the magnetic moments confined to the basal plane.
We show the specular and off-specular NSF and SF intensity maps in Fig. 8a, b for the multilayer with two different field histories, ZFC and FC and measured at H a = -0.01 kOe (-1 mT) and =-10.0 kOe (-1000 mT), respectively. The scenario for FC condition and measured at -0.01 kOe (-1 mT) also showed a similar intensity map (not shown). The measurement channels shown are the NSF channels I ++ , I −− and the SF channel I −+ . The intensity along the diagonal α i = α f is the specular reflection along the scattering vector Q ⊥ and has been discussed in the last section 43 . Apart from the multilayer Bragg peaks along the specular ridge depicting the periodicities of the multilayer, we could observe off-specular intensities in form of sheets of intensities perpendicular to the specular line (Bragg sheets) emanating across the Bragg peaks (I ++ and I −− channels) and intensities near the critical edge and along the Yoneda wings. Additional intensities along the Bragg peaks, parallel to α i,f (diffuse streaks in the I −+ channel) indicate configuration of a domain state for sample magnetization perpendicular to the H a direction.
Similar intensities of the specular Bragg peaks in the NSF channels indicate that they were measured either at -0.01 kOe (-1 mT) after ZFC (Fig. 8a) or at the coercive field of -10.0 kOe (-1000 mT) after FC (Fig. 8b). A slight imbalance in intensity along the specular ridge at = -10.0 kOe (-1000 mT) between I ++ and I −− channels indicate a preferential alignment of the magnetization, which is also in accordance with the fits to the specular profiles shown in Fig. 7b.
In the SF channel, there exists no off-specular Bragg sheet. However, off-specular diffuse streaks exist in the NSF channels as well as in the SF channels, indicating a domain state configuration of the sample magnetization also along and perpendicular to the H a direction (parallel and perpendicular components). The Bragg sheets in the NSF channels remain similarly intense with a change in the H a value while the diffuse streaks (indicated by arrows) are more intense in the I ++ channel for ZFC and in the I −− channel for FC. Similarly, in the SF channel I −+ , it is more intense along α i at -0.01 kOe (-1 mT) as compared to that along α f at = -10.0 kOe (-1000 mT), which obviously indicates that the diffuse streaks are of magnetic origin while the Bragg sheets, which remain similar at both fields, are not of magnetic origin. Such a temperature or field dependence of off-specular SF scattering from vertically correlated domains was reported earlier by Paul et al. in epitaxial Er/Tb system 3  with 2π-domain walls in both layers of Dy and Tb. Shown alongside are the nuclear (ρ n ) and magnetic (ρ m ) scattering length densities (SLDs) versus thickness of the multilayer certain field without any change in the field direction, they are supposed to remain asymmetric with respect to the diagonal intensities, i.e., I +− ≠ I −+ . Thus, the change in asymmetry of the scattered intensity in the respective SF channel I −+ for ZFC and FC in Fig. 8a, b, accounts for the fact that the spin states of the incoming and outgoing neutron waves are inversely populated indicating that they have been reversed by the field. This interchange of intensity obviously confirms their magnetic origin.
On the one hand, the Bragg sheets occur due to pronounced structural and/or magnetic vertical correlation along the stack of the multilayer. On the other hand, diffuse streaks run parallel to α i,f axes. Due to some instability in the system, one often witnesses fluctuations of the magnetization (Δϕ A ) from vertically uncorrelated domains around the mean magnetization in the form of diffuse scattering. Perpendicular fluctuations are characterized by the mean value sin 2 ðΔϕ A Þ, which are commonly seen in the SF channels below a saturation field 40,43 .
Distorted wave Born approximation (DWBA) has been applied in simulating the corresponding SF intensity maps in the SF I −+ channel (Fig. 8a, b). The parameters for the offspecular data simulations follow from the fit parameters obtained from the fitting of the corresponding specular data. Thus the inefficiencies of the polarizer and analyzer are intrinsically incorporated in the simulations. One can see an increase in the off-specular SF intensities near the critical edge in the simulated maps. These are typical signatures of random non-collinear arrangement of small-scale domains with sizes ranging from 1.5 ± 0.5 μm and 0.5 ± 0.5 μm accompanied by fluctuations of magnetization Δϕ A = 4 and 15, for ZFC and FC conditions measured at -0.01 kOe (-1 mT) and -10.0 kOe (-1000 mT), respectively. The domain sizes are typically smaller than the neutron coherence length along the x axis (few microns), i.e., as and when the coherence ellipse covers several domains 40,43 . The change in asymmetry for ZFC and FC has duely been accounted for in the respective simulations.
Noncollinear configuration of planar spin chain. Theoretically, Popov et al. 44 treated in-plane exchange interaction within mean field approximation, the energy of which can be described by a finite size 1D model for N planar classical spins Here, h A = H A /H EB and h = H a /H EB , where H EB , H A , and H a are the effective exchange field, effective uniaxial anisotropy field and external applied field, respectively. A spin on i th side of the chain forms an angle ϕ i A perpendicular to the chain axis. A number of collinear ground state and noncollinear equilibrium stable configurations were numerically calculated. Even though the theoretical discussion is meant for epitaxial films 9 , we could find it relevant for polycrystalline films as well emphasizing the robustness of spin configurations and thereby the generality of results.
A multilayer stack can relax into different equilibrium states depending upon the energy density introduced by an external magnetic field. With removal of the external field, the multilayer typically relaxes into a ground state as all magnetic moments remain aligned and attains a global minimum. However, for an intermediate field, the stack relaxes into a metastable state as the spins remain misoriented with the attainment of a number of local energy minima 5,8 . In the present context, Fig. 6b represents a scrambled ground state at -0.01 kOe (-1 mT) where all spins are randomly oriented while Fig. 7b represents a state of noncollinear equilibrium configurations at -10.0 kOe (-1000 mT) where spins are appreciably misoriented with respect to the hard axis as expected for an intermediate field. Supplementary Fig. 8a shows a schematic of the noncollinear spin configurations containing two π-DWs (topologically protected non-trivial structure) so that the overall rotation of the spins along the chain is 2π within the layers N ξ > 2π of Dy and Tb. Alternatively, a combination of π-DW and π-antiDW (topologically unstable trivial structure) of width ξ is possible. It may be noted that PNR alone cannot distinguish between the two configurational states. However, magnetic stability of the 2π-DWs or confined helices at -10.0 kOe (-1000 mT), exemplifies the fact that a non-trivial topological structure cannot be easily destroyed by a field. In order to reduce the cost of uniaxial anisotropy the spin orientations remain close to an integer multiple of π at the end points of the chain. Due to halved number of nearest neighbors, the spins at the surface/interface (centers of DW nucleation) were shown to be two times less, which reduces the energy ðE w ¼ 4a ffiffiffiffiffiffiffiffi ffi J ex K p Þ near the end points by a factor of ffiffi ffi 2 p . The DW width, which is given by ξ ¼ a 2 ffiffi γ p ' π ffiffi ffi A K q , leads to an uneven distribution of spins within an individual layer with relative compression at the interfaces and decompression in the inner layers. Here, γ ¼ K J ex as J ex , a and K are the exchange coupling constant, lattice constant and uniaxial anisotropy constant, respectively. Also, J ex ¼ aA 2S 2 where A is the exchange stiffness. Assuming the parameter A = 8 × 10 −8 erg/cm for Dy and Tb and K Dy = 7.5 × 10 6 erg/cm 3 and K Tb = 2.4 × 10 6 erg/cm 3 , the basal plane DW widths are ξ Dy~3 .2 nm and ξ Tb~5 .7 nm. Thus, a difference in the reduced number of spins at the end points due to unequal ξ along with structural interdiffusion ð¼ 2:0 ± 0:5 nmÞ at the interfaces of Dy and Tb layers would contribute to vertically uncorrelated domains when the two layers are stacked as a multilayer.
Based upon the temperature-dependent neutron-scattering intensity, a picture of the spin-glass state coexisting with the helical spin ordering was discussed earlier by Sato et al. 45 . An interdomain interaction should be pronounced when the competition of two types of modulated spin fluctuation becomes significant, particularly at the interfaces. A frustrated interdomain magnetic interaction would lead to the coexistence of spin-glasslike ordering and helical spin modulation.
Model interpretation of spin-glass as a short-range spin density wave (SDW) has been proposed earlier by Mydosh 17 . Therefore, magnetic configurations of topologically stable modulated helices can be visualized as a superposition of the double spin spirals or SDW within Dy and Tb with similar propagation wave vector q (=2π/N) and pitch, and separated by a phase factor. A resultant of superposition of two helices represented by the magnetic SLD profile is shown in Supplementary Fig. 8b, for clarity. The number of nodes/anti-nodes (modulation period) in the modulation of the envelope lines forming the double helices, which being commensurate with the chain's length, determines the number of twists in the chain. Increase in the number of twists increases the helix modulation, which in turn enhances the energy stored in the system 4,8 .

Discussion
We have demonstrated the realization of magnetic helices in form of 2π-planar DWs in an unprecedented multilayered system involving two hard rare-earth materials (Dy/Tb), instead of a more familiar soft-hard (Fe/Tb or Fe/Dy) combination. Using depth sensitive PNR, DWs have been shown to exist within both Dy and Tb. Since both materials are anisotropically hard and their respective helical and ferromagnetic ordering temperatures are also significantly different in bulk form, it was hardly expected for them to establish exchange coupled helices. Thus, our work realizes a novel state of spin configuration even within such unusual candidates.
The robustness of the helices at an applied field of -10.0 kOe (-1000 mT) ensures their topological stability. The helices within both RE/RE layers are stabilized by the strong in-plane magnetocrystalline anisotropy of the Dy and relatively softer Tb layers induced by the growth process. Exchange bias coupling of around -0.88 kOe (-88 mT) and -0.65 kOe (-65 mT) were revealed in two multilayers with different t Dy and t Tb thicknesses combinations, which can be owed to the enhanced pinning at spinfrustrated interfaces.
Rigorous magnetic data analysis suggests the existence of superspin-glass type of freezing of weakly interacting disordered spins. Additionally, such SSG type of behavior with T SSG ≃ 206 K, was found to be coexistent with the helical spin configurations. The dynamics in the superspin glass phase and its scaling behavior was further confirmed by TRM measurements. Such a coexistence of SSG and helical phases, is due to the spinfrustrated interdomain magnetic interaction as indicated by polarized off-specular neutron scattering.
We believe, that topological stability for the effective periodicity of the modulation and the internal effective field in the helix, which stores the magnetic energy density, can be interesting for exploiting their applicability in all-spin-based technology. The substrates were single crystalline Si wafers of 20 × 20 mm 2 , which were used after cleaning them in isopropyl alcohol. The targets were disks of 2 inch diameter and 0.125 inch thickness with a purity of 99.9%. They were bonded to a copper backing plate. The targets were cleaned by pre-sputtering for 10-15 min in Ar atmosphere. The deposition rates were about 0.215 nm/s for Tb and 0.282 nm/s for Dy and were precalibrated. During deposition, the Ar pressures in the magnetron sputtering chamber were 3.9 × 10 −3 mbar for Tb and 2.9 × 10 −3 mbar for Dy and the base pressure was 1.5 × 10 −8 mbar. The magnetron sputtered samples were expected to grow with a high degree of crystallographic orientation (texture) and in-plane easy axis as confirmed earlier from the XRD patterns (see Supplementary  Fig. 1) 8 .

Methods
X-ray characterization. XRR and XRD measurements were performed on an Empyrean diffractometer from PANalytical which provides information on the structure, thickness and interface roughness of the individual layers (see Supplementary Fig. 1). The individual layers and multilayers show broad polycrystalline peaks when investigated by XRD as observed earlier for similarly grown films 46 . Xray contrast between Dy and Tb is very poor due to their adjacent atomic numbers, hence a multilayer structure could not be properly studied using XRR.
Magnetometry. Conventional in-plane magnetization measurements were performed at various temperatures and fields using a superconducting quantum interference device (SQUID) from Quantum Design (MPMS-XL). Conventional AC field susceptibility measurements were acquired at various temperatures and frequencies using a physical property measurement system (PPMS) from Quantum Design.
Polarized neutron scattering. PNR measurements were performed at the reflectometer MARIA at MLZ, Germany 47 . In-plane magnetic fields of 10 kOe (1000 mT) and 0 Oe (0 mT) were used before the samples were cooled to 10 K separately in a cryostat and measured at applied fields H a = -0.01 kOe (-1 mT) and -10.0 kOe (-1000 mT). From the neutron polarization analysis we resolve different components of the magnetization within the film plane. The saturation field of approximately 50.0 kOe of the system exceeded the experimentally available field of 10.0 kOe (1000 mT) at MARIA. The SLDs of a specimen are given by the nuclear (ρ n ) and magnetic (ρ m ) components. Four different cross sections were measured, namely non-spin flip (NSF) scattering: (ρ n ± ρ m cos ϕ A ) represented by R ++ and R −− and spin flip (SF) scattering: (ρ 2 m sin 2 ϕ A ) represented by R +− and R −+ . Here, + and − signs are used to distinguish the intensity contributions R representing a polarization component parallel or anti-parallel to the guiding field, respectively. The angle ϕ A is between the magnetization M FM and the applied field H a . The spin dependent NSF reflectivities are given by while the SF reflectivities are given by Here, R + , R − are the Fresnel reflection amplitudes for neutron polarizations parallel and anti-parallel to magnetic induction of the sample. The amplitudes are given by the eigenvalues of the reflectance matrix which is diagonal at ϕ A = 0, π. The NSF reflectivities involve squares of the combinations of (1 − cosϕ A ) and (1 + cosϕ A ) terms. Thus, within the one dimensional analysis of the polarization vector it is not possible to discriminate the tilt angle ϕ A from (ϕ A + π) or the rotational sense of magnetization 40 . Note that by convention when the critical angle for R ++ precedes R −− along Q z , the ρ m value is regarded as negative.
Following the scattering geometry the momentum transfer vector Q can be written in terms of its various components: where the incident wave-vector defined by k i , makes an angle α i with respect to the x axis while the scattered wave-vector k f makes angle α f in the x-z plane and also χ in the x-y plane (relevant for diffuse scattering) and λ is the wavelength (here λ = 4.5 Å). The sample is located in the (x, y) plane. The Q y element can be ignored in this geometry as we do not resolve the intensities, but integrate, along that axis 40 . Figure 5a shows a sketch of the scattering geometry. The projection of the parallel component of magnetization onto the neutron polarization axis (y axis) is proportional to 〈cosϕ A 〉, while the projection of the component perpendicular with respect to the polarization axis onto the x axis is proportional to hsin 2 ϕ A i. Different length scales ξ ¼ 2π Q ranging from nm to μm can be accessed by using different scattering geometries in most practical cases.
The neutron coherence lengths 48 can be estimated from the uncertainties in the reciprocal space with the help of instrument resolutions. Thereby l x (along Q x ) and l y (along Q y ), which are contributions to the coherence length l || at the sample surface, turn out to be few μm and few å ngströms, respectively 40 . For domains (magnetic inhomogeneities) or roughness (nuclear inhomogeneities) larger than l || , the intensities reflected from different domains superimpose incoherently in the specular beam. In Fig. 5b, we have shown that if the domains are smaller than l || , then the neutron wave is scattered in specular and off-specular directions. If the domain correlations are larger ðξ k > l k Þ, then each domain reflect independently, which are then averaged over all possible domain orientations. Thus, the mean value of the orientation 〈cosϕ A 〉 can be obtained. The off-specular scattered intensity thus can be coherent (for periodic structure) or diffuse (for random distribution) with a spatial distribution of Δϕ A . This means that the domain magnetization is randomly tilted by angle Δϕ A with respect to the mean value averaged over the coherence volume along the perpendicular (〈sin 2 Δϕ A 〉) and parallel directions (〈cosΔϕ A 〉).

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The source codes used for neutron scattering fitting/simulation are available from the corresponding author upon reasonable request.