How to measure the local Dzyaloshinskii-Moriya Interaction in Skyrmion Thin-Film Multilayers

The current-driven motion of skyrmions in MnSi and FeGe thinned single crystals could be initiated at current densities of the order of 106 A/m2, five orders of magnitude smaller than for magnetic domain walls. The technologically crucial step of replicating these results in thin films has not been successful to-date, but the reasons are not clear. Elucidating them requires analyzing system characteristics at scales of few nm where the key Dzyaloshinskii-Moriya (DM) interactions vary, and doing so in near-application conditions, i.e. in systems at room temperature, capped with additional layers for oxidation protection. In this work’s magnetic force microscopy (MFM) studies of magnetron-sputtered Ir/Co/Pt-multilayers we show skyrmions that are smaller than previously observed, are not circularly symmetric, and are pinned to 50-nm wide areas where the DM interaction is higher than average. This finding matches our measurement of inhomogeneity of the magnetic moment areal density, which amounts to a standard deviation of the Co layer thickness of 0.3 monolayers in our 0.6 nm thick Co layers. This likely originates in small Co layer thickness variation and alloying. These film characteristics must be controlled with greater precision to preclude skyrmion pinning.

function from MFM measurements in saturating magnetic fields. The measurements show that the 0.6 nm thick Co layers vary on a 30-100 nm scale by up to ±1.2 ML (peak-to-peak), or equivalently vary with a 0.3 ML standard deviation per each Co layer in our multilayer structure. Our observation that the skyrmions are not circularly symmetric, are pinned, and have return point memory can likewise be traced back to this inhomogeneity. An important difference to other MFM work on skyrmions 23 is operation at room temperature on thin films, using calibrated tips which requires improved sensitivity and tip-sample distance control. It allows us to distinguish skyrmions from more extended skyrmionic bubbles and model them accurately.
The starting point for this work is the confirmation that our Ir/Co 0.6 nm /Pt-based multilayer (sample A; Fig. 1a) has a large net DM interaction, compared with a control sample B comprising repeats of Pt 1 nm /Co 0.6 nm (Fig. 1c), which on account of its symmetric interfaces is not expected to have a net DM interaction (though small values cannot be excluded 24 ). The as-grown and demagnetized states of sample A (Fig. 1b) consist of maze patterns of magnetic domains. They are narrower than the as grown pattern from sample B (Fig. 1d), which indicates a strong DM interaction.
The domain pattern contains information about the average DM interaction coefficient D, because it results from an equilibration of magnetostatic and wall energies. Through the latter the equilibrium domain configuration and the DM interaction are linked, a condition that can be exploited.

Average Dzyaloshinskii-Moriya interaction
For example, the width of regular stripe domains 25,26 that minimizes the total energy of the system (domain wallplus magnetostatic energy) has been used 19 to estimate D. Two conceptual difficulties with this approach make it inadequate in our case. The first is that the experimental domains patterns do not closely resemble regular stripes, and any statistic of the domain pattern chosen to represent the regular stripe domain width is therefore a priori unreliable. The second difficulty, of particular severity in multilayer films such as ours, concerns the account of magnetostatic energy and DM interaction when the system is not a single homogeneous layer, as closed analytical expressions assume 19,24 . Other work therefore calculates the energy of the multilayer relying on micromagnetic simulations of domains that resemble the measured ones 14 .
We use an alternative approach that relies solely on magnetometry data, the measured domain patterns obtained after demagnetization, and numerically calculated magnetostatic energies. Note that minimizing the total energy through variation of the width of domains in a regular stripe pattern 19 constitutes but one choice among many different variation paths in parameter space representing domain configurations. One could likewise select the spatial scale of a single arbitrary domain pattern for the variations. This approach is equivalent to varying the regular stripe domain width if and when the pattern of choice consists of regular stripes. More conveniently we can also choose a measured pattern, taking care it is close to the equilibrium. This avoids the problem of equating the experimental pattern to highly idealized stripes 19,24 , and of selecting which measure of the experimental pattern best represents the width of stripe domains. Once a domain pattern has been selected the energy is calculated numerically. For a given D, which determines the domain wall energy 27 , we find the measured pattern's scale for which the total energy is a minimum. In general, the minimizing scale will depend on the D. At one specific D value, however, the minimizing scale will coincide with the actual scale from the measurement. We take this D value to represent the average value of the DM interaction coefficient best. Thus for domain patterns obtained after demagnetization in an oscillatory field with decaying amplitude (applied perpendicular to the film plane) the average DM coefficient is D avg = 1.97 ± 0.01 mJ/m 2 ( Fig. 2a; and Supporting Material).
This value is higher than in structures with more repetitions of the Ir/Co/Pt unit 14 , which is likely related to the improved homogeneity attainable when fewer layers are used in the multilayer stack, as in our case. Conversely, our average D is significantly lower than the value we would have expected based on theoretical predictions by Yang et al. 21 , according to which |D Pt/Co + D Co/Ir | ≈ 5.65 mJ/m 2 for a Pt/Co(3 ML)/Ir system. The fact that we find D avg = 1.97 ± 0.01 mJ/m 2 could be an indication of interface roughness and intermixing (alloying) in our multilayer system. It is also noteworthy to point out that if we had modelled the system energy assuming regular stripe domains with a width given by the average domain size of 248 ± 26 nm obtained from the demagnetized domain pattern, we would have arrived at an effective D of 1.88 ± 0.03 mJ/m 2 (cf. Fig. 2b), a difference which our method can resolve clearly given the standard error attainable.
Based on the estimate of D avg the DM interaction is too weak to result in a negative domain wall energy. Hence we do not expect to observe arrays of skyrmions forming spontaneously, but isolated skyrmions remain a possibility 28 .

Room temperature skyrmions in thin film multilayers
The characteristics of the room temperature skyrmions of our Ir/Co/Pt-based multilayers become apparent in Fig. 3. After erasing the demagnetized domain pattern of Fig. 1b at −128 mT and then lowering the field amplitude to −1.1 mT we observe skyrmions (Fig. 3a) as small stable light spots with approximately circular geometry (up magnetization). As expected, no skyrmion-lattice but only a few isolated skyrmions exist in the 4.5 × 4.5 µm area of Fig. 3a. In increasing positive fields (+4 mT in Fig. 3b, and +8.5 mT in Fig. 3c,d), two of the skyrmions visible in Fig. 3a burst (strip out) into larger domains (examples indicated by the yellow arrows in Fig. 3b,c), and new spots of reversed magnetization with a non-circular geometry appear (examples given by red contours in Fig. 3d-f). These larger reversal domains are reminiscent of bubble domains characterized by a larger area of constant up magnetization at their centers.
Overall, these measurements reveal that skyrmions in our system are isolated, burst and become small non-circular domains in small positive fields (i.e. parallel to their core), and coexist with slightly larger reversal domains with non-circular geometry in stronger positive fields (see Fig. 3e,f). The data also show that skyrmionic bubbles (for instance areas in 3b to 3f enclosed by red ellipses) with larger diameters and non-circular boundaries can coexist with the smaller skyrmions. Additional data taken at a different location on the same film in small negative fields following saturation at −128 mT confirm that the skyrmions repeatedly appear at the same locations. Skyrmions are thus pinned at specific locations with physical properties different from those of most of the sample.
For detailed analysis we exclude those bright spots in Fig. 3a-f which on account of their distinctly irregular shape or larger size are likely beyond the bursting point, and thus inadequately analyzed in terms of skyrmion models. We acquire MFM images with 4.88 nm pixel width, shown in Fig. 3g-l. For example, Fig. 3g,j show two groups of two skyrmions each. The upper skyrmion in Fig. 3g has stronger MFM contrast and an almost circular geometry, whereas the lower one has a slightly elliptical appearance. Applying a weak field antiparallel to the core of the skyrmions erases them (the lower one at B = −2.9 mT, Fig. 3h, and the upper one at B = −3.5 mT, Fig. 3i), but both skyrmions re-appear if the field is set back to −1.1 mT. Two further skyrmions with a slightly elliptical shape are shown in Fig. 3j. One of those is erased by a small field of −2.1 mT, whereas the other one remains (a) For a measured (perpendicular) demagnetized state pattern near the energy minimum, the pattern scale is varied for a fixed parameter D (colored labels), and the total energy of the system is computed, accounting for discrete Co layers. The pattern scales resulting in an energy minimum for different D are connected with a broad red line as a guide. When the chosen parameter D coincides with the film's actual average DMI coefficient the broad red line crosses the scale 1.0, i.e. the minimum energy is found for the actually measured scale (modelled/ measured scale = 1.0). (b) Same as a except that a regular stripe pattern is assumed. Its width matches the 248 ± 26 nm average domain size measured from (a).
www.nature.com/scientificreports www.nature.com/scientificreports/ stable up to a field of −18.5 mT (Fig. 3k), becoming almost circular at this field. Figure 3l shows the MFM signal obtained at −80.4 mT after erasing the upper skyrmion of Fig. 3k. The fact that the skyrmions strip out at different fields is a further indication that these systems are not adequately described in terms of average quantities, for which that behavior ought to be the same.
The above observations imply that the values of the DM interaction coefficient D are dependent on position. Conversely, a detailed account of the magnetic structure's characteristics should provide, in principle, these position-dependent values. However, in practice, the insufficient resolution of the measurement method can be an impediment. That fact limited the information on local D that could be extracted from STXM data in prior work 14 . In this work we exploit the greater signal to noise ratio and spatial resolution of the MFM measurement to assess the position dependence of D. But we need to account for the fact that imaging at the limits of resolution, the details of the instrument's transfer function are important.

tip transfer function of the MFM/quantitative evaluation
To see this, consider that the image formation in an MFM system can be described with transfer functions in 2D reciprocal space. From a computational point of view, and in terms of the 2D reciprocal space vector k = (k x , k y ) and the spatial coordinate z, our instrument's transfer function ICF(k, z) enables us to relate the frequency shift map Δf(k, z) that our instrument reports to the sampled stray field (gradient) map H z (k, z): ICF(k, z) is a system-specific a priori unknown function, that is, it depends on cantilever stiffness, resonance frequency, oscillation amplitude and canting, tip-sample distance, and magnetic tip in a way that needs to be calibrated. The distance (and thickness-) loss factors 29,30 that describe the stray fields of sinusoidally varying magnetic charge distributions confer it a strong dependence on k (see Fig. 4a,b). Determining it for our system constitutes the calibration process 31 . It requires knowledge of Δf and corresponding H z for at least one case. We use sample B, which has perpendicular anisotropy and uniform magnetization through-thickness, and has adequately sized domains of 437 ± 64 nm in the in-plane demagnetized state. Note that a strong signal to noise ratio is needed in Δf(k, z) for an accurate calculation of ICF(k, z), which requires us to measure as close to the sample as possible, yet without altering the tip through contact. We accomplish this with a capacitive tip-sample distance control running on the second resonace of the cantilever 32 . In Fig. 4c we show a one dimensional representation of the result for our instrument's transfer function, more specifically the radially averaged amplitude.
Before making use of the calibration thus obtained to assess in detail the local skyrmion characteristics, we turn our attention to an important observation about the background in the images 14,19 , which our improved instrument sensitivity allows making.

Information contained in the magnetic background signal
The contrast of the 'background' , that is, the signal acquired far from the domain boundaries and skyrmions, is not uniform as one might expect at saturation. Here it amounts to an MFM contrast of ±0.5 Hz (see Fig. 5a showing MFM data acquired between the skyrmions shown in Fig. 3). These contrast variations occur with length scales of 100 nm and below. This scale is comparable to the skyrmions' width. From the relation between inhomogeneity www.nature.com/scientificreports www.nature.com/scientificreports/ in the domain walls' energy landscape and its effect on pinning 33 , we anticipate that if the background is magnetic it could have a large influence on the skyrmion stability. It is therefore important to establish the origin of this background inhomogeneity.
Because the background remains constant across all field levels we know that its structure is not caused by instrument noise. Considering potentially relevant interactions we conclude it could arise at shallow topographical features (roughness) by van der Waals interaction between tip and sample or result from an inhomogeneous sample magnetic moment areal density, the stray field of which would interact with the tip magnetization. Note that in the former case the interaction is independent of the magnetization, whereas in the latter it changes sign with the relative sign of tip and sample magnetization. This distinguishing feature allows establishing that the background is predominately of magnetic origin (see Supporting Material).
Recalling that the stray field of perfectly flat magnetic films of exactly uniform thickness and constant magnetization would vanish, we assess the possibility that the magnetic background signal is caused by departures from this condition.
First we consider roughness (Fig. 5b) in Co layers of uniform thickness. In our measurements the average tip-sample distance is actively controlled using slow feedback parameters such that the tip follows the average tilt of the sample, but does not compensate for local variations of the sample topography, which amount to ±0.48 nm RMS (not shown) and give rise to contrast. Using the calibrated ICF we can calculate the MFM contrast for a saturated film that has a periodical roughness of ±0.5 nm (Fig. 5d). We thus find that the contrast can be at most ±0.05 Hz, in effect ruling out the possibility of roughness as the cause for the contrast.
Next we consider the possibility a magnetic moment (or magnetization) inhomogeneity, which to a degree will arise naturally during the deposition process. It can be described with varying degrees of complexity, but for definiteness we assume here that it affects all Co-layers identically. The roughness of the magnetic surfaces that is necessarily implied can be disregarded on account of the negligible stray field associated with it, as we have argued above. Therefore, for the stray field calculation we take a flat layer of magnetization with inhomogeneous (in the plane) amplitude ΔM(k) corresponding to the local thickness-dependent moment of the Co layers. It gives rise to a contrast that can be calculated by convolution of the calibrated ICF(k, z) and the stray field of each layer of average thickness t Co at a reference distance, as per Eq. (1). In that equation, we can now express H z (k, z) as where z = 12 nm, t cap = 3 nm, t Pt = 1 nm and t Ir = 1 nm are the (average) thickness of the Pt and Ir layers. ΔM(k) can now be found from the measured background patterns (Fig. 5a, and also Fig. 3i,l) by deconvolution. For the background contrast (Fig. 5a) measured between the skyrmions shown in Fig. 5a (and Fig. 3i,l), we find a standard deviation of the magnetization of ±64 kA/m and a maximum magnetization amplitude of ±220 kA/m (Fig. 5e, and Fig. 3i,l), the latter amounting to a peak-to-peak of ±1.2 atomic monolayers, or equivalently a 0.3 atomic monolayer standard deviation per Co layer in our samples. Disregarding additional variations in the Pt and Ir layers, a peak-to-peak and RMS variation of the film thickness of 1.44 nm and 0.36 nm, respectively, can thus be expected. These values agree well with the film roughness values of 1.96 nm and 0.48 nm observed from topography measurements performed with a Bruker AFM in the peak-force mode. This suggests that the contrast from the ΔM(k) inhomogeneity can largely be attributed to thickness variations of the Co layers, but magnetic moment variations from alloying cannot be excluded. www.nature.com/scientificreports www.nature.com/scientificreports/ Note that according to the calculations of Yang et al. 21 we would expect to observe a local DMI above 5.65 mJ/m 2 for a perfect Pt/Co(2 ML) structure, and lower values for our multilayers with slightly diffuse or rough interfaces.

skyrmions and local values of DMI
Previous work was able to estimate a large net D in Ir/Co/Pt based multilayers by matching simulated magnetization profiles from skyrmions to the FWHM of the measured contrast spots 14 . The small skyrmions size, estimated to range between 30 and 90 nm, precluded a closer examination of variability in D with these measurements. Similarly, chiral domain walls with strong D were found in Pt/Co/MgO systems 18 , but the 130 nm wide chiral bubble domains stabilized through geometrical confinement could not be used to estimate local values of D.
To derive local values of D from the quantitative evaluation of high resolution MFM measurements of skyrmions, we model skyrmions following the original work from Bogdanov et al. 32 in which for a magnetization M saturating with M s , the reduced magnetization profile m = M/M s of the skyrmion is required to minimize the reduced energy functional Here h ext and h d are the external-and demagnetization fields 33 reduced with H D ≡ D 2 /AK (D is the DMI coefficient, A the exchange stiffness and K the anisotropy energy), β ≡  AK D / 2 is the reduced anisotropy constant, www.nature.com/scientificreports www.nature.com/scientificreports/ where M Co is the magnetization of a single Co layer containing a skyrmion formally taken to be uniform through thickness, H d is the local demagnetization or stray field, and M loc is the magnetization distribution of the stack of Co layers, i.e. it is equal to M Co inside the Co layers and is set to zero elsewhere in the film. Thus defined, α is a constant that allows approximating the magnetization-map dependent demagnetization field at each location with a quantity that depends only on the local value of the magnetization field. With the replacement h d = −α 2 m in Eq. (3) the corresponding Euler equation can readily be solved numerically (e.g. using Matlab). In our case we take the exchange stiffness to be A = 16 pJ/m, use the average uniaxial anisotropy K u = 414 kJ/m 3 and magnetization of the Co layer M Co = 653.6 kA/m (obtained from VSM). An initial value for α = + + ≈ . t t t t 5 /(5 4 4 ) 0 522 Co Co Pt Ir allows finding a first skyrmion profile, to be used in refining α according to Eq. 4 (see Fig. 6a).
Self-consistency is attained after one iteration. α varies slightly for different choices of D, as is apparent from the plot of Fig. 7a.
It could be argued that one should expect different values for A and K in this calculation, and interdependencies between A, D and K. The specific choice for the functional form of these relations is not necessarily clear, however. Initially, we use the same value of A = 16 pJ/m that we also used for the determination of the average DMI (see Supporting Material), and ask what would the coefficient D have to be in order to match the skyrmion profile with the average values of A and K, discussing further choices for these parameters subsequently. In this way we obtain the skyrmion magnetization orientation profiles e.g. for 3 mJ/m 2 < D < 3.5 mJ/m 2 , which we show in Fig. 7b. From them, we can directly produce a circularly symmetric magnetic charge pattern for each skyrmion (not shown) and then the corresponding z-component of the stray field (not shown) at the tip sample distance of 12 nm, used for the MFM measurements. This enables us to simulate the measured contrast by convolving the stray field map with the instrument calibration function, ICF(k, z), from where we can obtain the profiles in Fig. 7c.
As in most magnetic imaging methods, with the exception of sp-STM, in MFM the effective probe-sample interaction is a weighted average over a characteristic volume. Because of that, narrow magnetic objects can appear wider than they actually are. This means that small magnetic objects of arbitrary profile could result in contrast peaks resembling skyrmions. A particularly striking example of this effect is given by scaling the MFM profile of skyrmion with D = 3.0 mJ/m 2 to the same contrast as that with D = 3.4 mJ/m 2 (dashed blue line in Fig. 7c). Their apparent widths are roughly the same although the spin profile (Fig. 7b) differ by 30 nm. It becomes clear that calibrating the tip-and instrument-specific ICF(k, z) is necessary to remove this scaling degree of freedom. Figure 7d shows the peak contrast from the skyrmion simulations in Fig. 7c, and indicates that the maximum contrast is a proxy for D, provided the instrument calibration function is known. However, note that skyrmions with D = 1.97 mJ/m 2 would produce a contrast that is below noise level of our measurement (±0.02 Hz, cf. dashed horizontal line in Fig. 7d). www.nature.com/scientificreports www.nature.com/scientificreports/ The experimental data in Fig. 3j (See supporting material for a further example based on the data of Fig. 3g) contains a magnetic background, as argued, which we remove prior to comparison with simulations, i.e. we pointwise subtract the map of Fig. 3l (−80.4 mT) from that of Fig. 3j (−1.1 mT). This results in Fig. 6e out of which we extract the maximum contrast of the upper and lower skyrmion by averaging over several equivalent background subtracted images of the same area. The MFM contrast of the upper and lower skyrmion is 1.28 ± 0.04 and 1.10 ± 0.04 Hz respectively. With Fig. 7d D = 3.449 ± 0.020 mJ/m 2 and 3.426 ± 0.020 mJ/m 2 and for the upper and lower skyrmion respectively.
Note that the difference in skyrmion profile width is about 3 nm (Fig. 7f). The complete validation of the fit comes from the detailed comparison of the simulated skyrmion maps. These we plot in Fig. 7g. The sections across the skyrmions in images 7e and 7g can be seen in Fig. 7j-m. The black trace is the experimental result, highlighted in corresponding lines in Fig. 7e, whereas the superimposed red traces are the simulation result from Fig. 7g. Figure 7h shows the pointwise difference between experiment (Fig. 7e) and simulation (Fig. 7g). By construction, the simulation fits very well the center contrast of the skyrmions, and an approximately 30 × 30 nm 2 wide region surrounding it, but the quality of the fit degrades away from the center in some cases. Specifically, the deterioration is minimal for the cross-sections k and l but it is observable for the cross-sections j and to a lesser extent m (Fig. 7j-m). This indicates that the large values of D are circumscribed to small areas. Furthermore, matching the width of the skyrmion along the j-directions leads to a large overestimate of the center contrast (not shown). Consequently D must be changing along the j-direction (in particular) over length scales of the order of 30 nm which the isotropic model in Eq. (3) cannot capture in detail. Importantly, the large inhomogeneity of D, which www.nature.com/scientificreports www.nature.com/scientificreports/ we showed to be as much as 75% larger than the average, is a cause for pinning, of the kind we would expect to raise the current density threshold for motion, and consistent with skyrmions not moving when a field is applied. Figure 7i displays the calculated background magnetization variation (i.e. deconvolution of our ICF from the measurement in saturation), which is a proxy for the average layer thickness, and superimposes contours of the skyrmions. The figures show that the skyrmions are invariably localized in areas of relatively low relative magnetic signal (bright contrast), i.e. thin Co layers, which the zoomed-in insets show. This is consistent with the previously calculated strong film thickness dependence of D 21 . Regions in which all Co layers have the same constant small thickness are less likely to occur in our sputtered thin film, and therefore areas of high D are correspondingly sparse, consistent with our observations.
It is somewhat surprising that despite large significant differences between average and local values of D the analyzed skyrmions have a rather narrow distribution of D. The observation may be explained with the sparseness of large-D regions (and low thickness regions found by our analysis of the magnetic background contrast), and is connected with the observed absence of skyrmions with a weaker D than for the four assessed skyrmions. To understand why, consider that when we reduce the applied field from negative saturation to −1.1 mT the high D skyrmions, particularly for D > D κ=1 ≈ 3.28 mJ/m 2 (the skyrmion lattice threshold value) are the first to nucleate stably because their non-uniform magnetization lowers the energy (cf. the stability phase diagram in e.g. Kiselev et al. 27 ). As anticipated for our measured average D of 1.97 mJ/m 2 isolated skyrmions can also exist for D < D κ=1 , but their nucleation would require overcoming a barrier from the non-uniform magnetization in this case, and would be delayed to larger fields parallel to the skyrmion core (reversed or positive fields). Furthermore, such low-D skyrmions burst into large domains at relatively low reversed fields 27 . This is consistent with our experimental observation (c.f. Fig. 3). When lowering the field from negative saturation and reversing it to +1.1 mT, therefore, skyrmions at film positions with a high D nucleate first, and are then stabilized by surrounding areas of lower D, which hinder their expansion into bubble domains at sufficiently small reversed (positive) fields. Support for this hypothesis is provided by our simulations which reveal that skyrmions with such large values of D are rather unstable, i.e. a small further increase of D will lead to the expansion of the skyrmion into a (chiral) bubble domain.
As the reverse field is increased, skyrmions with lower D can nucleate, but may immediately expand into a (reversed) bubble domain (see e.g. red circle in Fig. 3b), or will be swallowed by an expanding reversed domains arising from bursted skyrmions with higher D.
An study of the skyrmion profile sensitivity with different choices of A, K u and D can be seen in Fig. 8.
In panel 8a we carry out a sensitivity analysis for the results of Fig. 7. At the center of the panel we show the simulated measurement result for a skyrmion calculated with K u and A kept as for the calculation of the average D. Along each of the three axes one parameter is varied independently (A: ±0.25 pJ/m, K u : ±5 kJ/m 3 and D: ±0.03 mJ/m 2 ) of the other ones, leading to skyrmion profile simulations that clearly differ from measurement. The discrepancy between the frequency shift simulations and the measured skyrmion is shown in panel 8b. It is clear that the solution found in Fig. 7 is stable for the independent parameter variations shown.
In a more general study of the stability of the local fit to the skyrmion we analyze the skyrmion profiles for A, K u and D on a grid in D − K u − A space near the solutions previously found in Fig. 7. Several parameter combinations lead to skyrmion MFM images that match the measure data reasonably well (RMS difference). They are indicated by black dots in Fig. 8c. An interpolating surface suggests there exists a continuous set of parameters leading to skyrmion profiles about matching the experiment. The local DMI could indeed be different from the values near 3.4 mJ/m 2 found for the skyrmions in Fig. 7. We also attempted an exploration of the parameter space into regions with smaller average D (but still significantly larger than the average D = 1.97 mJ/m 2 ) by considering much smaller A and K u . However, it was increasingly difficult to find stable skyrmion solutions and the required values for A and K u appeared rather unreasonable.
Departures from the average magnetization could also be investigated. The details of the corresponding simulations are provided in the supporting material. As one would expect, they confirm that variations of the magnetization (or areal moment density) affect the fitted value of the DMI, leading to a priori different values of the local D. We find that the quality of the 2D fit to the experimental skyrmion image improves slightly with decreasing M Co . However, below about 500 kA/m the matching of the measured MFM contrast would required even larger D-values for which the skyrmion solutions become unstable in zero-field. Between 653.6 kA/m and the stability limit, the value of D decreases slightly or increases rapidly in the simulations, depending on the path along which A and K u are varied together with D. We conclude that changes in the magnetization do not affect the final conclusion of a local value of D largely exceeding the average.

Conclusion
We observed skyrmions corresponding to DM-stabilized high-curvature core magnetization profiles in sputtered thin films multilayers. The DM interaction at the location of the skyrmions can significantly exceed the average value of 1.97 ± 0.02 mJ/m 2 reaching 3.45 ± 0.04 mJ/m 2 when the exchange stiffness and anisotropy retain average values. For smaller values of the DMI at the skyrmion location to be consistent with the data, the anisotropy or exchange stiffness would have to be significantly lower as well, in effect indicating that the DMI values significantly exceed the average. At these high levels of DMI the skyrmion lattice should be stable, suggesting that the skyrmion lattice could exist in sputtered thin films with additive DMI, which is critical for applications. Indeed, similar structures were recently shown 19 . But the inhomogeneity in the DMI (and A and K u , as discussed) is of the order of the skyrmion size, so that it will strongly affect the shape of the skyrmions and lead to strong pinning, as we found by MFM here. The underlying cause of the DMI inhomogeneity is likely the imperfect interlayer surface, characterized by possible alloying and Co-thickness inhomogeneity, revealed by the stray field at saturation. The magnetization variability in the Co layers is as large as 260 ± 64 kA/m, which would amount to a correlated thickness variation in the Co layers of up to ±1.2 atomic monolayers. The skyrmions are not circularly symmetric; (2019) 9:3114 | https://doi.org/10.1038/s41598-019-39501-x www.nature.com/scientificreports www.nature.com/scientificreports/ more extensive theoretical studies of skyrmion profiles in inhomogeneous D will benefit from comparison with these data. Overall, MFM offers a convenient way for studying the average and local D, and can even provide a measure of the subtle local magnetic moment variations in industrially relevant multilayer systems.

Methods
Film structures and preparation. We prepared two thin film magnetic multilayer structures by DC magnetron sputter deposition at room temperature, using Ar gas at 1.8 × 10 −3 mbar. The system's base pressure is 2 × 10 −8 mbar. Sample A is Si (nat.) \Pt 10 nm \Co 0.6 nm \Pt 1 nm [Ir 1 nm \Co 0.6 nm \Pt 1 nm ] x5 \Pt 3 nm where Si (nat.) stands for the Si (100) substrate with native oxide. Sample B is a control sample sputtered analogously to sample A but with symmetric interfaces, i.e. Si (nat.) \Pt 10 nm \[Co 0.6 nm \Pt 1 nm ] x5 \Pt 3 nm . compensated by the Co\Pt interface's at the other side. This symmetry is broken in sample A and a net DMI can arise. A first manifestation of this fact is the noticeably different size (>1000 nm in sample A and (246 ± 40) nm in sample S) of the magnetic domains in the as-grown state observed by MFM in either case, shown in Fig. 1b,h. Magnetic characterization. For the macroscopic magnetic characterization we used a vibrating sample magnetometer (Quantum Design PPMS). We obtain an anisotropy K u = 414 kJ/m 3 , magnetization of the Co layers of 653.6 kA/m, resulting in an effective anisotropy of K eff = 146 kJ/m 3 . For the microscopic characterization we used a room temperature magnetic force microscope operated in vacuum (10 −6 mbar) for better sensitivity