History-dependent domain and skyrmion formation in 2D van der Waals magnet Fe3GeTe2

The discovery of two-dimensional magnets has initiated a new field of research, exploring both fundamental low-dimensional magnetism, and prospective spintronic applications. Recently, observations of magnetic skyrmions in the 2D ferromagnet Fe3GeTe2 (FGT) have been reported, introducing further application possibilities. However, controlling the exhibited magnetic state requires systematic knowledge of the history-dependence of the spin textures, which remains largely unexplored in 2D magnets. In this work, we utilise real-space imaging, and complementary simulations, to determine and explain the thickness-dependent magnetic phase diagrams of an exfoliated FGT flake, revealing a complex, history-dependent emergence of the uniformly magnetised, stripe domain and skyrmion states. The results show that the interplay of the dominant dipolar interaction and strongly temperature dependent out-of-plane anisotropy energy terms enables the selective stabilisation of all three states at zero field, and at a single temperature, while the Dzyaloshinksii-Moriya interaction must be present to realise the observed Néel-type domain walls. The findings open perspectives for 2D devices incorporating topological spin textures.

linear, increase of K with decreasing temperature. In the mean-field simulations, we utilised a similar scaling of K with temperature to recover the temperature dependent size of the observed magnetic domains, which at low enough temperatures leads to uniform switching of the simulated sample, as was observed in the real-space measurements of the FGT flake.
Having characterised the bulk material, we proceeded to fabricate an FGT flake sample using the mechanical exfoliation technique. The optical microscope images in Fig. S3a

Supplementary Note 2: Flake thickness determination
We utilised a combination of AFM and x-ray absorption analysis to estimate the thickness of each region of the FGT flake. Fig S4a shows an AFM image of a corner of the investigated FGT flake.
Cuts across the steps of three thickness regions are labelled S1-3. The profile of each step is shown in Fig. S4b, revealing thicknesses of approximately 60, 30 and 15 nm, as well as a measure of the hBN capping flake thickness of roughly 15 nm. A STXM image of the flake is shown in Fig.   S4c, with a number of regions numbered 1 to 5 in order of decreasing thickness, including a 6th point passing through only the Si 3 N 4 membrane, which acts as a measure of the incident beam intensity, or I 0 signal. The ROI investigated in the main text is indicated, and the thicknesses of the regions measured in the AFM data are labelled. Absorption spectra measured locally in each thickness region, across the L 3 edge, are plotted in Fig. S4d. We extracted the final data point in these spectra, measured away from the resonant edge at 715 eV, and divided them by the I 0 value to determine the x-ray transmission, T , of each thickness region. These transmission values are plotted on a logarithmic scale as a function of sample thickness in Fig. S4e, where the solid line is a linear fit to Beer's law. Using this calibration, we determined that the three thicknesses featured in the ROI were approximately 60, 50 and 35 nm respectively.

Supplementary Note 3: X-ray absorption spectra results and oxide thickness
To investigate the oxidisation of our samples, x-ray absorption spectra of both the bulk and flake FGT samples were measured. In Fig. S5a we show the Fe L 3 and L 2 edge absorption of the bulk FGT sample measured using the surface-sensitive total electron yield (TEY) method (see methods). The spectra were processed by first rebinning the measured data to a 0.1 eV energy resolution. Then, the data was divided by the I 0 (the measured intensity of the synchrotron beam at that energy) measurement to isolate the signal from the sample. Next, a background was subtracted by performing a linear fit to the first 5 eV of the data set. Finally, the spectra was edge-normalised, as is standard practice, by dividing to set the after-edge value to unity. Spectra were measured for both an initially uncleaved surface, and for a surface cleaved under vacuum conditions. The large secondary peak in the uncleaved spectrum reveals the change in the valence state of the Fe 3+ to Fe 2+ due to reaction with oxygen under ambient conditions. The uncleaved surface showed no XMCD signal at 80 K, indicating oxidised FGT is non-magnetic.
We also measured absorption spectra at the O K edge for both the cleaved and uncleaved samples, as shown in Fig. S5b. The lack of absorption at the O K edge demonstrates that by cleaving the sample in vacuum, almost all oxidised FGT material can be removed. Comparison of the results in Fig. S5a and b, confirms that the large secondary absorption peak at the Fe L 3 edge in the FGT is due to oxidisation of the Fe.
In the cleaved sample spectrum, a strong XMCD signal was observed by measuring spectra at oppositely applied magnetic fields, as shown in Fig. 5c. Interestingly, there are still two clear peaks in the L 3 absorption edge, which we attribute to the two Fe ion species present in the FGT structure, commonly referred to as the Fe I and Fe II sites 1 . Fig. 5c also shows the ad hoc step function utilised to subtract the background from excitations into continuum states, where the height of the steps at the L 3 and L 2 edges was set at a ratio of 2:1, according to the occupation of the 2p 3/2 and 2p 1/2 core states. Having subtracting this non-resonant background, the resonant contribution of the absorption spectrum, and its integrated value, is plotted in Fig. 5d. Finally, the integrated XMCD spectrum is plotted in Fig. 5e. Using sum rules, the spin and orbital contributions of the magnetic moment, m s and m o , can be determined from the integrated resonant absorption and XMCD spectra 2 . We recovered values of 1.52 µ B and 0.05 µ B for m s and m o respectively, comparable to other theoretically and experimentally determined values found in the literature 2, 3 .
These corresponds to a value of approximately 375 kA/m in the bulk sample, which is comparable to that measured in the magnetometry data.
Spectra from ∼10 and 60 nm regions of the flake sample are shown in Fig. 6a and b, mea-sured at 100 K. The lack of XMCD effect in the 10 nm region illustrates the absence of magnetic order in this thickness, while the large secondary peak indicates significant oxidisation of the sample, as was seen in the uncleaved bulk spectrum. On the other hand, the spectra measured in the 60 nm region shows a strong XMCD signal, indicating magnetic order, and a significantly reduced second Fe absorption peak, suggesting less overall oxidisation of the sample. Given that all thickness regions of 15 nm and below were observed to be non-magnetic, and since the oxidisation of the sample is likely to be uniform across the entire surface of the flake, we can estimate that To directly investigate the oxidised FGT layer thickness, we performed standard transmission electron microscopy (TEM) to examine the thickness of the oxide layer in the FGT flakes fabricated by our exfoliation method in ambient conditions. Due to the challenge of acquiring a TEM lamella using focused ion beam from a flake sample stamped on a thin Si 3 N 4 membrane, we prepared an additional flake sample, capped by hBN, prepared on a SiO 2 substrate. Besides the change of substrate, the fabrication method was identical, with the flake being exposed to ambient conditions for about 30 mins. Focused ion beam was utilised to prepare a lamella of a cross section through this flake stack. Layers of carbon and Pt were deposited before the milling and lift out procedure to protect the underling structure of the sample from ion damage. The resulting TEM images are shown in Fig. S7a and b, where the square root of the intensities has been shown to improve the contrast of each material layer. By the lightening of the contrast, indicating lower average atomic number, and the breakdown of the clean vdW layers, within the FGT section, the oxidised FGT layer thickness can be estimated to be about 7 nm at the top and bottom of the flake. It is clear that the oxidised FGT does not form a clean interface with the pristine FGT. This can be expected, since one might assume that the penetration of the oxygen into the bulk of the FGT would be some form of exponentially decaying process.

Supplementary Note 4: Extended STXM imaging data
We performed STXM imaging following the zero field-cooling (ZFC) procedure to map out the phase diagrams presented in the main text. The results when ZFC to 150 and 125 K, and the subsequent decrease in field, are shown in Fig. S8a and b respectively. One can see that the stripe domains were frozen-in as the sample was cooled. After annihilating the stripes by saturating the magnetisation at these temperatures below 150 K, further field sweeps resulted in the entire sample switching simultaneously, as in the field-sweep protocol. We found it was possible to freeze-in domains at 0 mT in this way down to 30 K -the lowest temperature achievable in our STXM instrument. Beyond those presented in the main text, we acquired STXM images at further temperatures following the field sweep procedure, as shown in Fig. S9. In the case of the data at 186 K, only the 60 nm thickness region showed visible magnetic domain contrast, and so we focused our imaging measurements solely on this region.
In the main text, we present the majority of the STXM images after normalising the contrast within each thickness region. This allows the magnetic contrast to be discerned more clearly, by eliminating the structural contrast due to the different thickness layers. The contrast normalisation was performed by isolating each thickness region of the flake using a mask, as shown in Fig. S10ac. By summing these masked images together, we achieve the result shown in Fig. S10d, where the signal seen within each thickness region is primarily out-of-plane magnetic contrast.

Supplementary Note 5: Extended LTEM data
We performed Lorentz transmission electron microscopy (LTEM) on a separate FGT flake sample to investigate whether the magnetic textures were Bloch or Néel type. This FGT sample was prepared in a similar manner to those used for x-ray imaging, but the hBN capping layer was not used, since it introduced additional structural contrast to the images. The uncapped FGT sample was kept in an Ar atmosphere between preparation and the LTEM measurement to reduce its exposure to ambient conditions. The inset of where m z is the observed magnetic contrast, and c is the non-XMCD signal background count rate.

Supplementary Note 7: Field cooled skyrmion fitting
To extract the number and size of the skyrmions from the images obtained after FC, fitting was performed using the OpenCV Python Library 4 . Fig. S15a shows the contrast-normalised image after FC at 15 mT from above T C to 125 K. To begin, thickness regions were isolated in the raw data using an intensity threshold, as shown for the 60 nm region in S15b. Next, the edges of the skyrmions were located using a Canny edge detection algorithm, resulting in the output shown in Fig. S15c. These output fits were refined by noise filtering and manually excluding ill-fitting contours. The contour functionality of the OpenCV library was used to fit contours to the bubble edges, and extract their areas -the results for each thickness region are shown in S15d and e. We where a is the discretisation constant and m(T ) average magnetic moment of spins. We will take a to be equal to the Bloch parameter l w = A/K, which describes the smallest characteristic domain wall size in the material. It assumes small values for hard magnetic materials (e.g. l w = 1.3 nm for bulk Nd 2 Fe 14 B) and large values for soft magnetic materials (e.g. l w = 18 nm for bulk Fe) 6 .
Fe 3 GeTe 2 intuitively falls between these limits and we will take a = 5 nm as a reasonable estimate (similar to bulk Co).
The micromagnetic DMI energy constant was estimated in the main text based on experiments to be D = 0.12 × 10 −3 J/m 2 . Therefore, by using the first two relations in Eq. tion. Therefore, it is expected that they differ from the parameters utilised in the micromagnetic simulations, where we are simulating skyrmions at a higher temperature, and therefore require a smaller value of the anisotropy.

Supplementary Note 9: Temperature dependence of anisotropy in the mean-field model
To describe the temperature dependence of anisotropy in the mean-field model we first considered the classical Callen-Callen theory 7 , based on which we represented uniaxial anisotropy constant associated with individual magnetic moments as: Thus, Eq. (4) incorporates a piece-wise continuous linear trend in the multiplicative prefactor to enhance the rate of increase of anisotropy at low temperatures below threshold temperature β 1 , where β = 1/k B T , as demonstrated in Fig. S18. Although the form of Eq. (4) represents only a simple linear extension of Eq. (3) and is chosen purely on empirical basis, it reproduced the experimentally observed features very well as can be seen in Fig. 6 in the main text. To obtain Fig. 6, the free parameters in Eq. (3) were identified by a trial-and-error process to be equal to Finally, it turned out that to describe the observed experimental data trends it was also necessary to enhance the temperature-dependence of exchange in Eq. (1) as: To obtain Fig. 6, the free parameters required in Eq. (5) were again identified by a trial-and-error process to be equal to c 2 = 4.167J E (0)β C with J E (0) = 0.7.   ,d), to the c axis of the bulk Fe 3 GeTe 2 crystal. Measurements were performed either by fieldcooling (FC, orange) from 300 K to 5 K under an applied field, or by initialising the sample by zero field-cooling (ZFC, purple) and subsequently measuring from 5 K to 300 K under an applied field, as indicated by the coloured arrows. Data was measured for applied fields of both 30 Oe (a,c) and 300 Oe (b,d). e-h, Magnetisation versus field measurements carried out at a selection of temperatures, for both H ∥ c (purple) and H ⊥ c (orange) field alignments. i,j, Values of the uniaxial anisotropy, K, and the saturation magnetisation, M s , extracted from the magnetisation versus field measurements, plotted as a function of temperature, T .  Step height measurements in the sample locations S1, S2, S3 as labelled in d, and the thickness of the hexaboron nitride (hBN) layer. c, X-ray microscopy image of the FGT flake with the thickness regions measured by the AFM labelled. Further thickness regions are labelled 1 to 6. d, X-ray absorption spectra measured in each of the 6 thickness regions labelled in c across the Fe L 3 edge. e, The transmission of the flake is plotted on a log scale as a function of the determined sample thickness, and fitted with Beer's law, providing a calibration to determine the thickness of all regions. Figure S5 | Bulk sample x-ray absorption spectra measured with total electron yield. a, Xray absorption spectra measured at the Fe L 3 and L 3 edges on the uncleaved (orange) and cleaved (purple) Fe 3 GeTe 2 (FGT) bulk single crystal at 293 K. b, X-ray absorption spectra (XAS) measured at the O K edge on the uncleaved (orange) and cleaved (purple) FGT bulk single crystal at 293 K. c, X-ray absorption spectra measured at the Fe L 3 and L 3 edges on the cleaved sample at an applied field of ± 4 T (orange and purple) at 80 K. Subtraction of the two data sets gives the x-ray magnetic circular dichroism (XMCD) signal (green). d, Step-subtracted (black), and integrated (green), xray absorption spectrum. e, The integrated XMCD signal from c. Crosses indicate points in the integrated spectra utilised in the sum rules analysis. Figure S6 | Flake sample x-ray absorption spectra measured in transmission. a-h, X-ray absorption spectra (XAS) measured at the Fe L 3 and L 3 edges measured for the 15 and 60 nm thickness regions of the Fe 3 GeTe 2 (FGT) flake sample. Spectra measured at ± 250 mT are plotted (purple and orange) alongside the calculated x-ray magnetic circular dichroism (XMCD) signal (green) for the 15 and 60 nm thickness regions at a range of temperatures. At 125K and below, the applied field was no longer sufficient to switch the magnetisation of the sample, so the circular polarisation (C±) was selected instead. Figure S8 | X-ray microscopy imaging results following the zero field-cooling procedure. a,b, X-ray microscopy images measured following the zero-field cooling (ZFC) measurement procedure at 150 K and 125 K respectively. The green arrow indicates the measurement direction of the applied magnetic field. The scale bar is 1 µm. Figure S9 | Additional data sets for the field sweep x-ray microscopy measurements. a-o, Xray microscopy images measured following the field sweep measurement procedure at 160 K, 170 K and 186 K. The green arrow indicates the measurement direction of the applied magnetic field. The scale bar is 1 µm. igure S11 | Lorentz transmission electron microscopy images of a tilted flake. a, Lorentz transmission electron micrograph (LTEM) of an Fe 3 GeTe 2 (FGT) flake acquired at 92 K in zero applied magnetic field at a defocus of 5.11 mm. Instead of lying normal to the electron beam, the sample was tilted to 36°about the axis shown. The inset shows an optical micrograph of the flake mounted on a Si 3 N 4 membrane and the area from which the LTEM image was acquired is outlined in red. b, An enlargement of the image in a from the area outlined in green. c, An image from the same area as b acquired under the same conditions but the sample was tilted to 10°about the same axis.   Figure S16 | Formation of skyrmion after field-cooling. a-d, X-ray microscopy images of the FGT flake acquired after field-cooling (FC) the sample from 200 K to 125 K at different applied magnetic fields. The histograms plot the distribution of the skyrmion size in the 60 nm (purple), 50 nm (green) and 30 nm (yellow) thickness regions, as determined by image recognition software. The scale bars are 1 µm.
Figure S17 | Mean-field simulations without modified anisotropy scaling. a-e, Mean-field simulations of the magnetic domain structure evolution following the field sweep procedure for increasing applied field, H, at different temperatures, utilising the unmodified Callen-Callen anisotropy scaling model.  Figure S18 | Temperature dependence of the mean-field effective anisotropy. The effective anisotropy utilised in the mean-field simulations plotted as a function of β = (k B T ) −1 . The dashed line represents the temperature dependence as given by the Callen-Callen relation defined by Eq.
(3). The solid line represents the modified scaling to enhance the rate of anisotropy increase with decreasing temperatures (increasing β), defined by Eq. (4).