Imaging Nucleation and Propagation of Pinned Domains in Few-Layer Fe5–xGeTe2

Engineering nontrivial spin textures in magnetic van der Waals materials is highly desirable for spintronic applications based on hybrid heterostructures. The recent observation of labyrinth and bubble domains in the near room-temperature ferromagnet Fe5–xGeTe2 down to a bilayer thickness was thus a significant advancement toward van der Waals-based many-body physics. However, the physical mechanism responsible for stabilizing these domains remains unclear and requires further investigation. Here, we combine cryogenic scanning diamond quantum magnetometry and field reversal techniques to elucidate the high-field propagation and nucleation of bubble domains in trilayer Fe5–xGeTe2. We provide evidence of pinning-induced nucleation of magnetic bubbles and further show an unexpectedly high layer-dependent coercive field. These measurements can be easily extended to a wide range of magnetic materials to provide valuable nanoscale insight into domain processes critical for spintronic applications.


S1 Sample properties and preparation
An iron-enriched mixture consisting of Fe, Ge and Te at a 6:1:2 ratio was sealed in an ampule and held at 700 • in a tube furnace, yielding single crystals of Fe 5 GeTe 2 after slowly cooling down after a few days.An average ratio of Fe:Ge:Te = 4.6:1:2.1 was extracted for the stoichiometry, but it is noteworthy that the stoichiometry differs across the as-grown bulk crystal.The Curie temperature of the bulk crystal was determined to be T C = 260 K using vibrating sample magnetometry (VSM, see Fig. S2).While it is possible to exfoliate Fe 5 GeTe 2 down to about ten layers, it becomes increasingly difficult to obtain flakes thinner than that.Here, we exploit the interaction between gold (Au) and the top-most layer of the Fe 5 GeTe 2 crystal to exfoliate layers as thin as a monolayer.The sample was originally capped by a thin layer of amorphous selenium (Se) to avoid degradation.From experience, this Se cap typically also results in small islands, which can be moved easily with a scanning probe.To reduce the risk of contaminating the diamond tip with Se, the sample was additionally capped with a 20 nm thick and homogeneous flake of hexagonal boron nitride (hBN), indicated by the yellow outline in Fig. S1(b).In the 2D community, hBN is routinely used as a capping layer to avoid degradation.It is free of dangling bonds and very hard, which makes it an ideal surface for the diamond tip.Keeping the hBN thickness to 20 nm or below ensures that the NV's spatial resolution is not overly compromised.

S2 Physics and limitations of NV sensing
The NV center is an atomic-sized defect in diamond, consisting of a substitutional nitrogen atom next to a lattice vacancy.The ground state manifold of the NV forms a spin triplet with spin sublevels m s = 0, ±1 and the quantization axis (NV axis) along one of the four crystallographic axes of diamond.The crystal field splits the m s = 0 from the m s = ±1 by 2.87 GHz.
An additional Zeeman splitting between the m s = −1 and m s = +1 levels is induced by an external magnetic field along the NV axis, B NV .In the weak field approximation, the frequency splitting, ∆ν is linear and given as 2γB NV , with γ being the gyromagnetic ratio of 28 MHz/mT.Optical excitation results in a spin-state dependent fluorescence, depending on whether the NV spin population is predominantly in the m s = 0 (bright) or m s = ±1 sublevel (dark).Combined with ground state spin control via microwave excitation, one can optically detect the electron spin resonance, which is commonly known as optically detected magnetic resonance (ODMR).In order to retrieve the field polarity, we apply a small bias external field, B bias , along the NV-axis.Specifically, this allows us to extract the frequency splitting ∆ν to determine the total local field (including the bias field) projected on the NV-axis, given as |B total | = B NV + B bias .The actual field component due to the sample is then given as B NV = |B total | − B bias .The stray field maps in the main text are obtained by recording and fitting a full ODMR spectrum at one of the two NV transition frequencies at every pixel of the image.We adopt a pulsed ODMR protocol consisting of a microwave π-pulse and a subsequent laser pulse for combined readout and spin initialization.A wait time of 600 ns before the next microwave pulse ensures relaxation of trapped population towards the ground state. 1 Imaging based on ODMR is limited to magnetic fields below ∼ 10 mT due to spin-mixing induced by the orthogonal component of the magnetic field.This results in an inefficient spin polarization that leads to a decreased NV fluorescence as well as ODMR contrast. 2,3dard diamond probes typically have NV centers pointing at an oblique angle with respect to the out-of-plane z-axis, which renders the study of few-layer van der Waals magnets with strong perpendicular magnetic anisotropy (PMA) but low moments in a high out-of-plane (OOP) field environment extremely challenging.While prototype diamond probes with NV centers parallel to the z-axis are now available and could prove useful, 4 imaging at Tesla fields is still challenging from an instrumentation point of view.The microwave frequencies required for ODMR very quickly approach hundreds of GHz, requiring special hardware and microwave delivery solutions to the NV.

S3 Diamond quantum microscopy setup
The diamond quantum microscope (DQM) is an integrated confocal and atomic force microscope housed in a closed-cycle cryostat (attoDRY1000, Attocube Systems).The confocal optics is home-built and the atomic force microscope platform is based on an electrically readout tuning fork.All measurements are conducted at 4 K, unless specified otherwise.The NV center is optically excited and read out using a pulsed-ODMR protocol which tracks the spin resonances for improved sensitivity and reduces heat build-up.We utilize diamond scanning probes with a single NV center, implanted with an energy of 7 keV, at the apex (QZabre AG) for imaging.Microwaves for ODMR measurements are delivered via a 20 µm thick copper wire mounted close to the Fe 5 GeTe 2 flake.More details of the setup can be found in Ref. 5

S4 Sensor characterization S4.1 NV axis orientation
The orientation of the NV axis (the axis joining the vacancy and the nitrogen atom) with reference to the laboratory frame is characterized by monitoring the Zeeman splitting of the spin resonances as a function of applied magnetic field ⃗ B (B 0 , ϑ B , φ B ), using a three-axis Helmholtz coil.The maximum Zeeman splitting occurs when the field is parallel to the NV axis.Holding the field at a constant magnitude, we obtain ODMR spectra while varying the polar ϑ B (azimuthal φ B ) as shown in Fig. S3, resulting in ϑ NV = 120

S4.2 NV-to-sample distance
We infer the NV-to-sample distance d NV by measuring the magnetic field produced by the edge of a uniformly magnetized stripe. 6The calibration sample consists of a film of CoFeB patterned into magnetic wires, which remain saturated at remanence.A single linescan showing the position of the upper spin resonance extracted from ODMR measurements, along with the topography of the edge, is shown in Fig. S4(a).The red curve is a fit to the data, which directly gives the value for d NV .A statistical average value of (60 ± 5) nm is obtained from 20 linescans, as shown in Fig. S4(b).Note that an additional 20 nm needs to be accounted for due to the thickness of the hBN capping layer, yielding a total NV-to-sample distance of d NV ∼ 80 nm.
Figure S4: Zeeman shift of the upper spin resonance as a function of position above the magnetic wire, as indicated by the topographic height of the edge.An average NV-to-sample distance of (60 ± 5) nm has been extracted from 20 line scans across the edge.

S5 Reverse propagation 82
The magnetostatic problem of retrieving magnetization from a magnetic field distribution is in general ill-posed.However, in the limit of a thin magnet and a predominantly OOP magnetization (i.e., few layer Fe 5 GeTe 2 ), the problem is sufficiently constrained to produce a unique solution.Assuming a planar magnetization distribution M that is constant throughout its thickness t, the stray field B generated at a distance d from the surface is related to the magnetization in Fourier space (k x , k y ) as: with the dipolar tensor given as: Here, b(k x , k y , z) and m(k x , k y , z) are the spatial Fourier transforms of B(x, y, z) and M(x, y, z), respectively. 7Indeed, the non-invertible matrix in eq. ( S2) is a direct indicator that the pla-nar magnetization distribution cannot be uniquely reconstructed from the measured stray field.However, if we assume a strong perpendicular magnetic anisotropy in our system, such that M x,y = 0 and M z ̸ = 0, the problem is then well-posed and a unique solution can be found. 8,9The M z component can then be directly reconstructed from the z-component of the magnetic field: ) d+t) and the magnetic field projected on the NV axis with reference to the lab frame can be expressed as For the reverse propagation of the stray field maps in the main text, an additional Hann filter is used to suppress noise. 8,10From the reconstructed OOP magnetization, we can

S6 Size analysis of ZFC morphology
The size of domains in 2L and 3L Fe 5 GeTe 2 as observed after zero field cooling, was determined using two separate approaches, respectively.For the bubble morphology in 3L, a thresholded map was separated in regions, each containing a single bubble domain (see Fig. S6(a)).Since the shapes vary drastically between domains, we chose the major axis of each domain as a common denominator, yielding a mean value of 350 nm (Fig. S6(b)).
Rather than an absolute size, we extracted a domain period from 2L Fe 5 GeTe 2 using a 2D autocorrelation.The autocorrelation map of 2L Fe 5 GeTe 2 after ZFC is shown in Fig. S7, including the line profiles to extract the peaks, which in turn correspond to a high degree of similarity.For each line profile, the mean difference between peaks is extracted and plotted

S7 Domain pinning
We visualize the domain pinning observed in Fig. 3(b) of the main text by spatially quantifying the pinning probability.For the specific example of Fig. 3(b), an area is identified as pinned when it remains unpolarized (blue) while its surrounding is polarized (red) by the applied magnetic field.Since the sample undergoes a magnetic reset via a negative saturation field for each of the image in Fig. 3(b), we can repeat the analysis above on every image to collate the total number of times an area gets pinned.The probability of an area to get pinned is then extracted by normalizing the total pinned occurrences by the number of instances the area was identified polarized (red).Figure S6 shows the heatmap of pinning probability over the sampled region in Fig. 3(b).It strongly indicates that high pinned probabilities are correlated with areas of smaller size, and that pinning is distributed across the magnetic material.Here we outline the steps for the simulation of the magnetic bubbles as a function of helicity as well as the minimization of the cost function.Neglecting the contribution of the domain walls to the stray field, we can directly reconstruct M z from the stray field map in Fig. 1(c) in the main text.The average domain wall width is found by extracting line profiles across all domain walls visible in Fig. S9(a).The average normalized magnetization is then fitted to the standard analytical description of a domain wall: is then binarized via Otsu thresholding (Fig. S9(c)) and domain boundaries are detected using the Canny edge algorithm (Fig. S9(d)).Finally, the edge normal vectors are found and at every edge pixel, and a domain wall is created along each normal vector, with s 0 being the center of the domain wall.The magnetic field generated at a distance d above the simulated magnetization distribution can be calculated using the dipolar tensor introduced in section S5.We then minimize the sum of the squared differences between the experimentally obtained B exp NV and the simulated B sim NV with respect to the helicity ξ, which can be expressed with a cost function: where N x and N y are the number of pixels in x and y, respectively, while fixing the NVto-sample distance d NV = 80 nm, the domain wall width w DW = 137 nm, and the surface magnetization at I s = 66 µ B /nm 2 .

S9 Details on domain wall width extraction
We validate our approach of extracting the domain wall width from bubble domains by simulating m z with a fixed wall width of 135 nm based on the bubble morphology displayed in the main text and SI sections above (see Fig. S10(a).We extract a domain wall width of (126 ± 31) nm, close to the initial simulated value.The error of 31 nm can be attributed by some of the line profiles not being perpendicular to the wall.

S10 Robustness of domain wall analysis
In the domain wall analysis section we assumed fixed parameters for the NV-to-sample distance d NV (80 nm) and the domain wall width w DW (135 nm).Since both parameters are associated with a certain error, we show here that we can still reliably distinguish between Bloch and Néel walls even at the given error boundaries.To this end, we constructed a minimization matrix on the basis of the analysis explained in the main text and S7 and minimized the squared differences between experimentally obtained B exp NV and the simulated

Figure S1 :
Figure S1: Sample characterization.(a) Optical image of the exfoliated flake capped with hexagonal boron nitride (hBN), which is outlined in yellow.(b) AFM line scans across the edges, indicated by the colored dashed lines in (a), confirming the flake thicknesses of 1L, 2L and 3L.

Figure S2 :
Figure S2: Crystal strucutre and VSM of bulk Fe 5 GeTe 2 .(a) Side view of the average Fe 5 GeTe 2 crystal structure with Fe(1) and Ge being split sites.(b) M -vs-T measurement of the bulk Fe 5 GeTe 2 crystal using a vibration sample magnetometer, yielding a T C of 260 K.

Figure S3 :
Figure S3: Zeeman splitting of the upper NV spin transition.(a) As a function of φ B with ϑ B fixed at 90 • , and (b) as a function of ϑ B with φ B fixed at 96 • .The maxima in the curves determine the NV orientation, giving ϑ NV = 120 • ± 4 • and φ NV = 96 • ± 4 • .

Figure S5 :
Figure S5: Reconstructed OOP magnetization distribution of the remanent state.The panels represent the distribution after saturating with (a) 7 T, (b) −4 T, (c) 4 T and (d) 3 T.The red curves are Gaussian fits to retrieve the peak values.

Figure S6 :
Figure S6: 3L bubble domain size analysis.(a) Separation of thresholded map into regions, each containing a bubble domain.The red lines indicate half of the major axis.(b) Histogram of values extracted from (a) with a mean value of 350 nm indicated by the red dashed line.

Figure S7 :
Figure S7: 2L domain analysis.(a) Autocorrelation map of Fig.1(d) of the main text and line profiles used to extract the peak positions.(b) Autocorrelation line profiles corresponding to the lines in (a).(c) Single values for the domain period extracted from each line profile and the mean indicated by the blue line.
Figure S9(b) illustrates extracted line profiles (black line) and the fit to the profile average (red line) which gives a mean domain wall width of ∼137 nm.The map in Fig. S9(a)

Figure S10 :
Figure S10: Validating domain wall width extraction.(a) Simulated m z map based on the bubble domain morphology with a fixed wall width of 135 nm and the yellow line cuts to extract the width.(b) Resulting domain wall width of (126 ± 31) nm, the red curve is the average line profile.
Fig.S11, indicating an average helicity of π/2, hence suggesting that the resulting helicity is not greatly impacted by varying d NV and w DW in the parameter space associated with our data and analysis method.Specifically, with a spread of 70 nm to 100 nm for d NV and of 60 nm to 200 nm for domain wall widths, the resultant average helicity is bounded between 0.4 and 0.68.Here, helicity holds a value from 0 to 2 which wraps around to 0, and perfectly

Figure S11 :
Figure S11: Minimization matrix.Minimization of the squared differences between experimentally obtained B exp NV and the simulated B sim NV with respect to the helicity ξ, as outlined in S8 for fixed d NV and w DW .

Figure S12 :
Figure S12: Pure Bloch domain walls.Top: Experimental stray field map compared to the simulated stray fields generated by domains with pure Bloch walls (ξ = π/2).Bottom: Squared error between experimental and simulated field maps and the magnetization used to generate the simulated stray fields.

Figure S15 :
Figure S15: Domain formation after second cooling.Stray field map after second cool down (left) at the same location as after first cool-down (right) showing overlap of the down domains, strongly indicating domain formation due to pinning sites..

Figure S16 :
Figure S16: Stray field maps.Raw stray field maps used to reconstruct the m z maps in the AC demagnetization section of Fig.2(b)-(e).