Magnetostatic twists in room-temperature skyrmions explored by nitrogen-vacancy center spin texture reconstruction

Magnetic skyrmions are two-dimensional non-collinear spin textures characterized by an integer topological number. Room-temperature skyrmions were recently found in magnetic multilayer stacks, where their stability was largely attributed to the interfacial Dzyaloshinskii–Moriya interaction. The strength of this interaction and its role in stabilizing the skyrmions is not yet well understood, and imaging of the full spin structure is needed to address this question. Here, we use a nitrogen-vacancy centre in diamond to measure a map of magnetic fields produced by a skyrmion in a magnetic multilayer under ambient conditions. We compute the manifold of candidate spin structures and select the physically meaningful solution. We find a Néel-type skyrmion whose chirality is not left-handed, contrary to preceding reports. We propose skyrmion tube-like structures whose chirality rotates through the film thickness. We show that NV magnetometry, combined with our analysis method, provides a unique tool to investigate this previously inaccessible phenomenon.

The scale bar is 20 µm for panels a-e, 500 µm (500 nm) for panel f (its inset). For the regime t ≪ d we use equations (11) and (13), while the full expression for α x,y , α z is retained in all the other cases. We clearly see that the bubble-like feature has disappeared at this field. Note also that the stray field B z is qualitatively constant while moving along the boundary of the magnetic disc, supporting the assumptions that at these fields the magnetization is mostly out-of-plane. c Surface topography recorded by monitoring the piezo voltage V z . The dashed black boundary is the same as in a and b, stressing the fact that the region in which we observe a magnetic stray field is smaller than the actual physical disc size. The scale bar is 500 nm for all panels. the two different bias fields of 9.5 (11.8) mT. As the NV depth increases the noise in m z also increases due to the exponential prefactor in eq. (7). Scale bar is 400 nm.  where fractional values are obtained by linear interpolation of the data in d. f Field required to stabilize a R = 50 nm radius skyrmion in the multilayer as a function of DMI. For large DMI, i.e., once all layers have a left-handed chirality, the field scales linearly with D i , as indicated by the red line. For D i < 1 mJ/m 2 , the required field is larger than expected from a material with uniform magnetization along z.

Supplementary Note 1. Sample preparation
In our experiment, magnetic discs were deposited at the top of cleaved quartz fiber ∼ 50 µm in diameter. We employ Sutter Instruments quartz rods with the initial diameter of 1 mm. The fiber diameter is controlled via a laser-based puller (P-2000 Sutter Instruments). Fabrication proceeds with a controlled mechanical cleaving of the pulled fiber after the introduction of an intentional crack by means of a diamond scribe. The resulting tip, shown in Supplementary  Our experiments were performed with a type IIa diamond grown by chemical vapor deposition by Element 6 measuring 4x2x0.05 mm 3 . We studied NV centres formed by 15 N ion implantation at an energy of 18 keV and a density of 500/µm 2 and subsequent annealing for 2 hours at 800 • C. This implantation energy is expected to yield NV centres at an estimated ∼30 nm depth from the diamond surface 1 .
For the diamond pillars we first prepared an etch mask, patterned on the diamond via electron beam lithography using a FOX 16 flowable oxide resist from Dow Corning. Adhesion of the resist was guaranteed by a very thin, ∼ 10 nm layer of Titanium deposited via electron beam evaporation. The exposed etch mask pattern was transferred onto the diamond via a conventional top down anisotropic plasma etch performed in a Unaxis Shuttleline ICP reactive ion etching (RIE) system. An initial Ar/Cl 2 plasma etch was used to remove the Titanium adhesion layer, while O 2 plasma was used to etch the diamond. A typical ∼ 1.5 µm tall and ∼ 200 nm wide diamond pillar, imaged via electron microscopy, is shown in the inset of Supplementary Figure 1f. We then defined a set of Ti:Au (5:100 nm) coplanar waveguides via photolithography. The gaps between the central conductor and the ground plate of the waveguides were aligned (see Supplementary   Figure 1f) with the pillar rows via alignment markers defined on the diamond during the same O 2 etching described above.

Supplementary Note 2. Principles of stray field magnetometry
Because of the spatial confinement of its local spin density, to a volume below 1 nm 3 (see Ref. 2), a nitrogen-vacancy centre in diamond can be well approximated as a point-like sensor of magnetic fields.
Following other works 3,4 , we will consider the following Hamiltonian for this spin-1 defect: where D is the zero-field splitting, γ = 2.8025 MHz/G is the NV gyromagnetic ratio and , ⊥ indicate the directions parallel and perpendicular to the spin quantization axis of the colour centre.
In this work, our [100]-cut diamond hosts NV centres forming an angle θ NV = arccos(1/ In our experiments, both the upper (ω + ) and lower (ω − ) NV resonance frequency was measured at each point in space. The values for B and B ⊥ were then obtained using the following expressions 4 : In our experiment we measure the stray field B (x, y) in a plane at a distance d from the magnetic surface. We can call such quantity B (ρ, d), with ρ = (x, y).
At the probe position the stray field is curl free and it is therefore possible to define a magnetostatic potential φ M such that the vector field can be written as 5,6 : As pointed out in the past 5,7 , the previous relation implies that the stray field components are not independent. We will now derive the relation among the different stray field components and obtain further insight by starting from the following 2D Fourier transform of the quantity B(ρ, d), The vectors k and ρ are 2-dimensional vectors in reciprocal and real space, forming an angle φ k and φ with the x-axis (see also left of Supplementary Figure 2). As derived in Ref. 4, when the stray field is produced by a sheet of magnetic dipoles distributed over a thickness t and with , the stray field can be written in momentum space as 4 : where the expression for the traceless symmetric D(k, d) kernel matrix reads: Note that with respect to Ref. 4, we have included the finite size film thickness t by assuming the local magnetization vector m(ρ j ) to be constant through the magnetic film thickness and therefore integrating that dimension out. M s is the nominal, space-independent, saturation magnetization of the magnetic film.
From eq. (7) we can immediately realize that the rows of the matrix D(k, d) are not independent and for this reason D(k, d) is not invertible. In general, in eq. (6) it is impossible to obtain m by simply measuring all the components of the vector B.
In momentum space, the algorithm relating the stray field component along the z-axis B z (k, d) (see Supplementary Figure 2) to B (k, d), for an NV lying in the zx-plane, can be simply written as: In a similar way, all vector field components can be reconstructed without singularities from a single measurement of B (ρ, d) within a whole plane, provided θ N V = π/2.
Finally, we point out that the expression in eq. (7) includes also an analogy with the Huygens principle in optics. In particular, in order to reconstruct the field within a plane at a different dis- Once B z is known, all the other field components can be reconstructed according to (7). For Using full knowledge of all the stray vector field components, in Fig. 2 of the main text we have therefore reconstructed the expected magnitude of the stray field transverse to the NV axis and originating from the magnetic disc, using the expression: Note that the reconstructed map of B ⊥,r (ρ, d) will exactly match the map extracted from the spin level mixing given knowledge of the uniform bias field, which in our case is known up to the direction of a small perpendicular component.
As we shall see, eq. (7) allows for an intuitive real-space interpretation.
We define the real-space expression for the resolution function α x,y (d, t) of the in-plane magnetization as: Note that if t ≪ d the previous resolution function can be simplified as: In the same way, we define the real-space expression for the resolution function α z (d, t) of the out-of-plane magnetization as: Once more, if t ≪ d the previous resolution function can be simplified as: The spatial dependence of these resolution functions or effective point spread functions (PSFs) for magnetometry is plotted to the right of Supplementary Figure 2 for different sets of parameters.
With such notation, no approximations, and using the convolution theorem we can rewrite the real space expression for B(ρ, d) as: A single component of the stray magnetic field vector carries all the information, as all the other components are fixed given the first one. Due to the symmetry of our problem it is particularly illuminating to consider the B z component: with m x,y = (m x , m y ) the in-plane magnetization vector. It's clear that (15)  in far field such sheets compensate each other much faster than ∼ 1/r and more like an isolated dipole of the form ∼ 1/r 3 . It is evident from (11) and (13) that α z and α x,y indeed scale as ∼ 1/r and as ∼ 1/r 3 for very thin films.
Finally, in our experiments we have considered a stack of N = 10 magnetic thin films separated by a distance s. As the magnetization M s m is assumed to be constant through the film thickness, the stray field for the N = 1 case would read exactly like (15) with the difference that the PSFs are replaced by:

Supplementary Note 3. Magnetization reconstruction in the Bloch and Néel gauge
In Supplementary Note 2, eq. (15) provided us with a real space interpretation of stray field magnetometry. Since convolutions commute with derivatives, we can reformulate the problem of reconstructing the underlying magnetization pattern from the stray field measurements starting from Gauss's equation: where the two-component vector field F(ρ, d) plays the role of an effective electric field and the function B z (ρ, d) describes the charge density.
The effective electric field can be written down as: where m z (ρ) and m x,y (ρ) play the role of an effective scalar and vector potential, respectively. A solution to eq. (17) is defined up to a divergenceless term, which in our case can be written as: where C z (ρ, d) is an arbitrary function of space, a priori undetermined, and u z a unit vector perpendicular to the surface. We choose C z (ρ, d)u z to point in the z direction because F is oriented in the (x, y) plane. Our derivation diverges from classical electromagnetism (EM) 8 . In particular, in EM the curl of the vector potential is determined by a magnetic field measurement.
In our effective problem we only have access to B z meaning that ∇ × m x,y , and in turn C z , is fully undetermined.
Even if we had full knowledge of F, a second degree of arbitrariness in the knowledge of the vector and scalar potential comes, as in EM, from the following gauge-like degree of freedom: where Λ(ρ, d) is an arbitrary function of space.
As explained in the main text, in order to fix the arbitrary functions Λ(ρ, d) and C z (ρ, d) and therefore classify the different spin structures producing the measured stray field, we proceed in analogy with EM. Each physically distinct configuration of the spin texture is obtained after making local assumptions about the vector field m, with a procedure that resembles standard gauge fixing in EM 8 .
Two of these possible assumptions, motivated by the spiral (cycloid) nature of Bloch (Néel) domain walls 9 and the resulting partial differential equations that need to be solved in order to determine m are reported in the next Sections. By a solution in the Bloch gauge to the stray field equation, we mean a solution to (17) for m in which we make the local assumption: whose physical justification has been given in the main text.
Since m x,y plays the role of an effective vector potential, the condition (21) reminds us of the Coulomb gauge in EM. Exactly as in EM, in the Coulomb gauge the equation providing us with the scalar potential is the Poisson one: easy to solve in Fourier space. We also know that the solution to (22) is unique once boundary conditions are fixed 8 . Once m z (ρ) is found, we can then obtain m x,y by solving ∇ · m x,y = 0. The complete partial non-linear differential equation in the azimuthal angle φ(ρ) reads as: Eq. (23) takes normalization to the space-dependent saturation magnetization m s (ρ) into account (see also Supplementary Note 6 for more details on this last point). We obtain a solution to (23) variationally, by minimizing the following cost function with respect to φ: A very brief reminder of the popular steepest descent method we have used for minimizing the quadratic form in (24) is presented in Supplementary Note 4. By a solution in the Néel gauge to the stray field equation, we mean a solution to (17) for m in which we make the local assumption: Fixing the curl of m x,y is equivalent to fixing the curl for the effective electric field F or, equivalently, the function C z in (19). The vector field F becomes therefore conservative and it can be obtained explicitly from: At this point, an explicit solution to the stray field equation is still not possible as we retain the degree of freedom given by the arbitrary function Λ(ρ, d) (note that a transformation like the one in (20) preserves the curl of the vector field m x,y ). In order to further reduce the manifold of possible solutions, we introduce the normalization of the vector field m in the form of: with u φ the unit vector (cos(φ), sin(φ)). Eq. (27) represents two coupled non-linear partial differential equations in φ and m z . In order to produce Fig. 3e of the main text we have solved it by minimizing (with respect to φ and m z ) the following cost function, variationally: In Supplementary Note 5 we discuss the degeneracy of the solution once normalization and curl have been fixed.

Supplementary Note 4. Steepest descent minimization
The minimization of a quadratic form using numerical, iterative steepest descent procedures is reported in several textbooks 10 , for instance in the context of energy functionals.
In general, it is well known that a quadratic function C ({x α }) of N -variables x α , α = 1 . . . N can be minimized starting from the guess x α,0 , by iteratively moving antiparallel to the gradient direction, e.g. 10 : where λ is a constant and the i = 0, · · · N s -index refers to the iteration number. The previous follows from the fact that gradients are orthogonal to isolines directions: When C becomes a functional, N → ∞, {x α } → φ(α) and the functional increment upon a change φ(α) to φ(α) + η(α) can be written in first order as 11 : The analogy between (31) and (30) allows to rewrite the update in (29) in the continuous limit as 10 : where the derivative with respect to C is a functional one 11 . The functionals that we have to minimize in this work have, like (24) and (28), the general form: The functional derivative in (32) can therefore be computed using chain derivatives as: We then numerically implement (32) in order to obtain the function φ that minimizes C .

Supplementary Note 5. Uniqueness of the solution at fixed gauge
Fixing the curl of m x,y and locally imposing a normalization of the ordered moment does not guarantee that the solution reproducing a given target stray field will be unique.
In this Section we first discuss which kind of transformations would preserve the curl and normalization of the vector field and will finally briefly comment on the uniqueness of the solution. As known from standard EM, a gauge transformation like (20) would preserve the curl of the magnetic structure and the field it produces. On the other hand, in general that same transformation does not preserve the normalization of the vector field. We start rewriting (20) in the following form: If we want to locally have ||m 1 || = ||m 2 || we find that the following must hold: The previous equation fixes a condition for the norm of the vector ∆m = (∇Λ ′ , Λ). In particular, one can see that if we assume ||m 1 || = ||m 2 || = 1, then ||∆m|| = 2 cos(θ 1 − θ ∆m ), where θ 1 − θ ∆m is the angle between the m 1 and the ∆m vector. Eq. (36) can be easily solved for Λ in special cases, which allow us to prove that in general fixing the curl of m x,y and locally imposing a normalization of the ordered moment does not automatically guarantee a unique solution for the magnetic pattern.
Assume for instance that θ 1 = 0 everywhere in space, meaning we are considering a ferromagnetic pattern, which clearly has zero curl. Eq. (36) then reduces to: where u φ = (cos(φ ∆m ), sin(φ ∆m )) contains the azimuthal angle of the ∆m vector. Finding a solution to (37) is complicated by the presence of convolutions.
We make the only assumption that the sinusoidal functions in θ ∆m have a Fourier spectrum centred aroundk, such that the convolutions can be approximated with multiplications: We can now evaluate the ratio α x,y (k)/α z (k) using the results in Supplementary Note 2 and obtain: So, in essence, good solutions for ∆m are those in which its polar angle has a gradient with a constant norm. Functions linear with the spatial coordinate will be solutions, such as plane waves θ ∆m = (k · r)/2 with u φ ∆m k or radial waves θ ∆m =kr/2. All these solutions are therefore cycloids, in the definition given in the main text.
We found that cycloids with a single wave vector and constant ordered moment can actually produce no stray field, as ferromagnetic states also do not according to (15).
On the other hand, we also note that if we fix boundary conditions, e.g. assume the magnetic moments to be in a ferromagnetic state for r → ∞, then in order to avoid having B z = 0 at the boundary between the ferromagnetic and the plane wave region we need to selectk → 0.
We therefore argue that fixing the curl of m x,y , then locally imposing a normalization and boundary conditions selects a unique solution to the stray field equation.

Supplementary Note 6. Calibrations
In our work we have carried out a reconstruction of the underlying local magnetization configuration starting from eq. (15). The unknown parameters in the equation are the film thickness t, the NV depth d and the local value of the saturation magnetization M s · m s (ρ), which has been used for the reconstruction in eq. (23) and (27).
In order to calibrate these values we start from a simultaneous measurement of the magnetic disc topography and stray field map at saturation, as shown in Supplementary Figure 3. We first compare the stray field maps in Supplementary Figure 3a  If now we assume that a certain magnetic structure m produces a stray field B z (m), then a solution for θ(x, y) in (42) can be obtained by minimizing with a numerical variational analysis (see Supplementary Note 4) the following functional: The resulting function θ is found to be continuous and it's used to compute the solutions in Fig. 4 of the main text.

Supplementary Note 8. Micromagnetic simulations: stability of the Néel caps
To understand the origin of the observed right-handed chirality in a material with left-handed DMI, we performed full 3D micromagnetic simulations of a skyrmion in a multilayer system with ten magnetic layers. Here, we assume a CoFeB-based system with a saturation magnetization of c ν α x,y (d + ν · s, t).
Opposite chirality between, say, layer i and j can be imposed by simply setting c i = −c j . This results in a new α x,y (d, t). In this way, the magnetization pattern m(ρ) can still be considered as layer-independent but the in-plane magnetization will be summed up oppositely between the layers i and j; within the minimization process, this effectively inverts the relative chirality between the layers. In order to account for the intermediate layers hosting Bloch-like skyrmions, we have decided to run the minimization process still within the Néel gauge, but setting for the intermediate i-layer the coefficients c i = 0. This condition implies that the term ∇ · m x,y will not contribute, for those layers, to B z , as it should be for a real Bloch solution. The last method allows us to obtain a layer-independent m z profile, which we plot in Fig. 5c of the main text.