Introduction

The appearance of staircase-like structure is a fascinating phenomenon that is observed in a variety of condensed matter systems. In 2D electron gas, quantized conductance is manifested as a step feature in the Hall effect measurements1. In quantum materials, interplay of competing interactions with multiple periodicities in a system can give rise to a ground state whose length scales are defined by the modulation of the original periodicities. Example of such modulated periodicities includes commensurate and incommensurate phases, such as density waves in solids2, stripes and charge density waves in cuprate superconductors3,4,5, charge-ordered state in manganites6, and helical spin structure in magnetic systems7. A well-known staircase structure is the Devil’s staircase which appears when a system goes through numerous phase-locked modulated periodicities8,9,10. Devil’s staircase has been observed in magnetic systems8,11,12,13,14, liquid crystals15 and in ferroelectrics16. Apart from fundamental science, the staircase structures have potential technological applications such as in metrology, sensing devices etc17.

Interesting staircase structures in domain size and in magnetoresistance have been observed in Dzyaloshinskii-Moriya interaction (DMI) based solitonic system18,19 Competition between symmetric Heisenberg exchange interaction and the antisymmetric Dzyaloshinskii-Moriya interaction (DMI) can give rise to interesting magnetic textures such as a helix and skyrmion lattice phases11,20,21,22. DMI-based chiral magnetic order in a helimagnet is called Dzyaloshinskii type helimagnet structure, while a helical magnetic order due to competition between ferromagnetic and antiferromagnetic exchange interaction is known as Yoshimori-type helimagnetic structure23. The chiral magnetic structures in a helimagnet exhibits solitons that can be manipulated by an external magnetic field24. More specifically, the soliton periodicity changes in a stepwise manner which is attributed to the discrete changes in the soliton number because of confinement at the grain boundaries22,24. Field evolution of confined helicoids has also been shown to occur via discrete steps in helical magnet MnSi25. The thin film structure of MnSi accommodates a finite number of turns and the jumps are explained by the annihilation of individual turns of the helicoid.

In this article we report the observation of a staircase-like structure in the field-evolution of the scattering wave vector Q which emerges from the stripe phase of an amorphous and centrosymmetric Fe/Gd system. In contrast to bulk DMI magnets, previously described, the Fe/Gd system is a perpendicular magnetic anisotropy (PMA) system that exhibits dominant dipolar interactions and negligible DMI26,27. Given the centrosymmetric nature of the Fe/Gd system, magnetic phases will exhibit an equal distribution of opposite helicity magnetic spin textures that on average make the dipolar magnet globally achiral. We performed resonant coherent soft X-ray scattering to study the evolution of the stripe periodicity as a function of applied perpendicular magnetic fields at various temperatures. Under applied perpendicular fields and temperature conditions it is possible to obtain a hexagonal skyrmion lattice phase that consists of an equal population of left- and right- chirality27. We observed that the scattering wave vector Q changes in step-like fashion with no well-defined step height and width. However, when Q is resolved into components Qx and Qy, the step heights along Lx = 2pi/Qx, were found to be in integer multiple of 7 nm, which is close to the exchange length of the system. At higher temperatures, the steps were smeared due to thermal fluctuations.

Our X-ray scattering studies have been complemented by spin dynamics calculations that take into account the net achiral nature of the material. We have simulated an experimentally observed (non-equilibrium) process where a global versus local phenomena delicately balances each other. On one hand, the total periodicity of the stripes increase with increasing magnetic field due to an enhancement in the majority stripe width (spins aligned along the field direction). On the other hand, this is being counter balanced by the local minority stripe width which cannot fall below a certain size governed by the exchange length of the system. Thus, instead of a bulk macroscopic motion of domain walls over the entire sample, these competing tendencies cause the local forces and energetics experienced by the minority domains to locally annihilate some of these half-periods. This in turn leads to (as observed experimentally and verified theoretically) a local readjustment of the domain sizes. By defining a magnetic domain length scale ratio, T, we have performed a Landau-Lifshitz (LL) simulation to generate steps as a function of the applied magnetic field and show the important role that anisotropy plays in generating the steps in these systems. We have developed a theoretical model for the appearance of steps using exchange interaction, dipole-dipole interaction, anisotropy, and external magnetic field. Our calculations indicate that the origin of the steps lie in the anisotropy term. Even if exchange and dipole interactions are present, the absence of anisotropy does not produce steps. Although the appearance of steps do look similar in single crystal DMI material and the amorphous Fe/Gd system under investigation, the physical origin of the steps in the two systems is different.

Results

Resonant scattering due to stripes

The scattering geometry for the experimental set-up is shown in Fig. 1(a). X-ray beam whose energy is tuned to the Fe L3 edge (707 eV) is incident normally on the sample. A pinhole was placed on the beam path upstream 5 mm from the sample to establish transverse coherence of the beam. In this geometry the X-ray photons are sensitive to the magnetization (mz) along the beam direction. The scattering pattern was collected on a charge coupled device camera (CCD) placed about 0.5 m away downstream. Resonant X-ray scattering measurements are sensitive to static magnetic structure (S(q)) and spatial correlation length (ξs). From the position and intensity of the Bragg peaks it is possible to extract information about the periodicity and strength of the magnetic order. Fig. 1(b) shows the full field X-ray microscope images (top panel) and X-ray resonant scattering pattern (bottom panel) of the sample. We observed the presence of three distinct magnetic phases, namely, disordered stripe, ordered stripe and skyrmions, obtained by either varying the temperature or applied perpendicular magnetic field. The X-ray microscopy images were obtained by varying the applied magnetic field at 300 K, while the resonant X-ray scattering data was measured from LN2 temperatures to room temperature as a function of the applied magnetic field. In the ordered stripe phase (T = 239K) the domain periodicity (2π/Q) at remanence is (119 ± 5) nm. The stripe pattern persists as the field is increased from zero to around 170 mT, when new peaks in the form of a distorted hexagonal pattern start to appear indicating a transition to the skyrmion phase. These observations are consistent with previous findings26,27,28.

Fig. 1: Experimental set-up and magnetic phases.
figure 1

a Schematic of the coherent magnetic X-ray scattering geometry. b Real space (top panel) and reciprocal space (lower panel) of the different magnetic phases present in a Fe/Gd system. The scattering images were taken at H = 0 mT; T = 85 K (disorder stripes), H = 0 mT; T = 225 K (order stripes) and H = 190 mT; T = 239 K (skyrmions). A small residual in-plane field is present during ramping down of field from saturation to zero. c Variation of the 1st and 2nd order magnetic diffraction peak with field at 239 K. The Inset image shows the appearance of both 1st and 2nd order diffraction peaks.

Full size image

Figure 1 (c) shows the field evolution of the integrated intensity of 1st and 2nd order diffraction peak from the ordered stripe domains. Starting from the zero magnetic field condition the 1st order peak is a maximum and the 2nd order peak is a minimum (~zero). At remanence, the average value of the out-of-plane component of the magnetization is zero. Therefore the widths of the up and down domains are equal. This results in odd-order diffraction spots. Increasing the applied perpendicular field breaks the symmetry between the diffraction peaks, from the stripe phase, causing the even order diffraction peaks to appear29. Around 170 mT, the intensity of both the 1st and the 2nd order peaks start to diminish and eventually, new peaks in the form of a hexagonal ordering pattern appear (see Fig. 1(b), bottom right panel). It is interesting to note that in the hexagonal phase, we observe two relatively strong intensity spots along the same direction where the stripe peaks are present which would indicate that somehow the original direction of the stripes is retained even in the hexagonal phase.

Staircase structure of Q-vector

The evolution of the stripe-diffraction peak in Q-space is shown in Fig. 2(a) as a function of the applied perpendicular field at T = 230 K. At the start of the field cycle, the momentum transfer vector, Q1 is (= 0.052 nm−1) of the magnetic Bragg peak. As the field increases, the magnetization increases and the size of the favorable perpendicular domains (spins along the field directions) also increases leading to an increased domain periodicity which results in a decrease in the Q-value 2π/L0, where L0 is the periodicity with respect to zero field Q-value. Interestingly, we observed that the Q-value corresponding to the magnetic Bragg peak decreases in discrete steps as a function of applied magnetic field giving rise to a staircase-like structure.

Fig. 2: Evolution of stripe diffraction peak, periodicity and correlation with field.
figure 2

a Plot of the q-vector of the satellite peaks as a function of the applied out-of-plane (OOP) magnetic field at 230 K as the system transitions from magnetic stripe phase to skyrmion phase. Dotted arrows indicate the Q-value of the Bragg peak positions (purple symbol) starting from at different fields. b Evolution of the stripe-periodicity with field at various temperatures showing discrete steps like feature. c Correlation coefficient values with respect to the remanent state (0 mT) for increasing magnetic field at 85K.

Full size image

The evolution of domain periodicity obtained from scattering data happens in several steps that involve sudden jumps and appearance of modulated periodicities. We find that along with the main magnetic Bragg peak, a much weaker satellite peak forms at a smaller Q-value, and both the peaks evolve in an interesting way as the field is changed. The increase in field leads first to the appearance of an initially weaker intensity satellite peak at Q2 (at a smaller Q-value than zero field Bragg peak at Q1). With further increase in the field, the main Bragg peak (Q1) suddenly merges with Q2 giving rise to a step-like feature in Fig. 2(a). Since the position and intensity of the Bragg peak gives the periodicity of the stripe domains and density of domain scattering respectively, we can conclude that the number of domains with periodicity P1 = 2π/Q1 decreases with increasing field while the number of domains of periodicity P2 = 2π/Q2 starts to increase and finally all domains suddenly transform to the periodicity P2. This sequence of events, changing Q from Q1 to Q7 with a similar mechanism of peak shifts (Q1 → Q2; Q3 → Q4; Q5 → Q6) was observed throughout the stripe phase (see Fig. 2(a)). In some cases, a direct change in Q-values corresponding to the Bragg peaks without any satellite peaks (Q2 → Q3; Q4 → Q5) was also observed.

In Fig. 2(b) we convert the wavevector into real space periodicity (2π/Q) and plot it as a function of applied field at different temperatures. At higher temperatures the total number of steps increases which results in the appearance of the first step at much lower fields for higher temperatures than the lower ones. The plot of the correlation values of the stripe-diffraction spot at different fields with respect to the one at remanence for increasing magnetic fields is shown in Fig. 2c. Any subtle changes in the speckle pattern between two frames taken at 0 mT and H mT will result in a value of the correlation coefficient (CC) which is defined by

$${{{\boldsymbol{CC}}}}=\frac{{\sum }_{m}{\sum }_{n}({A}_{mn}-\bar{A})({B}_{mn}-\bar{B})}{\sqrt{({\sum }_{m}{\sum }_{n}{({A}_{mn}-\bar{A})}^{2})}\sqrt{({\sum }_{m}{\sum }_{n}{({B}_{mn}-\bar{B})}^{2})}},$$
(1)

where A and B correspond to the two images taken at two different field values. Amn denotes the intensity value of the pixel position at mth row and nth column of the 2D scattering image. \(\bar{A}\) is the mean value of the 2D image. If CC = 1 then the two images are perfectly correlated, CC = 0 means completely de-correlated and CC values lying between 0 and 1 mean partially correlated. Thus the variation of the correlation-coefficient can be attributed in the real space as either a change in magnetization or density or periodicity of the stripes or any combination of these factors with the applied field.

In Fig. 2c the CC for the stripe phase observed at 85 K is calculated between the scattering image taken at remanence (zero-field) and another image taken at high fields. The field-dependent variation of the CC is plotted in Fig. 2c (blue color line). The correlation coefficient also exhibits distinct steps and most interestingly the horizontal portion of the steps are populated with small step-like features identical to a self-similarity devil’s staircase like behaviour. As a measure of the stability of the entire set-up during the measurement we also calculated the correlation coefficients for the Airy pattern (green color line), which remains fairly close to unity at all the fields.

Resolving staircase along Qx and Qy direction

A typical diffraction pattern consisting of the centrosymmetric first order peaks in the stripe phase is shown in Fig. 3(a) along with their in-plane Q-vectors. The enlarged image of the diffraction spot in Fig. 3(b) exhibit modulation with speckles indicative of heterogeneity in the ordering of the magnetic domains in the real space. The diffraction spots appears at about 45∘ to the beam propagation direction (see Fig. 3(a)), meaning the stripe-domains are oriented 45∘ to the X-ray propagation direction. We note here that a small in-plane field of magnitude ≈ 1mT is present along y-direction (see Fig 1(a)) in the magnet during rampdown from saturation to zero fields. We resolved the resultant Q-vector into Qx and Qy components, to get information about the stripe periodicity along real space X and Y direction thereby obtaining real space value Lx and Ly as shown in the schematic representation in Fig. 3(c). From this we find that the steps along Lx are significantly distinct compared to that in the Ly direction (see Fig. 3(d and e)).

Fig. 3: Stripe orientation and staircase-like behaviour.
figure 3

a A typical scattering pattern of the stripe lattice along with the projection of the in-plane Q-vectors. b Enlarged image of the stripe-diffraction spot in Q-space. c Schematic real space view of stripe-domain orientation according to scattering image of Fig. 3a, where blue circles with dot resemble the spin along the field direction while the red small circles with cross resemble the spins opposite to the field direction and \({{{\rm{L}}}}={{{{\rm{L}}}}}_{x}^{2}\)+L\({}_{y}^{2}\) corresponds to the periodicity of the stripe domains. Plot of the evolution of d Ly (= 2π/Qy) at 85 K and e–g Lx (= 2π/Qx) as a function of the applied magnetic field at T = 85 K, 183 K and 236 K.

Full size image

Interestingly, we found that the steps in the plot of Lx change in multiples of 7 nm. That is, the minimum change in periodicity along Lx is 7 nm. No such relationship exists for the Ly evolution. The magnitude of the change in Ly is typically smaller than 7 nm and random compared to Lx (Fig. 3(d and e)). A schematic of a possible stripe domain arrangement is shown in Fig. 3(c). The blue (red) domains are majority (minority) domains. The stripes are slanted with respect to the applied perpendicular field direction (z). We know from the experimental results that as the applied field is increased, the Q-vector of the magnetic Bragg peak moves to a lower value, but maintains its orientation of 45∘ with respect to beam direction. This indicates that the majority domain expands but the stripe-domains maintain the 45∘ orientation. Step changes in multiples of 7 nm along Lx imply that the x-component of the periodicity (2π/Qx) changes in units of 7 nm.

Interestingly, this value matches with the exchange length (Lex) of the studied Fe/Gd system. One of the ways to think about this behavior is that as the minority domains shrink (the horizontal region of the steps), there is a minimum distance between successive domain walls below which there cannot be a smooth deformation of the spin texture, as a result, a sudden jump happens. In the theoretical section we will show that indeed by defining a term that signifies the ratio of spin kink to the spin chain, it is possible to predict jumps. At higher temperatures, we observed an increase in the number of steps in the average-periodicity curves (see Fig. 3(f and g) as a function of the applied perpendicular field. This is due to the fact that thermal fluctuations aid in a faster transition from one step to the other as a result we obtain more steps at 236 K than 85 K even though the field range over which such steps occur is much higher at lower temperatures.

The existence of steps in solitonic systems with DMI has been observed experimentally and explained theoretically19,22. The presence of DMI introduces a topologically protected kink in the spin texture. The topological protection of the kink means that there is an energy cost to kink annihilation. Different topological sectors have different energy which is the reason for step-like features. In contrast, in the Fe/Gd systems, the dominant interactions are exchange, dipole, and anisotropy. This supports an achiral magnetic structure. So far there have been no theoretical studies of the step-like behaviour on dipole interaction dominant achiral spin-structures in an amorphous system. In the theoretical model presented in the next section, we have mimicked the experimental conditions by investigating a one-dimensional dipolar mediated spin chain which is achiral in nature. We have numerically solved the LL equation of motion to understand the magnetization dynamics observed in the Fe/Gd system experiment. Based on our calculations we show that the origin of the step-like behaviour under the application of an external perpendicular magnetic field could be explained by the spin dynamic behavior of an achiral spin chain.

Model and theory

The spin kinks caused by long-range dipolar interaction in the Fe/Gd system can be classified by a number n. In Fig. 4 we show the local spin arrangement in a finite-size chain under zero applied magnetic field with fixed boundary condition on both ends. We consider a N-site 1D chain where spins interact with exchange interaction, dipolar interaction, anisotropy, the in-plane, and the out-of-plane (perpendicular) magnetic field. The spin on each site is parameterized as

$${{{{\boldsymbol{S}}}}}_{i}=(\sin {\theta }_{i}\cos {\varphi }_{i},\sin {\theta }_{i}\sin {\varphi }_{i},\cos {\theta }_{i}),$$
(2)

where the site spin angle \({\varphi }_{i}=2\pi ni/{{{\mathcal{N}}}}\) and \({\theta }_{i}=\frac{\pi }{2}\). Here \(i=0,1,2,...,{{{\mathcal{N}}}}\) where \(N={{{\mathcal{N}}}}+1\). The kink sectors are classified by n which indicates the number of domains existing in the chain. The Hamiltonian for our Fe/Gd system is

$$H={H}_{J}+{H}_{D}+{H}_{K}+{H}_{h},$$
(3)

where the meaning and expression of each term is given by

$${H}_{J}=-J\mathop{\sum}\limits_{i\in {{{\mathcal{N}}}}}{{{{\boldsymbol{S}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{i+1}({{{\rm{exchange}}}}),$$
(4)
$${H}_{D}=D\mathop{\sum}\limits_{i,j\in {{{\rm{sc}}}}}{{{{\boldsymbol{S}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{j}{\Pi }_{ij}\,({{{\rm{dipolar}}}}\,{{{\rm{interaction}}}}),$$
(5)
$${H}_{K}=-{K}_{U}\mathop{\sum}\limits_{i}{({{{{\boldsymbol{S}}}}}_{i}\cdot {{{\boldsymbol{x}}}})}^{2}\,({{{\rm{anisotropy}}}}),$$
(6)
$${H}_{h}=-g{\mu }_{B}{H}_{x}\mathop{\sum}\limits_{i}{S}_{i}^{x}-g{\mu }_{B}{H}_{y}\mathop{\sum}\limits_{i}{S}_{i}^{y}\,({{{\rm{magneticfield}}}}).$$
(7)
Fig. 4: Achiral spin chain arrangement.
figure 4

Achiral magnetic order is generated due to the competition between exchange interaction and dipolar interaction. The magnetic texture shown in the cartoon depicts an achiral spin arrangement where the spins rotate 180∘ out of plane (say in the positive screw direction, z) and then rotate back to the original position (with an opposite negative screw rotation, −z). The spatial distance over which the achiral twist occurs can be captured by defining a local coefficient T = Ld/Lc, where Ld is the length of the achiral domain area and Lc is the chain length. The cartesian coordinate system shows the definition of the angle φi and the angle θ used in the simulation.

Full size image

In the above i either denotes the lattice site in the 1D chain or the location of a spin site inside a supercell (sc). The exchange interaction strength is given by J > 0, the dipolar interaction coupling by D, the anisotropy by KU, and the out- and in- plane magnetic field is given by Hx and Hy, respectively. The symbol g denotes the gyromagnetic ratio and the μB is the Bohr magneton. The Πij in the dipolar interaction term is the Ewald coefficient which captures the long-range nature of the dipolar interaction. Using the angular representation of the spin Si we can write the total energy as

$$\begin{array}{l}\frac{H}{J{S}^{2}}\,=\,-\mathop{\sum}\limits_{i\in N}\cos ({\varphi }_{i+1}-{\varphi }_{i})+{J}_{d}\mathop{\sum}\limits_{i,j\in \,{{\mbox{sc}}}\,}{\Pi }_{ij}\cos ({\varphi }_{i}-{\varphi }_{j})\\ \qquad\quad\,\,-K\mathop{\sum}\limits_{i}{\cos }^{2}{\varphi }_{i}-{h}_{x}\mathop{\sum}\limits_{i}\cos {\varphi }_{i}-{h}_{y}\mathop{\sum}\limits_{i}\sin {\varphi }_{i},\end{array}$$
(8)

where we have now introduced the scaled variables \({J}_{d}=\frac{D}{J}\), \(K=\frac{{K}_{U}}{J}\), \({h}_{x}=\frac{g{\mu }_{B}{H}_{x}}{J}\) and \({h}_{y}=\frac{g{\mu }_{B}{H}_{y}}{J}\). In all our figures we will report the scaled fields in milli-units, that is, hx = 1 stands for 10−3 scaled field units.

We implement the local achiral spin structure N = 216 sites, shown in Fig. 4, to perform the LL simulation. To mimic the finite size of the experimental sample and to allow for the domains to grow and collapse as observed experimentally, we utilized an embedding trick to simulate the LL equations-of-motion (EOM). To capture experimentally realistic sample conditions, from a computational perspective, we introduced the concept of a local coefficient T. From a physical perspective, T represents the ratio of the length of the achiral structure (which contains the twist sectors solely) over the length of the 1D chain. Thus, the achiral structure is embedded within a uniform ferromagnetic background spin texture. The computational embedding trick allows us to capture the spontaneous rearrangement of the twist sectors in the chain configuration, thereby simulating the growth and collapse of the achiral domain walls. Our numerical simulations indicate that the eventual fate of the twist sectors and subsequent realization of jumps (as observed experimentally) is a subtle balance between Jd, K, and N. We compute the minimum energy Emin using Eq. (8). The magnetization M is calculated using

$$M=\frac{1}{N}\mathop{\sum }\limits_{i=0}^{{{{\mathcal{N}}}}}\cos {\varphi }_{i}.$$
(9)

We present the energy and corresponding magnetization response of the local achiral state in Fig. 5. When the anisotropy is absent, we observe that the energy is degenerate for different twist sectors and no jumps are created by enhancing the dipolar interaction (see Fig. 5(a)-(c)). Moreover, it indicates that larger dipolar parameters induce a downshift in energy with no visible effects on the magnetization behavior in the local achiral state. In the presence of anisotropy, we keep the dipolar interaction constant and increase the K parameter as shown in Fig. 5(c)-(e). We compute the LL dynamics on a chain of local achiral state with different anisotropy parameters. We find that upon enhancing anisotropy in the presence of a magnetic field, the energy degeneracy of the different twist sectors is broken with a simple upshift. With a relatively small \(T=\frac{1}{4}\) and a strong enough anisotropy K = 0.2, we observed jumps in both energy and magnetization in response to magnetic field as shown in Fig. 5(e).

Fig. 5: The energy and magnetization response for different T values.
figure 5

In the upper panel, dipolar parameter is Jd = 0.00934 and the anisotropy parameter is K = 0.2 (c). Jd = 0.00962, 0.01004 in b and a while anisotropy parameter K remains the same as the one in c. Anisotropy parameter K = 0.1, 0.2 in d and e while dipolar parameter remain the same as the one in e. In the lower panel, Jd = 0.00962 and K = 0.2. The local coefficients from the left to the right are \(T=\frac{7}{8},\frac{7}{9},\frac{2}{3},\frac{1}{2},\frac{1}{3}\), respectively. Red arrows in h–j represent the first jump for energy curve with n = 5. The corresponding magnetic fields are hx = 10.3, 9.2, 6.9. All results are calculated with the number of sites N = 216.

Full size image

In Fig. 5f–j we show our calculations of energy and magnetization response as T is varied. With decreasing T, jumps begin to happen in energy response with higher twist sector. When \(T > \frac{1}{2}\), jumps happen in energy curves with n = 6 as shown in Fig. 5(f)-(h). However, jumps happen in energy curves with smaller twist sectors n and lower magnetic field intensity hx. In both Fig. 5i and Fig. 5j, jumps happen when twist sector is n ⩾ 4. And the critical magnetic field intensity for the first jumps to happened decreases as T decreases. It is found that energy response is more powerful to show the disappearance of kinks while the jumps in magnetization response might be caused by the position shifting of the kinks.

We have considered a chain with larger number of sites. When the number of sites is N = 432 and the dipolar parameter Jd = 0.00916, jumps can be observed in energy curve with twist sector n = 4 and local coefficient \(T=\frac{1}{4}\). However, no jumps can be observed with \(T=\frac{1}{3}\) and \(\frac{1}{2}\) (plots not shown). This behavior can also be seen in a system with N = 864. The result that the decreasing T contribute to the jumps, is also consistent with the N = 216 system. Moreover, when the Jd increases, more kinks are able to establish and more jumps are observed. Thus, we draw a conclusion that not just the declining local coefficient T, but also the rising dipolar parameter Jd results in the jumps happening in energy curve with smaller twist sector n and weaker magnetic field hx. Note, that for particular parameters J, D, K and the number of sites N, the value of the ground state energy remains almost the same when the local coefficient T changes.

Discussion

In this work, we have shown experimentally that in an amorphous and centrosymmetric Fe/Gd magnetic thin film that exhibits stripe and skyrmion lattice phases, the stripe-domain periodicity changes in steps because of the abrupt disappearance of stripe-domains. This result is interesting in itself because similar to a DMI-based solitonic system, exchange-dipole mediated Fe/Gd system also shows similar step-like behavior even though a sole global chirality is absent in the system. Since the presence of a global DMI can be ruled out in the Fe/Gd systems26,27, we can conclude that the predicted spin twists are formed due to the competition between exchange and dipolar interactions, and the spin twist sectors can be smoothly transformed to the uniform phase by any number of finite deformations.

Intuitively, due to the achiral spin texture of the domains, as the magnetic field is increased, the minority domains start to shrink, resulting in two “like-domains" to come closer. The minimum distance between the two domains is guided by two spin-kinks on either side which should be equivalent to the length of two domain walls. Using the well-known formula \({l}_{w}=\sqrt{J/K}\), where lw is the domain wall width, the domain wall width for Fe/Gd comes out ≈3.2 nm, twice of which is 6.4 nm, which is in close agreement to the experimental value of 7 nm. Thus the minimum distance between the two like-domains comes out to be equivalent to the exchange length (Lex) of the system from our experimental study. The above explanation also points to the existence of a “global" and “local" length scales in the system which will give rise to two energy scales. It is these competing energy scales that give rise to steps. Our system is reminiscent of the case of modulated periodicities.

Our 1D model suggests the existence of discrete magnetization steps in the REXS experiment results from magnetic spin textures exhibiting achiral spin-twist characteristics. Recently, a transport and micromagnetic study30 of a patterned Fe/Gd specimen, with the same material composition as the one addressed in this work, suggests the domain walls of the stripe domains undergo a local chirality spin rearrangement from chiral to achiral, under similar applied perpendicular field conditions, which results in stripe domains exhibiting achiral spin-twist characteristics as the one addressed in the 1D model of this paper. Although Ref. 30 addresses the formation mechanism of skyrmion lattice phases in the Fe/Gd system specimen, close inspection of the field-dependent micromagnetic vector magnetization (mx, my, mz) evolution, in their work, shows the magnetization components attributed to the domain wall (mx, my) undergoes abrupt/sharp changes as the perpendicular field is swept from zero-field towards magnetic saturation which is likely correlated to the collapse and local rearrangement of magnetic spin textures with achiral spin-twist characteristics.

We take the achiral nature of the stripe spin structure as an important point in our theoretical development and show that the magnetization steps can indeed be observed in a dipolar magnet with net global achirality. The variations and interplay of the length scales is captured in the parameter T. Analysis of the energy expression with different values of T suggests that in an achiral spin arrangement staircase structure can be observed only under certain specific ratios of 1D spin-to-spin-twist length scale. Although simplistic, our LL calculations using local achiral spin structure shown in Fig. 4 is able to capture the essential feature that the system has jumps in response to an external magnetic field. The jumps happen only when anisotropy is present. The absence of anisotropy leads to a degeneracy of energy response for different twist sectors, meaning the absence of jumps in the system. Our study provides evidence and further impetus to study magnetic spin textures in a centrosymmetric magnetic, both from an experimental and theoretical viewpoint.

Methods

Sample details

The samples studied were nominally [Gd (0.4 nm)/Fe(0.34 nm)] × 80 multilayers deposited using DC magnetron sputtering with 20 nm Ta seed and capping layers. The samples were deposited on 50-nm or 200-nm thick Si3N4 membranes to allow for transmission RSXS experiments, respectively. Non-resonant 12 keV x-ray diffraction indicated strong intermixing of Gd and Fe layers thereby forming an amorphous structure rather than a multilayer. Given the centrosymmetric nature of the Fe/Gd system, negligible Dzyaloshinskii-Moriya exchange interactions are expected27,28.

Experimental details

The coherent X-ray magnetic scattering measurements were performed at beamline 12.0.2.2 of the Advanced Light Source, LBNL. The incident beam was tuned to the Fe L3 edge (707 eV). Transverse coherence of the X-ray beam was established by inserting a 10 μm pinhole in the beam path before the sample. The scattering experiment was done in the transmission geometry at temperatures ranging from 40 K to 300 K as a function the perpendicular magnetic field from 0 mT to 500 mT. (Fig. 1(a)). The sample was subjected to the following initial magnetic field protocol. First the field was raised to 500 mT, then lowered to −500 mT and finally to zero before taking the measurements. The field ramp rate for the first two legs is 13 mT/sec, while the final drop of field from -500 mT to 0 mT took place at a rate of 380 mT/sec. We start our measurement at this zero-field condition and proceed to measure the diffraction signal as a function of applied magnetic field at a constant rate of 1.575 mT/s. A Charge Coupled Device (CCD) camera placed at about 0.5 m downstream of the sample was used to record the scattered intensity patterns.

Ewald method

In the Hamiltonian calculation, Ewald summation is applied, which is given by

$$\begin{array}{l}{\Pi }_{ij}=\sqrt{\frac{2}{\pi }}\frac{1}{3{\sigma }^{3}}\mathop{\sum}\limits_{n}{e}^{-\frac{{\left\vert {r}_{ij}-{\mathsf{n}}L\right\vert }^{3}}{2{\sigma }^{2}}}+\frac{4\pi }{\Omega }\mathop{\sum}\limits_{k\ne 0}{e}^{-\frac{{\sigma }^{2}{k}^{2}}{2}}\cos (k{r}_{ij})\\ \quad\;\;-\,\sqrt{\frac{2}{\pi }}\frac{1}{3{\sigma }^{3}}{\delta }_{ij},\end{array}$$
(10)

where rij represents the distance between two spin sites, L is the size of the supercell, n is the supercell label, σ is the real-space cut-off, k is the momentum space label, and Ω is the volume of the supercell which in our case is equal to L. The Ewald parameter will be redefined as Πij = Π∣i−j∣ ≡ Πm, where the symbol m tracks the number of Ewald parameters for the specific supercell size choice. The values of Πm are shown in Table 1.

Table 1 Ewald parameter symbols and the corresponding Ewald parameter values given by Eq. (10).
Full size table

Landau-Lifshitz (LL) equation of motion

We perform a LL EOM spin dynamics calculation on Eq. (3). We obtained an iterative equation which can be used to calculate the angle φi of each spin on the chain. Based on our computations, we are able to generate a stabilized spin order along the chain. Next, we computed the twist angle Δφ of the ground state in the absence of an external magnetic field using the energy minimization condition for the supercell given by the expression

$$\frac{E}{J{S}^{2}}=-\cos \Delta \varphi +{J}_{d}\mathop{\sum }\limits_{m=1}^{L-1}(L-m)\cos (m\Delta \varphi ){\Pi }_{m}.$$
(11)

The angle φi was analyzed to obtain the relationship between the range of Jd and the number of kinks for a given lattice size of site N. The relationship between the number of sites N, dipolar parameter Jd and the maximal sector nmax is shown in Table 2. Note that to perform the calculation, one needs to choose a supercell size that stabilizes the ground state and ensures that there will be minimal to no numerical oscillations in the computed result due to convergence issues. We found that L = 6 is the optimal supercell size which yields numerically stable results for our LL analysis. Using the numerically stable data, we computed the minimum energy Emin (scaled relative to JS2) and magnetization M (scaled relative to S). To compare our numerical results with the experimental setup of the Fe/Gd system, we need to mimic the experimental conditions. Therefore, all the results are calculated by applying a tiny in-plane field hy.

Table 2 Maximum allowed number of kink sectors nmax under different scaled dipolar parameter Jd for a given number of lattice sites N.
Full size table

The two angular variable EOMs are given by (for ℏ = 1)

$$S\sin {\theta }_{i}{\partial }_{t}{\theta }_{i}=\frac{\partial H}{\partial {\varphi }_{i}},\quad S\sin {\theta }_{i}{\partial }_{t}{\varphi }_{i}=-\frac{\partial H}{\partial {\theta }_{i}},$$
(12)

where only the first equation is required because the angle θ is held constant. Using Eq. (12) we can obtain the following expressions

$$\begin{array}{ll}0=\frac{1}{2}\left[\sin ({\varphi }_{i}-{\varphi }_{i-1})-\sin ({\varphi }_{i+1}-{\varphi }_{i})\right]\\ \quad\;\;\;+{J}_{d}\mathop{\sum}\limits_{j\,\in \,\,{{\mbox{sc}}}\,}\left[\sin ({\varphi }_{i+| i-j| }-{\varphi }_{i})-\sin ({\varphi }_{i}-{\varphi }_{i-| i-j| })\right]{\Pi }_{ij}\\ \quad\;\;\;+2K\sin {\varphi }_{i}\cos {\varphi }_{i}+{h}_{x}\sin {\varphi }_{i}-{h}_{y}\cos {\varphi }_{i}.\end{array}$$
(13)

The above can be split further into a form convenient for a numerical iterative self-consistent approach to solve for the angle ϕi. Hence, we write

$${A}_{i}=\frac{1}{2}(\sin {\varphi }_{i+1}+\sin {\varphi }_{i-1})-{J}_{d}\mathop{\sum}\limits_{j\,\in \,sc}(\sin {\varphi }_{i+| i-j| }+\sin {\varphi }_{i-| i-j| }){\Pi }_{ij}-K\sin {\varphi }_{i}+{h}_{y},$$
(14)
$${B}_{i}=\frac{1}{2}(\cos {\varphi }_{i+1}+\cos {\varphi }_{i-1})-{J}_{d}\mathop{\sum}\limits_{j\,\in \,sc}(\cos {\varphi }_{i+| i-j| }+\cos {\varphi }_{i-| i-j| }){\Pi }_{ij}+K\cos {\varphi }_{i}+{h}_{x},$$
(15)

with the site angle φi is defined as

$$\sin {\varphi }_{i}=\frac{{A}_{i}}{\sqrt{{A}_{i}^{2}+{B}_{i}^{2}}},\quad \cos {\varphi }_{i}=\frac{{B}_{i}}{\sqrt{{A}_{i}^{2}+{B}_{i}^{2}}}.$$
(16)

The LL equation is solved with the boundary condition φ0 = 0 and \({\varphi }_{{{{\mathcal{N}}}}}=2\pi n\). To obtain the final spin structure, we need to do the following things. First, we set the spin structure as a local achiral structure in which ratio of the length of the domain area over the chain length equals to T (see Fig. 4). Second, we setup the dipolar and anisotropy parameters and run the LL simulation program to compute the stabilized spin structure. Finally, we change the magnetic field to compute the system’s magnetization and energy response to the magnetic field.