The evolution of skyrmions in Ir/Fe/Co/Pt multilayers and their topological Hall signature

The topological Hall effect (THE) is the Hall response to an emergent magnetic field, a manifestation of the skyrmion Berry-phase. As the magnitude of THE in magnetic multilayers is an open question, it is imperative to develop comprehensive understanding of skyrmions and other chiral textures, and their electrical fingerprint. Here, using Hall-transport and magnetic-imaging in a technologically viable multilayer film, we show that topological-Hall resistivity scales with the isolated-skyrmion density over a wide range of temperature and magnetic-field, confirming the impact of the skyrmion Berry-phase on electronic transport. While we establish qualitative agreement between the topological-Hall resistivity and the topological-charge density, our quantitative analysis shows much larger topological-Hall resistivity than the prevailing theory predicts for the observed skyrmion density. Our results are fundamental for the skyrmion-THE in multilayers, where interfacial interactions, multiband transport and non-adiabatic effects play an important role, and for skyrmion applications relying on THE.


Supplementary Note 1. EXPERIMENTAL METHODS
Multilayer films consisting of Ta(3)/Pt(10)/[Ir(1)/Fe(0.5)/Co(0.5)/Pt(1)] 20 / Pt(2) (layer thickness in nm, in parentheses) were deposited on thermally oxidized 100 mm Si wafers by DC magnetron sputtering at RT, using a Chiron TM UHV system manufactured by Bestec GmbH. The base pressure prior to film deposition was 1 × 10 −8 Torr, and a working pressure of 1.5 × 10 −3 Torr was maintained during the deposition by Ar sputtering gas. Buffer layers of Ta(3)/Pt (10) were deposited before the active stack for better adhesion and film texture. Additional capping by Pt (2) was added to protect the stacks from unwanted oxidation. Film texture, surface roughness and interface quality were verified by X-ray diffraction (XRD), atomic force microscopy (AFM), and X-ray reflectivity (XRR) [1]. Magnetization The magnetization [M (H)] measurements were carried out using a superconducting quantum interference device (SQUID) with a magnetic property measurement system (Quantum Design MPMS XL). Magnetotransport measurements were carried with a lock-in technique using small current densities (as low as 10 5 A/m 2 ) to avoid current induced changes in the magnetic textures. The data was acquired at various temperatures (T ) between 5 K to 300 K through a full hysteresis cycle, with 5 mT steps between the negative saturation field (H − S ) and the positive saturation field (H + S ) after saturation at large fields (+4 T). Measurements were carried out using an physical property measurement system (Quantum design, PPMS 6000), complemented by custom built variable temperature insert (VTI) housed in a high field magnet [1].
The magnetic imaging experiments were conducted using a cryogenic frequency modulated MFM system [2]. The sample was first stabilized at a given temperature (5 − 200 K) and then magnetized in the out of plane direction by applying µ 0 H = −0.5 T, well beyond the saturation field (H OP s ). After saturation, MFM images were acquired at field increments as H was swept from H − S to H + S . Images were acquired by scanning the magnetic MFM tip parallel to the film surface at a constant scan height h without feedback, while recording changes of the resonant frequency (∆f ) of the cantilever from its natural value (f 0 ) due to magnetic forces acting on the tip [3]. These shifts provide information on magnetic texture through the magnetic forces acting on the tip. In the limit of small oscillation amplitude and small frequency shift [2]: ∆f ≈ f 1 −(f 0 /2k 0 )(∂F z /∂h), where k 0 is the spring constant of the cantilever [4], and f 1 is a constant offset. It is possible to account for the finite oscillation amplitude [5], but as we have demonstrated in recent work [6], this is not needed to describe the signal phenomenologically. It is important to note that the signal magnitude strongly depends on the tip-sample separation height h. Therefore, comparing the variation in ∆f between scans is only meaningful if they were obtained at very similar values of h.  [7]. The details of magnetic parameters such as D, A, K eff and skyrmion properties for various thicknesses of Fe/Co have been reported previously [1].

Supplementary Note 3. MAGNETIC FIELD CALIBRATION
The estimated topological Hall effect (THE) in these multilayers is 5% − 10% of the total Hall signal, and any field offsets between the experimental setups can often appear as residual signals. In recent work, we have established the robustness of the residual signal at RT, and the results were found to be consistent across different measurement setups [1].  Temperature dependence of the effective out-of-plane anisotropy and the saturation magnetization.
Error bars in K eff represents the systematic error in determining the magnetic saturation field and the magnetization.
relevant to the samples studied here. To avoid any field offset discrepancies resulting from the variation in the field protocols, the calibration runs were performed using the same protocols used for measuring the actual samples.

A. Topological Hall Analysis
In case of B20 skyrmion lattices, such as in bulk crystals of MnSi and FeGe, the measured raw Hall resistivity (ρ yx ) shows a characteristic hump [8,9] which correlates with the emergence of a skyrmion lattice. This characteristic feature in B20 clearly differentiates the transport data from magnetization. However, as mentioned in the main text, this feature does not always appear. In our multilayer case the characteristic features in the raw Hall data are subtle and need more careful inspection. As described in the main text, the topological Hall signal is extracted as the residual of the fit of ρ yx to ρ f it yx = R 0 H + R S M (H). It is thus important to establish the validity of this residual signal and the analysis employed. Here, by comparing the first derivatives of the measured data −dρ norm yx /dH and dM norm /dH, we show attributed to the chiral domain walls stabilized by DMI [10]. However, the impact of such chiral domain wall resistance on ρ yx is not clearly understood yet. Supplementary Fig. 4 clearly shows that ρ xx (H) has no effect on our extraction of THE.
The anomalous Hall effect in B20 systems, like FeGe, is well described by linear [ρ xx (H), skew] and quadratic [ρ 2 xx (H), side jump and intrinsic] terms [11,12], with transport usually dominated by the side jump mechanism for pure FeGe, and skew scattering when doped [11]. However, the scaling of ρ yx with ρ xx is not well understood in multilayers [13] due to the inhomogeneous structural and magnetic environment experienced by charge carriers.
The scaling behavior in multilayers is expected to be controlled by the relative thickness of the layers and the mean free path of the charge carriers (l). When the layer thickness is less than l, transport is dominated by skew scattering [13]. This is consistent with the description of AHE by R S M (H) , Eq. 3 of the main text, generally employed for magnetic multilayers [14,15].
Supplementary Figs. 5 and 6 show the variation of ρ yx and ρ xx with temperature, field and the scaling behavior for different field values corresponding to skyrmions, chiral domains (at H < H S ) and a fully polarized ferromagnetic state (H > H S ). The figures show that within the field and temperature regimes explored in this study, the scaling between ρ yx and ρ xx can be approximated by a linear dependence. This is consistent with the dominant skew scattering expected in these multilayer systems [13]. A quadratic dependence is followed for 1 9 . 6 6 2 1 9 . 6 6 5 4 5 0 ( 1 )    The asymmetry in peak magnitude (H < 0 vs H > 0) results from hysteretic domain formation as the magnetization changes from negative saturation to positive upon field reversal [1].
As indicated in Supp. Fig. 7(b), the high field slope of the transport data representing the ordinary Hall coefficient, R 0 , changes systematically as the temperature changes from 5 K to 300 K. Supplementary Fig. 7(d) summarizes the fit parameters R s (T ) and R 0 (T ), which is also reproduced in Fig. 2(b) in the main text together with max (|∆ρ yx |). Importantly, R 0 (T ) changes sign near T = 85 K, while R S (T ) and max (|∆ρ yx |) show that they scale differently with T . This suggests that the anomalous and topological components originate from different Berry phase mechanisms and that ∆ρ yx is not a spurious signal associated with magnetization or the AHE component.  because in both cases it is parallel to the external field directly under the tip, where it is strongest. This indicates that the perturbation from the tip drives the magnetic texture to a lower energy state, implying the metastability of these textures.

A. Simulated MFM Images
The MFM images presented in the main text and in Supp. Figs. 10, 11 were obtained by convolving the transfer function of the tip with the calculated stray field from the simulation magnetization [19]. The stray field was calculated using [20]: where H z ( k) is the Fourier transform along the x − y directions, with k = (k x , k y ), k ≡ k 2 x + k 2 y , h is the distance between the MFM tip and the surface, d i is the distance between the i-th layer and the surface and m i ( k) is the Fourier transform of the simulated magnetization in the i-th layer. To account for the blurring effects of the tip and amplitude we have used the model and tip parameters from [6]. Since the MFM does not measure in 4 × 4 nm 2 resolution, we have down-sampled the images so that each pixel covers a 24 × 24 nm 2 area, which reflects the sampling intervals of the experimental measurements.  Therefore, for each temperature we took care to obtain an MFM scan at those peaks.

B. Topological Charge of Simulated Worms
Supplementary Fig. 12 shows the magnetic texture measured at the negative peak and Supp. MFM tip, the tip experiences a repulsive force giving rise to a strong positive signal. Such a strong signal from worms/skyrmions makes the identification easier. After labeling, the second step involves the fitting of each worm with a train of skyrmions. For each connected group of pixels labeled as "worms", we perform several iterations and attempt to find the optimal number of skyrmions which best describes the data.

A. Labeling
We start by applying the K-means clustering algorithm [21] on the ∆f values of a scan in order to divide the data into two groups -worms (including skyrmions) and background.
We denote the mean value for worms (the stronger signal) as c 1 and the background as c 2 (weaker signal). Each pixel in the image is thus attributed a value c α with α = 1 or α = 2, according to the group it belongs to.
Due to the effects of tip convolution and scan height, the magnetic features measured by the MFM are larger than their true size. Therefore, the smallest domain sizes measured by MFM are limited by the tip curvature radius of ≈ 30 nm, which is approximately the size of a pixel in our images 25 × 25 nm 2 . Also, based on our previous work [1,6], skyrmions with a radius of 40 nm are a reasonable lower limit. We feed this spatial information into our labeling analysis to avoid features smaller than two pixels, which are not physical.
The optimal labeling, given these two conditions, minimizes the following cost function: to find the optimal matrix of labels f . Here i, j are pixels in the image, i is a sum over all pixels labeled c 1 , i is a sum over all pixels labeled c 2 , D α i = (∆f i − c α ) 2 is a measure for how well the label α fits the ith pixel, V ij = λ if f i = f j and 0 otherwise, and i,j is a sum over nearest neighbors. The minimization of J(f ) is performed using graph cuts with a piecewise constant prior [22]. The value of the parameter λ is chosen as the smallest parameter for which the optimal labeling no longer includes single-pixel domains. The resulting labeling f is a binary image which is true for pixels positively identified as worms, and false for those in the background.

B. Fitting
For each connected group of pixels labeled as worms, we now fit the number of skyrmions.
We perform several iterations, assuming a different minimal spacing between skyrmions in each one. Similarly, K-means was employed to iteratively divide the worm pixels into groups, where the mean of each group is assumed to be the location of a skyrmion. If the mean of the group lies outside the worm area (as would happen with a small number of groups or banana shaped worms), we choose the closest point within the worm to be the new mean for the group. The worm is divided into an increasing number of groups until at some number N , the spacing between groups is smaller than the minimal spacing. We then fit N − 1 skyrmions, with the group means used as initial skyrmion positions. The fitting is done using a multipole expansion for the magnetic field from skyrmions (MEFS) model [6], in which we approximate the stray field from a skyrmion as the field of a magnetic dipole and quadrupole. Since we are only interested in the number of skyrmions which fits the data, we employ the truncated cone model [6] for the tip because it is computationally cheap.
Skyrmion-skyrmion repulsion grows exponentially on the scale of the skyrmion radius [23], therefore we reject fit results in which the distance between neighboring skyrmions is smaller than a skyrmion radius, which we previously estimated to be approximently 40 nm [6].
In order to provide some constraint on the amplitude, we first consider the case of a single skyrmion; as we have stated before, a single skyrmion will be captured by several pixels. Given the sharp shape of the skyrmion signal, the peak will be larger than the mean value of all the pixels which constitute the skyrmion. Extending this condition to worms, we constrain the amplitude of each fitted skyrmion to be at least the mean value of all pixels labeled as worms.
The fitting procedure that we establish therefore assigns the minimum number of skyrmions which best fits the data by constraining the amplitude per skyrmion. As a lower bound we take the numbers of skyrmions which can fit in a worm, and are similar in width to skyrmions observed in Fig. 1(j) (≈ 100 nm) since if their densities were lower we should able to resolve them. For an upper bound we take the number of skyrmions which can fit into a worms so that their centers are still 40 nm apart which is the typical radius [6]. The lower bound thus corresponds to the width of skyrmions which we have measured, and the upper bound corresponds to our estimate of the actual skyrmion radius without tip convolution [6].  Fig. 7(c). We would like to obtain some insight on whether or not this contribution to ∆ρ yx comes from worm-like features, as was clearly the case for images corresponding to fields between the peaks in ∆ρ yx and saturation. Our hypothesis is that as the dense stripe domain pattern evolves upon reducing the field towards µ 0 H = 0 T, some of the worm-like features observed at higher fields persist and coexist with stripe domains and contribute to THE. In order to extract the density of worms and hence topological charge from these images, we need a vigorous method for classifying pixels to worms and dense domain background.
Here, classification which relies on the ∆f value is problematic, as both up-domains and worms exhibit a similar positive signal. Establishing such a classification based on a quantitative metric which will be reliable for a range of scans obtained over a wide variety of temperature, field and scans height using a set of predetermined features can be remarkably difficult, as such general features which distinguish a worm from a dense domain pattern are very hard to define. Here, we pursue a more data-oriented approach in which an algorithm extracts relevant features from supplied examples. We assume that a feature contributes to THE signal if it resembles a worm [such as in Supp. Fig. 16(d)] more than it resembles dense domains [such as in Supp. Fig. 16(g)] and hence the parameters corresponding to worms in Supp. Fig. 12(a) act as reference for classification of features in dense images.
Deep learning models, specifically convolutional neural networks, are especially well suited for this task [24], as they can learn to distinguish complex features accurately using a large number of successive convolution layers. In addition to classification, we require the model to output spatial information -specifically which pixels in the image correspond to worms and otherwise dense domains. As shown by Long et al. [25], it is possible to add several upsampling layers to a classification convolutional network in a directed acyclic graph (DAG) architecture, and thus obtain spatial information in addition to classification. In the data analysis for this paper, we have used the MatConvNet toolbox [26], which provides an implementation for DAG neural networks that can be used in Matlab. The classification network we used was the pretrained VGG-VD16 network, downloaded from the MatConvNet website. We adapted the algorithm for our purpose by changing the last fully connected layer to have 2 classes only, and adding upsampling layers with 2 classes for the FCN8 architecture [25].
In order to minimize the training time, we try to make maximal use of existing knowledge, and only fine-tune the network to our needs. The VGG-VD16 network was trained on RGB images. Since the measurement noise in our MFM is at least 10 mHz, and the maximum values we measure are smaller than 100 Hz, we require only four significant digits to accurately describe our data. We can therefore convert our MFM data (single channel in double-precision floating-point format) to RGB images (three channels of unsigned 8-bit) without any loss of accuracy. We achieve this by creating a histogram of the data with the number of bins adjusted so that the bin centers are 10 mHz apart (noise limit). Each in this paper. Each image in the intermediate field range was then fed through the network. The output is a binary image which is true for worms and false for the background. An example for the network application on data can be seen in Supp. Fig. 15. We then copied the pixels classified as worms (with a margin of two pixels) from the original image and placed them in a background of zeros. A Gaussian filter of 5 × 5 pixels with σ = 10 was then used to smoothen the edges. The result is an estimate of a scan which has only the worms in it Given the true negative and positive rates of 0.98 and 0.76 respectively, we proceeded to estimate the classification errors. Our estimate is based on the probability that the network will make a classification error for a single pixel. For simplicity, we will assume that the prediction of a pixel class is independent of other pixels, and that it follows a binomial distribution. For convolutional networks that assumption is not necessarily accurate, but it is sufficient for determining error bounds.
We can calculate an explicit upper bound for how many worms the network missedlet us assume a worst case scenario where the network has identified zero worms, therefore there are N = 256 2 pixels labeled as dense domains. The probability that the network will miss a worm is p = 0.02. The probability that the network missed more than 3000 pixels who should have been labeled as worms is thus negligible. From Supp. Fig. 17(a), ≈ 3000 pixels were labeled as worms, and this was the lowest density we measured. We therefore conclude that the density which can be attributed to worms the network missed is negligible for our data.   Fig. 9 where taken from.