Do Images of Biskyrmions Show Type‐II Bubbles?

Abstract The intense research effort investigating magnetic skyrmions and their applications for spintronics has yielded reports of more exotic objects including the biskyrmion, which consists of a bound pair of counter‐rotating vortices of magnetization. Biskyrmions have been identified only from transmission electron microscopy images and have not been observed by other techniques, nor seen in simulations carried out under realistic conditions. Here, quantitative Lorentz transmission electron microscopy, X‐ray holography, and micromagnetic simulations are combined to search for biskyrmions in MnNiGa, a material in which they have been reported. Only type‐I and type‐II magnetic bubbles are found and images purported to show biskyrmions can be explained as type‐II bubbles viewed at an angle to their axes. It is not the magnetization but the magnetic flux density resulting from this object that forms the counter‐rotating vortices.

Polycrystalline samples of MnNiGa with a target composition of (Mn 0.5 Ni 0.5 ) 0.65 Ga 0.35 were prepared by arcmelting pellets of stoichiometric amounts of manganese (99.99%, Sigma-Adrich) and nickel (99.9%, Alfa Aesar) powders and gallium shots (99.99%, Alfa Aesar) in an arc furnace under a high-purity (5N) argon atmosphere. A 2 mol% excess of Mn was added to compensate for losses during the arc-melting procedure. The sample buttons were melted and flipped several times to ensure homogeneity. As-prepared buttons were annealed in vacuum at 800 • C for 1 week to stabilise the structure by removing strains and potential disorder and to enlarge the grains before being quenched in water to retain the high temperature structure.
The phase-purity was determined by x-ray diffraction using a Panalytical X-Pert Pro powder x-ray diffractometer operating in Bragg-Brentano geometry. This showed that the crystal was hexagonal P 6 3 /mmc (space group 194) with lattice parameters a = b = 4.15125(3)Å and c = 5.32774 (8)Å which is in good agreement with the structure reported in refs. 1 and 2. According to both references, Ga atoms occupy the 2c sites with Mn preferentially occupying the 2a and Ni preferentially occupying the 2d sites (see Figure S1).
The composition was determined to be Mn 0.325(1) Ni 0.324(1) Ga 0.350(1) by energy dispersive x-ray spectroscopy conducted using a Zeiss SUPRA 55-VP scanning electron microscope equipped with large area SDD EDX detector. (The number in brackets indicates the error in the last decimal place.) A Quantum Design Magnetic Property Measurement System, MPMS-5S, superconducting quantum interference device (SQUID) magnetometer was used to investigate the magnetic properties of the polycrystalline samples. The Curie temperature was T C = 340±5 K and the * j.c.loudon@gmail.com saturation magnetization at 5 K under an applied field of 5 T was 1.66 ± 0.02 µ B per formula unit, both in good agreement with previously reported values [2].

S2. SAMPLE PREPARATION
MnNiGa samples were prepared for electron microscopy and x-ray holography using an FEI NanoLab-600 Helios Dual-Beam focused ion beam (FIB) microscope equipped with an Omniprobe-200 micromanipulator. Samples of approximate size 10 × 5 µm were cut from a single grain using gallium ion milling, lifted out with the micromanipulator and mounted on Omniprobe grids. As the orientation of the grain was unknown, two initial samples were cut parallel to the grain boundaries of the material and their orientation found ex-situ using the transmission electron microscope. The first sample was close to the [100] zone axis and showed no magnetic contrast, the second was close to the [011] and showed only striped magnetic contrast. From these, the tilting angles required to cut a sample in the [001] zone axis were calculated and the final samples were cut within 7 • of the [001] zone axis. Two samples were cut from the same grain and thinned to about 200 nm. One sample FIG. S1. Crystal structure of Mn 1/3 Ni 1/3 Ga 1/3 . was used for x-ray holography and the other for transmission electron microscopy. Both striped magnetic domains and bubbles were observed in this orientation in agreement with the findings of Wang et al. [2].

S3. X-RAY HOLOGRAPHY
MnNiGa samples were prepared as described above and mounted on 3 µm diameter circular apertures patterned on Si 3 Ni 4 membranes which had been coated with a 600 nm thick gold mask. Reference slits were cut alongside these holes, 6000 × 10 nm in size. The membranes were mounted in a vacuum chamber, containing an Halbach array, capable of applying a magnetic field of up to 0.5 T in any direction. A thermal resistor in contact with the sample was used to control the temperature.
X-ray holography was undertaken on the COMET instrument at the SEXTANTS beamline at SOLEIL, using holography with extended reference by autocorrelation linear differential operator (HERALDO) [3]. In this technique, a reference wave is produced using an extended slit instead of the more common pinhole aperture, leading to an improved signal-to-noise ratio and spatial resolution [4]. Magnetic contrast was generated by exploiting x-ray magnetic circular dichroism (XMCD) at the L 3 absorption edge of Mn at 637.5 eV. Far field diffraction patterns formed via the interference of the sample and reference aperture were recorded using a 2048 × 2048 pixel charge-coupled device (CCD). By taking the difference between patterns recorded with opposite x-ray helicity then applying a linear differential filter and a Fourier transform, reconstructed images of the local magnetization in the samples were obtained. The x-ray polarization was controlled with an HU44 APPLE II undulator [5].

S4. ELECTRON MICROSCOPY
Electron microscopy was undertaken using an FEI Titan 3 80-300 transmission electron microscope operated at an acceleration voltage of 300 kV and equipped with a high-brightness XFEG electron gun which increased the electron count in our images by a factor of about ten times compared with microscopes equipped with conventional field-emission guns. This allowed us to acquire meaningful images at low defocus despite the low contrast.
In normal operation, the electromagnetic objective lens applies a 2 T field to the specimen which would force it into the saturated state. Images were instead acquired in low-magnification mode where the image is formed using the diffraction lens and the objective lens was weakly excited to apply a small magnetic field parallel to the electron beam. The applied field corresponding a given objective lens current had been calibrated to within 1 mT using an FEI Hall probe holder.
Images were filtered using a Gatan 865 Tridiem energyfilter so that only electrons which had lost between 0 and 10 eV on passing through the specimen contributed to the image and an aperture subtending a half-angle of 0.14 mrad was used to ensure that only the 000 beam and those associated with magnetic scattering contributed to the image. Images were recorded on a 2048 × 2048 pixel charge-coupled device (CCD).
The defocus and magnification were calibrated by acquiring images with the same lens settings from Agar Scientific's S106 calibration specimen which consists of lines spaced by 463 nm ruled on an amorphous film. The defocus was found by taking digital Fourier transforms of these images and measuring the radii of the dark rings that result from the contrast transfer function of the lens using the method described in ref. 6 (ch. 28).
The specimen thickness was measured by taking infocus images of the same area with and without energy filtering. Dividing the unfiltered image by the filtered image and taking the natural logarithm gives the thickness as a multiple of the inelastic mean free path of the electrons -a 't-over-lambda' map [6] (ch. 39.5). The absolute specimen thickness was measured at 15 positions using the two-beam convergent-beam electron diffraction technique described in ref. [6] (ch. 21.2) and the inelastic mean free path calibrated from this, giving a map of the absolute thickness across the whole specimen to within 10 nm accuracy.

S5. MICROMAGNETIC SIMULATIONS
We performed micromagnetic simulations using the Object Oriented Micro-Magnetic Framework (OOMMF) code [7,8]. From the domain wall width of δ ≈ 47 nm we estimate a ratio between the exchange constant and uniaxial anisotropy of the system as A/K u = 2.3×10 −16 m 2 . In a recent study, Zuo et al. [9] reported a magnitude for the saturation magnetization of µ 0 M s = 0.547 T. In order to specify a exchange constant magnitude A, we did a parameter sweep starting with a hexagonal array of bubbles in 1 µm × 1 µm × 200 nm samples and varying M s , A and the applied field. From the results, we observed the stabilization of bubbles around fields of 0.2 T for µ 0 M s = 0.550 T and for A between 10 and 20 pJm −1 . By decreasing M s , bubbles are observed when decreasing the exchange constant but at very weak fields. In addition, setting µ 0 M s ≤ 0.550 T, it was not possible to stabilize an isolated bubble at fields close to 0.2 T. From the simulations, a larger M s magnitude allows the stabilization of both bubble arrays and isolated bubbles around 0.2 T. Thus, we chose an exchange field of of 20 pJm −1 , and a saturation magnetization of µ 0 M s = 0.648 T. Any magnitude of A and M s around these values should also be valid for the stabilization of the bubbles. An optimal choice of micromagnetic parameters is beyond the scope of this work, since the observed experimental structures are in agreement with our simulations. With our choice of parameters, we estimate a quality factor Q of about 0.52.
Using the magnetic parameters discussed previously, we analyzed the possibility of generating bi-skyrmion configurations through simulated hysteresis-like processes. We defined cuboid samples of 1 µm × 1 µm × 100 nm size considering exchange, dipolar and uniaxial anisotropy interactions. We discretized the system using cells of 4 nm × 4 nm × 4 nm volume, whose dimensions are well below the exchange length constant of the material l ex = 2A(µ 0 M 2 s ) −1 of 10.94 nm magnitude. The systems are initialized at zero field with randomly oriented spins. After energy minimization we found spiral configurations with varying chirality, as shown in Figure S2 and Video S1.
Across the thickness of the sample the helicities of the domain walls change from a Bloch-like profile at the center of the sample to a Néel-like profile towards the surfaces like that shown in Figure 4(a) (main text). Interestingly, we notice the presence of adjacent domain walls with the same chirality. Increasing the field up to 0.2 T, we observe the formation of the elongated spiral domains. Among them we notice elongated spiral domains surrounded by domain walls of equal chirality. This occurs by the presence of Bloch lines at the extremes of the longest domain dimension. Furthermore, we observe the emergence of different types of bubbles which coexist with the spiral domains: ordinary bubbles, which are of Néel type at the top and bottom surfaces of the sample and Bloch type at the sample center (across the thickness). We also find type-II bubbles of ellipsoidal shape, which are characterized by two Bloch lines at the extremes of the bubble long axis, that propagate through the sample thickness. Above fields of 0.2 T, only isolated bubbles remain in the system. Type-II bubbles are observed for fields below 0.26 T. From these simulations we did not directly find a magnetic ordering resembling a bi-skyrmion, however we noticed that the type-II bubbles have a profile that, averaged through the sample thickness, can generate a bi-skyrmionic state. Hence we focused our attention to isolated bubbles.
Simulations of isolated type-II bubbles were performed in samples of L µm × L µm × 200 nm volume, with L varying from 600 nm up to 1.2 µm. As a initial state we specified a Bloch-like bubble propagating through the thickness of the sample, with the upper and lower half of the bubble (along the y-direction) with opposite chiralities. After energy minimization, at applied fields below 0.9 T we find elongated domains surrounded by domain walls of the same chirality. At the surfaces of the sample the helicity of the domain wall defining the stripe domains vary as in the hysteresis-like simulations. At sufficiently strong applied fields we find isolated type-II bubbles of ellipsoidal shape with the same varying helicity profile of the stripe domains across the thickness. With increasing fields the bubble profile becomes more spherical and decreases in size, as shown in Figure S3. This figure also indicates that the stability of these bub- S2. Micromagnetic simulation of a field sweep process on a MnNiGa sample of 1 µm size and 100 nm thickness. The initial state was generated from randomly oriented spins. The hysteresis-like process was simulated using steps of 20 mT. Snapshots are taken from a slice at the center of the sample, perpendicular to the z-direction. Video S1 shows the same process as a movie.
bles is affected by varying the side lengths of the sample, thus isolated type-II bubbles are not observed in infinite systems. Moreover, they are only observed up to a field of 0.26 T, which might change by varying the saturation magnetization or the exchange constant values. In contrast, in Figure S2 we observe that type-I bubble arrays survive in a larger field range. When decreasing the system thickness, we saw from simulations that type-II bubbles are stable in a more restricted range of applied field magnitudes. Figure S4 shows the effect on a simulated magnetic bubble of changing the specimen thickness. It can be seen that the bubbles in thinner specimens are smaller and more elliptical with their major axis following a line connecting the internal domain walls. All the bubbles share the feature that the magnetization points radially at the top and the bottom layers and the magnetization twists throughout the specimen thickness, showing the classic type-II configuration only in the middle layer.
All the simulation data can be reproduced from the repository given in Ref. [10].  Figure S5(a) shows an image from a defocus series with stripe domains and magnetic bubbles coexisting. The sample was cut to two different thicknesses with a thickness step running nearly vertically down the middle. The thickness was measured as described in section S4. The thinner region on the left appears lighter and has a gradual increase in thickness from 110 nm at the top of the image to 180 nm at the bottom. The thicker region on the right is 200 nm thick at the top and increases to 230 nm at the bottom. The straight lines parallel to the thickness step in the thinner region are slight thickness undulations of about 1 nm caused by the focused ion beam milling used to prepare the specimen.
The bubbles have a mean spacing of 220 nm in the thinner region and 310 nm in the thicker region, both with a standard deviation of 25 nm. They have a mean diameter of 72 nm in the thinner region with a standard deviation of 13 nm and a mean diameter of 117 nm with a standard deviation of 6 nm in the thicker region. This follows the trend shown by the simulations that the bubbles in thicker specimens tend to be larger.
The bubbles in both regions exhibit a similar Yin-Yang appearance although the black-white contrast of the bubbles in the thinner region is rotated with respect to those in the thick region. Simulations show that this rotation reflects the positions of the internal domain walls. A closer inspection shows that within the same region, neighboring bubbles show small rotations with respect to one another. The transport of intensity equation was used to reconstruct the projected B-field. The B-field from the bubble in the thinner region outlined in blue is shown in (b) and the B-field from a bubble in the thicker region outlined in red is shown in (c). The bubble from the thinner region is from a thickness of 180 nm and has a diameter of 87 nm while the bubble from the thicker region is from a thickness of 220 nm and has a diameter of 117 nm. It can be seen that the pattern of B-field lines is very similar in both cases and closely resembles the B-field in figure  3  type-II.

A. Simulation of Defocused Images
Once a model for the magnetization in a specimen has been obtained, electron micrographs can be simulated as follows. The first step is to find the magnetic vector potential. This can be done by splitting the specimen into an infinite number of magnetic dipoles and summing the vector potentials from each. Thus the vector potential A(r) can be related to the magnetization M(r) via: (1) where r and R are 3 dimensional position vectors and µ 0 is the permeability of free space. As this is a convolution, it is treated most simply as a multiplication in Fourier space [11]. Taking the 3-dimensional Fourier transform gives: Taking Cartesian coordinates in which the initial electron beam direction is parallel to z, the phase shift the electron beam has experienced when it exits the specimen is related to the vector potential via the Aharanov-Bohm formula [12]: In Fourier space, this becomes: and so the phase can be calculated directly from the magnetization using the formula given by Beleggia and Zhu [11]: where d is the specimen thickness, k ⊥ ≡ (k x , k y , 0) and M ⊥ (k x , k y ) is the two-dimensional Fourier transform of the magnetization averaged along z through the thickness of the sample. The in-focus wavefunction is then ψ 0 = exp[iφ] and it can be seen immediately that an in-focus image of a purely magnetic specimen will be featureless as I 0 = |ψ 0 | 2 = 1. One way to obtain an image showing magnetic contrast is to take images out of focus. This technique is called defocused or Fresnel imaging and is often referred to as Lorentz imaging or LTEM although the latter term should encompass all the magnetic imaging techniques available with an electron microscope rather than this specific method. The wavefunction for a defocused image ψ ∆f can be calculated using the Fresnel-Kirchoff integral [13,14] which describes the propagation of waves through free space and can be written as a convolution: where r 2 = x 2 + y 2 and * indicates a convolution. The calculation of the defocused image is usually performed in Fourier space as the convolution becomes a multiplication. Thus The defocused image is then I ∆f = |ψ ∆f | 2 .

B. Projected B-field
It is often of interest to calculate the projected B-field, B ⊥ . This is the component of the magnetic flux density normal to the electron beam averaged along the trajectory of an unscattered electron beam. If Stokes' theorem is applied to the Aharanov-Bohm equation (Eqn. 3) and the magnetic flux density substituted for the curl of the vector potential, the resulting integral can be inverted to give the projected B-field as:

C. Analytic Model for a Type-I Magnetic Bubble
Thus, if we know the magnetization, the method described above allows us to simulate images that would be acquired using transmission electron microscopy. Here we describe the analytic model for the magnetization of a magnetic bubble used for comparison with experimental images. In this model, it is assumed the magnetization does not to vary through the specimen thickness and so represents the magnetization averaged through the thickness along the electron beam direction. Using a Cartesian coordinate system in which the axis of the bubble is parallel to the z-axis, a standard model [11,15] for a circular magnetic domain wall has in-plane magnetization: and z-component where r is the distance from the bubble's axis and θ is the azimuthal angle. The constants are: M s the saturation magnetization, R the radius of the magnetic bubble and δ the domain wall width. The domain wall width is given by δ = π A/K with A being the exchange energy constant and K the uniaxial magnetocrystalline anisotropy. If the magnetization direction does not follow the domain wall but deviates from it by a constant angle θ d , the z-component remains the same whereas the in-plane component becomes: A purely Bloch bubble is produced when θ d = 0 and a Néel bubble is produced when θ d = 90 • .

D. Field Lines for a Magnetic Bubble
It is beneficial to calculate the shape of the field lines for the magnetization of a magnetic bubble as an internal domain wall will likely follow a field line of M to avoid generating magnetic poles. Coey [16] (p 28 eqn. 2.11) defines a field line as a line which solves Expressing Eqn. 11 in polar coordinates gives So the magnetic field lines are the solution to Solving this differential equation gives the equation of a logarithmic spiral [17] so that a field line going through a point (R, θ i ) has equation:

E. Model for a Modified Type-II Magnetic Bubble
A type-II bubble can be modeled by adding a reversing function θ rev to Eqn. 11: The reversing function adds π to the angle the magnetization makes with the radius between the two internal domain walls which occur at angles θ 1 , θ 2 . If the internal domain walls run radially, the reversing function can be written as: where the arctangent of the hyperbolic sine is the function describing the magnetization direction in a domain wall given by Coey [16] and δ i is the internal domain wall width.
The magnetization of each of the internal domain walls can point towards or away from the center of the bubble. In combination with this, the magnetization of the bubble can circulate in one of two directions leading to eight possible domain wall arrangements. Each of these can be selected by reversing the sign of M s in Eqn. 16 in combination with choosing the sign of the arctangent functions in Eqn. 17.
When θ d = 0, this model represents a conventional type-II bubble. When θ d is non-zero, it may be that the internal domain wall follows a field line of magnetization to avoid generating magnetic poles and thereby minimize stray field energy. The equation for the field lines (Eqn. 15) can be used so that if the walls go through the points (R, θ 1 ) and (R, θ 2 ), the reversing function becomes: θ rev = ± atan sinh πr δ i (θ − θ 1 − cot(θ d ) ln(R/r)) ± atan sinh πr δ i (θ − θ 2 − cot(θ d ) ln(R/r)) (18)

F. Modified Type-II Magnetic Bubble with an Outer Domain Wall
A comparison with experimental images indicated that the bubble had an outer domain wall. To accommodate this, we allowed the deviation angle θ d to vary between two limits: θ i near the center and θ o far from the center with a domain wall at radius R o and width δ o . Thus we have the following model for a magnetic bubble which is the analytic model used in the main text for comparison with the experimental images: where θ rev = ± atan sinh πr δ i (θ − θ 1 − cot(θ d ) ln(R/r)) ± atan sinh πr δ i (θ − θ 2 − cot(θ d ) ln(R/r)) (20) and This model was compared with the experimental images of 19 bubbles as described in the main text, yielding the parameters given in Table S1.