Control of stripe, skyrmion and skyrmionium formation in the 2D magnet Fe3−xGeTe2 by varying composition

Two-dimensional (2D) van der Waals (vdW) magnets have recently emerged as novel skyrmion hosts. This discovery has opened a new material platform for tuning the properties of topological spin textures, such as by exploiting proximity effects induced by stacking of 2D materials into heterostructures, or by directly manipulating the structural composition of the host material. Previous works have considered the effect of varied composition in the bulk crystals of the vdW magnet Fe3−xGeTe2 , but so far the effects on the hosted spin textures have not been thoroughly investigated. In this work, real-space x-ray microscopy is utilized to image magnetic stripe domain, skyrmion and composite skyrmion states in exfoliated flakes of Fe3−xGeTe2 with varying Fe deficiency x. In combination with supporting mean-field and micromagnetic simulations, the significant alterations in the magnetic phase diagrams of the flakes, and thus the stability of the observed spin textures, are revealed. These arise as a result of the varying temperature dependence of the fundamental magnetic properties, which are greater than can be explained by the removal of spins, and are consistent with previously reported changes in the electronic band structure via the Fe deficiency.


Introduction
Magnetic skyrmions and related spin textures have received significant interest as potential data storage elements in future spintronic devices due to their topological and transport properties [1,2].The formation of monochiral skyrmions is induced by the Dzyaloshinskii-Moriya interaction (DMI), which arises due to the broken inversion symmetry of the host system [3].The particular underlying symmetry class may stabilize different forms of skyrmions: Néeltype skyrmions in ferromagnet/heavy metal multilayers with interfacial DMI [4][5][6], or in bulk polar magnets [7] (symmetry C nv ); Bloch-type skyrmions in bulk chiral magnets such as the B20 materials (symmetries T and O) [8,9] or antiskyrmions in bulk systems such as the Heusler alloys [10] or the recently discovered Fe 1.9 Ni 0.9 Pd 0.2 P [11] (symmetries D 2d and S 4 , respectively).
The properties of these skyrmions are further defined by the interplay of the DMI with other magnetic parameters, including the exchange and dipolar interactions, magnetocrystalline anisotropy, as well as extrinsic factors such as the shape of the sample, and the presence of structural defects and disorder.The balance of these contributions determines the size of the skyrmions [12], the extent of their stability in temperature and applied magnetic field [13], the exhibited type of skyrmion lattice ordering [14], and their transformation into other non-trivial spin textures [15].Moreover, skyrmion dynamics such as the lifetime of metastable states [16], response to microwave resonance excitations [17], and currentinduced mobility [18] may all be affected.Methods of controlling these fundamental material parameters are therefore highly valuable for engineering desired skyrmion properties.In bulk magnet skyrmion hosts, such methods have included modifying the stoichiometry of the basic material, which, for example in the Co 10−x Zn 10−y Mn x+y alloys [19], resulted in the alteration of both T C and the skyrmion size, as well the discovery of different skyrmion lattice orderings (triangular and square) [20][21][22].Alternatively, control has been achieved by partial or complete substitution of particular elements in the material composition.Well-studied examples include Mn 1−x Fe x Ge [12,23], and Fe 1−x Co x Si [24,25], where modulation of both T C and the skyrmion size was achieved by adjustment of the DMI and exchange interactions, or (Cu 1−x Zn x ) 2 OSeO 2 [26], where the increased disorder from the substitution resulted in a large enhancement of the metastable skyrmion lifetime via pinning effects [27].In thin-film systems, refinement of the skyrmion properties has largely been realized by varying the composition, thickness, and number of repetitions of each material layer in the heterostructure, enabling control of skyrmion density, size, and stability [28][29][30].
Control of the material parameters of FGT, and thus manipulation of the skyrmion properties, has now begun to be investigated.In particular, exploitation of the proximity effects when stacking FGT with other 2D materials has resulted in initial prototypes such as spin-valve and spin-orbit-torque devices [55][56][57].Moreover, the stability of skyrmions was found to increase when FGT flakes were combined with additional layers such as WTe 2 or Co/Pd multilayers [58,59].More fundamentally, compositional effects of FGT, such as altered stoichiometry [60][61][62] or chemical substitution [63][64][65], have previously been studied by bulk measurements.In particular, a previous work has highlighted the possibility to greatly modulate the magnetocrystalline anisotropy of FGT via hole doping, achieved by changing the Fe composition of the underlying crystal [62].Similar studies have been carried out with the related compound Fe 5 GeTe 2 [66][67][68].However, the direct effects of these alterations on the topological spin textures in FGT flake samples have yet to be observed.Such investigations may further clarify the origin of skyrmion formation in FGT, as well as provide a pathway for tailoring the functionality of proposed FGT-based devices [69,70].
In this work, we study both bulk and exfoliated flake samples of Fe 3−x GeTe 2 , with Fe deficiency x between 0.03 and 0.37.Using a combination of magnetometry and real-space x-ray microscopy, we investigate the effect that this altered composition has on the formation and stability of all observed spin textures, including stripe domain, skyrmion and skyrmionium states.Supporting mean-field and micromagnetic simulations confirm the vital role that the temperature varying uniaxial anisotropy plays in altering the magnetic phase diagrams of the material.The results provide a foundation for future works seeking to achieve fine control over topological spin textures in Fe 3−x GeTe 2 and other 2D magnets.

Fe 3−x GeTe 2 flake samples
We selected four bulk single crystals of Fe 3−x GeTe 2 to investigate [61], with measured Fe deficiencies of x = 0.03, 0.10, 0.27 and 0.37 (see methods).Results for magnetization M versus temperature T measurements performed with a magnetic field of 30 mT applied along the c axis of each bulk crystal are shown in figure 1(a) (see methods).In agreement with previous studies [60,61], the Curie temperature T C is reduced with higher Fe deficiency.We acquired values of 215, 202, 187 and 161 K for the x = 0.03, 0.10, 0.27 and 0.37 compositions, respectively, by finding the temperature point exhibiting largest dM/dT.Further M versus T measurements are plotted in figure 1(b), this time with a magnetic field of 3 mT applied perpendicular to the c axis of each crystal.In addition to the typical sharp increase of M at T C , there is a further step in each curve at lower temperature, labelled as T * .Previously, features observed at similar temperature ranges in FGT have been argued to be evidence for both antiferromagnetic ordering [71,72], or the onset of heavy fermion behaviour [73].Given that we observed that this additional step feature is most prominent at very low magnetic fields applied in the ab plane, and it is not observed for fields above 100 mT (see further M versus T data in supplementary data figure S1), we suggest that this feature may be related to the interlayer ordering between the individual vdW layers occurring at a lower temperature than the intralayer ordering at T C .
The measured values of T C and T * for each sample are plotted as a function of x in figure 1(c), alongside a measurement of the saturation magnetization M S extracted from M versus H measurements performed at 5 K.Such a reduction in M S can be expected by the removal of magnetic Fe atoms from the crystal lattice.As shown in figure 1(d), the crystal structure of FGT is composed of layers separated by a vdW gap, and exhibits two distinct Fe atomic sites: the Fe I site located closer to the outer Te ions; and the Fe II site, located in the central layer of the slab along with the Ge ions.Previous studies of Fe deficient samples have found that there is typically a lower occupancy of the Fe II site, while the occupancy of the Fe I site remains complete.In addition, higher Fe deficiency typically increases and decreases the c and a lattice parameters, respectively [60,61].
From these bulk crystals, we prepared three exfoliated flakes of FGT with x = 0.03, 0.27 and 0.37 compositions on Si 3 N 4 membranes for scanning transmission x-ray microscopy (STXM) measurements.The FGT flakes were each capped with flakes of protective hexaboron nitride (hBN) to prevent further oxidation under ambient conditions (methods).Optical micrographs of the prepared samples are shown in figures 1(e)-(g), highlighting the regions of interest (ROI) where subsequent STXM measurements were performed.From the optical contrast, we selected flakes with thicknesses as similar as possible (around 80 nm), to allow for a meaningful comparison of the spin texture formation in each composition (atomic force microscopy data shown in supplementary data figure S2).By tuning the incoming x-rays to the Fe L 3 absorption edge and exploiting the x-ray magnetic circular dichroism (XMCD) effect, STXM imaging yields a magnetic contrast signal proportional to the out-of-plane magnetization averaged through the thickness of the flake sample, m z (schematic illustration and x-ray absorption spectra shown in supplementary data figure S3).The images in figures 1(h)-(j) are x-ray micrographs of the ROI, revealing the formation of a dense array of skyrmions after field-cooling (FC) each flake sample from above T C to the indicated temperatures under an applied out-of-plane magnetic field of 20 mT.Such a method of skyrmion formation has been seen for the x = 0 FGT composition before [37,38,42], and here we demonstrate that similar skyrmion states can be readily realized across the full range of investigated FGT compositions.Previous observations have demonstrated FGT flakes within this thickness range exhibit Néel-type domain walls, indicating a source of interfacial-like DMI, although the origin of this effect remains under debate [38][39][40]42].

Magnetic phase diagrams
To investigate the formation of skyrmion and stripe domains across the composition series, we performed extensive STXM measurements on all three flake samples following a field-sweep (FS) procedure: at a range of temperatures, each sample was initialized in the saturated, monodomain state at −250 mT, and imaging was performed as the applied field was increased stepwise up to +250 mT.A selection of x-ray micrographs acquired following this procedure for different temperatures is shown in figures 2(a)-(o) (further images of each sample presented in supplementary data figures S4-S6).Alongside each row, we have labelled the fraction of T C of each sample at which the measurement was performed, to facilitate comparison between each composition.The T C of each flake sample was determined as the temperature at which real-space magnetic contrast could no longer be observed, and the values were found to be slightly lower than in the bulk: 207, 180 and 146 K for the x = 0.03, 0.27 and 0.37 compositions, respectively.The reductions of T C are consistent with the decreases seen in thinner flakes of 2D magnets [47], but we also cannot discount a temperature offset of a few Kelvin in the thermocouple measurement.The results for the x = 0.03 flake acquired at 203 K, or 0.98 T C are shown figure 2(a).The images reveal the formation of a dense disordered array of skyrmions (labelled Sk) for both positive and negative applied OOP fields, as well as stripe domains (labelled SD) at 0 mT.At decreasing temperatures, shown in figures 2(b) and (c), the characteristic size of the stripe domains increases, and skyrmion formation is no longer observed.Finally, at 150 K (0.72 T C ) and below, only uniform switching between the positive and negative monodomain states (labelled MD±) is observed.
The high temperature results for the two flakes with greater Fe deficiency, x = 0.27 and x = 0.37, shown in figures 2(f), (g) and (k), (l), reveal a similar behaviour, with dense skyrmion formation only observed close to T C , and the characteristic stripe domain size increasing with decreasing temperature.However, at lower temperatures, there is a significant difference in the crossover to monodomain switching behaviour.In comparison to the x = 0.03 flake, this occurs at a comparatively lower temperature in the x = 0.27 sample (0.61 T C ), revealed in figures 2(h)-(j).Furthermore, in the x = 0.37 sample, monodomain switching behaviour was not observed down to the base temperature of the STXM instrument at 30 K-instead, we observed the formation of stripe domains across the full investigated temperature range, as shown in figures 2(m)-(o).
The overall behaviour is visible from the magnetic phase diagrams of each flake sample presented in figures 3(a)-(c), which plot the observation of each magnetic state as a function of the applied magnetic field at each temperature when following the FS procedure.In cases where both stripe and skyrmion states coexisted, these were included in the skyrmion regions for clarity.The temperature is plotted on both an absolute scale, and as a fraction of T C .The phase diagram of the x = 0.03 sample, with the crossover to the monodomain switching behaviour, is similar to those previously reported for FGT flakes with x close to 0 [38].However, the results of the x = 0.37 flake, where real-space spin textures were observed across the full temperature range, is more reminiscent of typical dipolar-stabilized skyrmion bubble systems [74], as well as the multilayer skyrmion hosts [75].
To explore the stability range of the skyrmion state in each flake sample, we performed further xray imaging measurements following the FC procedure: the sample was initialized above T C , cooled down to the target temperature under an applied field of 20 mT, and images were acquired as a function of both increasing and decreasing applied field (with two separate cooling procedures).In all three samples, the FC procedure resulted in the formation of disordered arrays of skyrmions at low temperature, which were effectively quenched from the high temperature skyrmion pockets close to T C , as shown in figures 1(h)-(j).The summary phase diagrams in figures 3(d)-(f) show the field and temperature extent of the resulting metastable skyrmion state in each sample, which is greatly increased in comparison to the equilibrium skyrmion phase revealed in the FS measurements.This large extent of the FC skyrmion phase was observed for all three compositions, with the main difference being a slight change in the field range of the skyrmion stability-for the x = 0.37 sample, the stability of the skyrmions was both higher in the positive field direction, and lower in the negative field direction, in comparison to the x = 0.03 sample.
Note that these results show the phase diagram for a FS procedure from negative to positive magnetic field, and for a FC procedure with a positive applied field.In both cases, with a change in sign of the applied field, each phase diagram would be mirrored in the x axis.In all three samples, for increasing applied field, the skyrmions reduced in size and number, while for decreasing field their size dramatically increased below 0 mT (corresponding microscopy images shown in supplementary data figure S7).As discussed in previous works, this highly fielddependent skyrmion size is a strong indication of the predominantly dipolar-stabilized character of the observed skyrmion states, rather than DMI-stabilised [38,76].These results demonstrate that large differences in the low temperature behaviour of the sample are accessible through a modulation of the Fe content of FGT samples, while formation of skyrmions is maintained across the full investigated compositional range.Interestingly, we did not notice any different behaviour in the flake samples associated with T * , suggesting this may only be relevant for the bulk crystal samples.

Interplay of magnetic energy terms
We sought to perform simulations in order to reproduce the experimental results, and investigate the origin of the observed differences in the magnetic phase diagrams.In order to acquire reasonable model parameters, we performed extensive magnetometry measurements of the bulk single crystals.Figures 4(a)-(c) plots M of each FGT composition measured at 5 K as a function of H with the field applied both parallel and perpendicular to the c axis.Immediately, the plots reveal a significant change in the saturation field between the three compositions when the field is applied in the ab plane.This indicates a decrease in the uniaxial anisotropy, K U , responsible for the easy axis along the c axis, as a function of increasing Fe deficiency, as has been observed previously [60].
By measuring further M versus H loops at a range of temperatures, we extracted both M S and, using an integration of the hysteresis loops, estimated values of K U for each sample (further data shown in supplementary data figure S8, see methods).The calculated parameters for each composition are plotted as a function of temperature in figures 4(d) and (e).Note that in plotting K U , we estimated and subtracted the shape anisotropy contribution of each plate-like bulk sample from the measured effective anisotropy K eff (details shown in supplementary data figure S9).The results show both the expected decrease in M S as well as the strong change in K U with increasing x.The inset of figure 4(e) plots the extracted value of the anisotropy field H K = 2K u /M S , as a function of T for each sample, demonstrating that the change in K U between samples cannot be accounted for by the change in M S .
This implies that the magnetocrystalline anisotropy contribution is tuned by a factor of 4 for a change in x of ~0.35, in line with previous estimates [60].Notably, this change in anisotropy with x is larger than the change in M S , indicating a significant change in the magnetic behaviour of the system.Based on angle-resolved photoelectron spectroscopy and density functional theory calculations, this has been argued to be linked to fundamental changes in the electronic band structure of FGT [62].Due to the resulting hole doping at higher values of x, the majority of bands that are strongly affected by spin-orbit coupling are shifted away from the Fermi energy with increasing Fe deficiency, thus reducing the resulting out-of-plane magnetocrystalline anisotropy [62].
Considering this strong change in the uniaxial anisotropy with the composition of the samples, we modelled the experimental systems using a temperature-dependent computational model based on the mean-field approximation of a classical spin Hamiltonian with exchange, DMI, uniaxial anisotropy and dipolar interaction terms (see methods), similar to our previous work [38,77].We noticed that the M S and K U values could be collapsed upon one another via a linear scaling of both axes, indicating a homogeneity of the M S (T) and K U (T) functions across the composition series (see supplementary data figure S10).Based on a dimensional analysis of the model Hamiltonian (see supporting note A and B), it is possible to reduce the complexity of the simulation to only rescaling a few model parameters while keeping the temperature dependence of all energetic terms the same.For example, to obtain the exchange parameter J E of a system with Fe deficiency x 2 from an initial set of parameters in the system with x 1 , one can utilize ) .
Thus, starting from a set of model parameters for the x = 0.03 system, we calculated parameters for each composition system, scaled by the experimentally measured M S (T) and K U (T) values shown in figure 4. A more detailed and thorough description of these scaling arguments is included in the methods and the supporting material.We performed simulations of each system, following a field-sweep procedure starting at negative applied fields, with a summary of the results shown in figure 5 (further data shown in supplementary data figures S11-S13).Simulations 1, 2 and 3 correspond to parameters selected to model the x = 0.03, 0.27 and 0.37 experimental systems, respectively.The visualizations of the simulations in figures 5(a)-(c) show a good agreement with the experimental behaviour in figure 2. The formation of skyrmions at high temperatures is reproduced, while the observed crossover from stripe domain formation to uniform magnetization switching is seen in the case of simulations 1 and 2. The results are better visible in the simulated magnetic phase diagrams in figures 5(d)-(f).Comparison to figures 3(a)-(c) shows a reasonable agreement to the experimentally determined phase diagrams.In particular, specific features such as the applied field asymmetry of both the high temperature skyrmion pockets and the stripe domain formation (due to the field asymmetry) are reproduced.The main features that are poorly replicated are the switching fields at low temperatures, which might be explained by the presence of thermal fluctuations or defects allowing easier switching in the experimental system.Note that we did not attempt reproduction of the experimental phase diagram following the FC procedure, but expect that similar agreement would be achieved.The results demonstrate that the balance of the dipolar energy from M S and the anisotropy from K U are responsible for the magnetic phase diagrams of these FGT flakes.

Composite skyrmion formation
In a previous work, we reported the observation of composite skyrmions in exfoliated FGT flakes [51], with topological charge |Q| ̸ = 1 (the skyrmions shown so far have Q = −1).Their formation was realized by the seeding of loop-like states within the stripe domain state when following a zero-field cooling (ZFC) procedure, from which skyrmionium (Q = 0) [78][79][80][81], skyrmion bag (Q > 1) [52,53,[82][83][84] and skyrmion sack (Q < −1) [53,54,85] states emerged upon increasing an out-of-plane magnetic field.To investigate the formation of these magnetic states in the Fe 3−x GeTe 2 flakes, we performed imaging following the ZFC procedure to a range of temperatures.As exemplified in figures 6(a)-(c), in all three flakes, loop-like seed states were observed within the stripe domains upon ZFC (further data shown in supplementary data figures S14-S16).In many cases, isolated skyrmioniums emerged upon increasing the applied out-of-plane magnetic field, and occasionally more complex composite skyrmion states with |Q| > 1 were created.
The number of skyrmionium states (N SkM ) observed in each FGT flake is plotted as a function of temperature in figure 6(d).For all three samples, there appears to be a higher probability to form skyrmionium states at higher temperatures, close to T C .For the x = 0.03 and 0.27 compositions, the probability to form skyrmionium states decreases significantly for lower temperatures.On the other hand, in the x = 0.37 FGT flake, the number of observed skyrmionium states increases for lower temperatures.Figures 6(e)-(g) plots the observation of skyrmionium states as a function of both applied field and temperature, for the x = 0.03, 0.27 and 0.36 samples, respectively.The spot sizes indicate the average size of skyrmionium states observed at each temperature and field point.These phase diagrams highlight a dramatic change in skyrmionium formation characteristics across the composition series.In general, once formed, an individual skyrmionium state will decrease in size for increasing applied magnetic field [51].However, as revealed by figures 6(e)-(g), most often the skyrmionium states with larger size emerge at, and survive to, higher applied magnetic fields (see supplementary data figure S17 for more details on skyrmionium sizes).Note that due to the limited time for central facility-based measurements, we have limited statistics for these measurements.Due to the random nature of the seeded formation mechanism, more thorough analysis would require repeated imaging runs to improve the averaging process.
Nevertheless, it is interesting to speculate on the observed differences between FGT compositions.For example, the increased chance to form composite skyrmions at lower temperatures in the x = 0.36 flake may be related to the lack of low temperature monodomain switching behaviour observed in the FS measurements presented in previous figures.One might expect that the observed increased stability of stripe domains, likely due to the decreased uniaxial anisotropy lowering the cost to form magnetic domain walls, would be extended to the skyrmionium states.Alternatively, we can expect that nonstoichiometric samples exhibit increased disorder and defect density within the crystal structure, which may act as pinning sites and in turn stabilise the observed magnetic textures.As further evidence for this possible pinning effect, we note that the composite skyrmion states observed in the x = 0.36 flake often exhibited irregular shapes, while in the x = 0.03 sample they were often more close to circular, as one would expect from an unpinned domain structure.

Micromagnetic simulations
To investigate the stability of composite skyrmion states for varying magnetic properties, we performed multiple micromagnetic simulations with a range of saturation magnetization, M S (100 to 400 kA m −1 ) and out-of-plane anisotropy field H K (100 to 400 mT), including terms for the demagnetizing field and a small interfacial-like DMI, with D = 0.14 mJ m −2 , which reproduces the Néel domain walls observed experimentally in similar FGT flakes.We initialized the simulation while including temperature fluctuations, and then gradually reduced these while minimizing the energy of the system in an attempt to model the experimental ZFC procedure (see methods, supplementary data figure S18).Once at zero fluctuations, we increased the applied field in a stepwise fashion.A summary of the results is shown in figures 7(a)-(p), where we plot the field state at which the system exhibited the most composite skyrmion spin textures.It is evident that a variety of loop-like composite skyrmion states are realized for a wide range of micromagnetic parameters as shown by figures 7(q) and (r).
First we shall consider the results for M S ⩾ 200 kA m −1 in more detail, shown in figures 7(e)-(p).The domain sizes in these sets of results exhibit typical behaviour for a magnetic film with out-of-plane anisotropy, where the domain formation is dominated by the dipolar interaction [86].For increasing M S , the typical domain sizes become smaller, while for increasing H K the domain wall energy cost increases, leading to larger domains [87].For smaller H K values, we also observed the possibility to stabilise Bloch point-like structures within the magnetic domain walls, which have been termed chiral kinks [54], and further modify the topological charge of the After ZFC, the temperature is effectively 0 K. (q)-(t), Zooms into regions containing a skyrmion and skyrmionium (q), skyrmion bag (r), and objects with Bloch points (or chiral kinks) embedded within their domain walls (s), (t).The colour map indicates both the out-of-plane magnetization component, mz (black and white), and the in-plane components, mx and my (rainbow colour wheel).
object by introducing additional winding, as shown in figures 7(s) and (t).We note that in simulations performed with a slower cooling rate the number of such Bloch points were reduced (see supplementary data figure S19), tentatively indicating that fast cooling rates (exceeding 10 K min −1 ) may be required to realize them in experimental samples.The trend in domain size is reversed when considering the results from the M S = 100 kA m −1 simulations, where the domain sizes are considerably smaller.We argue that this is caused by the onset of negative domain wall energy owing to the dominance of the DMI contribution for sufficiently small M S [76] (see supplementary data figure S19 for simulation with zero DMI).In this DMI-dominated regime [88], the diversity and density of composite skyrmion states is increased.Due to the large M S of the experimentally investigated FGT flakes, it is likely that the investigated samples are in the dipolar-dominated regime.
In this work, our x-ray microscopy method is insensitive to the in-plane magnetisation, we cannot determine nature of the magnetic domain walls.However, in previous Lorentz transmission electron microscopy measurements of the flakes from the x = 0.03 and x = 0.37 crystals, we saw magnetic contrast consistent with Néel-type domain walls [38,51].One could approach this from the simulation angle, utilziing micromagnetics to estimate the size of the DMI.However, one signficant challenge is that typical micromagnetic simulation methods utilize an isotropic exchange interaction.In layered magnets such as FGT, we can expect the intralayer and interlayer exchange couplings are different.Therefore, further work is needed to develop truly quantitative micromagnetic simulations of the spin textures in 2D magnets.Nevertheless, our simulations show that in both dipolar and DMI dominated systems, the general seeding of composite skyrmions is similar, demonstrating that the zero-field cooling seeding mechanism may be generally applicable to a wide range of materials, and this should be a fruitful avenue for future research.

Conclusion
We have shown comprehensive x-ray microscopy results of the magnetic spin textures hosted in exfoliated flakes of Fe 3−x GeTe 2 with Fe deficiency between x = 0.03 and 0.37, observing the formation of stripe domain, skyrmion and composite skyrmion states.By mapping the magnetic phase diagrams, we identified significant differences between the compositions.Specifically, close to x = 0.00, FGT exhibits a monodomain switching behaviour at low temperatures, while in the x = 0.37 sample, the formation of stripe domains was maintained to the base temperature of our instrument (30 K).Supporting mean-field simulations indicate that this behaviour is due to the decrease of the uniaxial anisotropy, which could be realized by an alteration of the underlying band structure with higher Fe deficiency.
Moreover, we observed that skyrmionium and other composite skyrmion states were realizable across all compositions, emerging from seeded states in the stripe domain after a zero field cooling procedure.The density and temperature-field extent of the formation was much greater in the x = 0.37 sample, although it is unclear whether this is due to the change in uniaxial anisotropy, or due to pinning effects.Accompanying micromagnetic simulations demonstrated that skyrmionium states should be realizable across a wide range of anisotropy and saturation magnetization, in both the DMI and dipolar dominated regimes.The results provide a foundation for future compositional studies of spin textures in FGT flakes, and emphasize that composite skyrmion states may exist in a wide range of magnetic systems.The relationship between changes in the electronic structure and the properties of the hosted spin textures indicates that one should expect interesting results from other methods capable of manipulating the band structure, such as electrostatic or ionic gating.In the future, the combination of FGT flakes with additional 2D materials will allow for the tuning of the relevant magnetic energy terms by proximity effects, enabling greater control of the magnetic spin textures.

Sample preparation
Bulk crystals of Fe 3−x GeTe 2 were grown via the chemical vapour transport technique.Stoichiometric quantities of Fe (STREM Chemicals, Inc. 99.99%), Ge (Acros Organics, 99.999%), and Te (Alfa Aesar, 99.99%) powders, along with 5 mg/cm 3 of iodine as the transport agent, were sealed within evacuated quartz ampules.The single crystal growth was performed using a two-zone furnace by holding the source and sink parts of the tube at 750 • C and 675 • C, respectively, for 2 weeks, before cooling to room temperature.The growth procedure resulted in the formation of silvery metallic platelet single crystals with areas of ~2×2 mm 2 .Energy dispersive x-ray spectroscopy measurements using a scanning electron x-ray microscope determined the composition of the four chosen single crystal samples to be x = 0.03, 0.10, 0.27, and 0.37, as reported previously [61].
To prepare the flake samples on the Si 3 N 4 membranes for the x-ray transmission microscopy measurements, we utilized an all-dry viscoelastic transfer method.Flake samples were prepared from crystals with compositions x = 0.03, 0.27, and 0.37.The bulk FGT flakes were mechanically cleaved and exfoliated onto a PDMS stamp using Nitto scotch tape.By inspecting their contrast under an optical microscope, flakes with a thickness of around 80 nm were selected, and stamped onto a Si 3 N 4 membrane.Subsequently, a similar exfoliation and stamping transfer procedure was performed to cap the FGT flakes with larger hBN flakes with thickness ~15 nm.All fabrication steps were performed under ambient conditions, such that each side of the FGT flakes was exposed to the atmosphere for ~20 min.Our previous measurements indicate that this results in an oxidized FGT layer with an approximate thickness of 6 nm on both sides of the FGT flake [38].The thicknesses of the resulting FGT/hBN heterostructures was measured using atomic force microscopy measurements with a Bruker Dimension Icon, as shown in supplementary data figure S2.

Magnetometry measurements
Magnetometry measurements were carried out using a Quantum Design MPMS3 vibrating sample magnetometer.Bulk crystals of each FGT composition were fixed to a quartz rod sample holder using GE varnish.Measurements were performed for each sample with the field aligned with both the c crystal axis (H ∥ c) and in the ab plane (H ⊥ c).The sample temperature and applied magnetic field value were controlled by the built-in helium cryostat.The raw data, which yields values for the magnetic moment, was converted to SI units of magnetization using the measured the measured mass of each crystal and assuming a density of FGT of 7.3 g cm −1 for all samples.Measurements of the sample magnetization were acquired as a function of both sample temperature and applied magnetic field.
An approximate value of the uniaxial anisotropy at each temperature was estimated from the magnetization versus field data measured along both crystal orientations with the following process.A linear trend was acquired from the highest field (in the saturated state) data points of the measured hysteresis loop at each temperature, and the resulting fit was taken as a diamagnetic and paramagnetic background subtraction.This allowed a determination of the saturation magnetization, M S , at each temperature.From this data, the effective anisotropy, K eff of the bulk sample was estimated from the work done, W, to reach M S , along each crystal orientation Ĥ, (2) Thus, by integrating over the magnetization versus field curves for each orientation, the measured effective anisotropy could be calculated, K eff = W H∥c − W H⊥c .Part of this measured anisotropy will be the result of the shape anisotropy, which is particularly strong in the plate-like FGT bulk samples.This shape anisotropy K sh , or demagnetization effect, can be calculated from M S and the demagnetization factor N d = 0.93 estimated from the dimensions of the platelike samples, by The negative sign highlights the fact that the shape anisotropy for a plate-like sample seeks to align spins in the plane.Thus, the magnetocrystalline uniaxial anisotropy K U can be approximated by K U = K eff − K sh .Finally, the anisotropy field H K can be extracted as These calculations were performed for all measured bulk samples as shown in figure 4 in the main text.

Scanning transmission x-ray microscopy
STXM measurements were performed with the MAXYMUS endstation at the UE46 beamline at the BESSY II electron storage ring operated by the Helmholtz-Zentrum Berlin für Materialien und Energie.Each Si 3 N 4 membrane holding the FGT/hBN stacked flake samples was fixed to the sample holder using GE varnish.The sample holder was placed into the vacuum chamber of the instrument, with cooling achieved by a liquid He cryostat.Applied out-of-plane magnetic field control was achieved by varying the arrangement of four permanent magnets within the sample chamber.The incoming x-ray beam was focused to a small spot size using a Fresnel zone plate with central stop, and the remaining zeroth order light was isolated by positioning the order selection aperture.To acquire an x-ray micrograph, the sample was rastered through the focused beam using a piezoelectric motor stage, with the x-ray transmission measured pixel by pixel with an avalanche photodiode.Magnetic contrast was achieved by exploiting the XMCD effect, where magnetization along the x-ray propagation direction, m z alters the absorption of circularly polarized x-rays (see supplementary data figure S3).Thus, an image of the m z component of the present magnetic spin texture can be acquired.An XMCD signal can be found at the resonant L 3 and L 2 absorption edges of the magnetic element within the sample; in this case we selected the Fe L 3 edge with a nominal x-ray energy of 707.5 eV.All presented images of the magnetic domain structure were recorded using a single circular x-ray polarization (positive).In some cases, such as in figure 4 of the main text, the magnetic contrast was isolated from the structural contrast by subtracting an image of the sample saturated at an applied field of +250 mT, leaving only contrast associated with the magnetic domains.

Mean-field simulations
To perform temperature-dependent simulations we developed a mean-field model based on the standard classical spin Hamiltonian [77] with spins distributed on a 2D hexagonal lattice of 30 × 30 spins (approx.150 nm × 150 nm) with periodic boundary conditions in the plane.The mean-field energy reads: where µm represents the vector moment of meanfield spins, with µ being the zero-temperature moment magnitude and m i the temperaturedependent magnetic moment vector of length normalized to take values between 0 ⩽ |m i | ⩽ 1.The individual energy terms starting from the left correspond to exchange, DMI, uniaxial anisotropy, Zeeman energy, and finally the dipolar energy.The symbol ⟨•⟩ appearing in the first two sums in equation ( 5) implies the nearest neighbour coupling, and mi in the anisotropy term is a vector of unit length.The interaction constants J E , J D , and J DP are in the units of Joule ×µ −2 .The exchange and anisotropy energy constants J E and J K are temperature-dependent (see supplementary data for details).The temperaturedependence of J D is weak and will be ignored.Thus the DMI energy term depends on temperature only through the moments m i .The DMI vector d ij is oriented in the lattice plane and perpendicular to the line connecting two neighbouring spins i and j.The anisotropy unit vector n is oriented along the ẑ-axis perpendicular to the lattice plane.The spin moment m i has temperature dependent magnitude following the expression: where L(x) = coth x − x −1 is the Langevin function, k B is the Boltzmann constant, and T the temperature.The vector B e i represents the effective field acting on the spin m i and |B e i | is the magnitude of B e i .The expression for the effective field can be obtained from equation ( 5) by calculating the variational derivative: Thus according to equations ( 6) and ( 7), the meanfield spin m i depends on the exchange and DMI energy couplings, material anisotropy, external and dipolar field, and also on the temperature T. The magnetization structures for a given set of parameters are evaluated by minimizing equation ( 5) during the field or temperature evolution starting from uniquely defined initial states, such as fully saturated uniform state in case of the field sweep computations, or random vector distribution in the temperature cooling computations starting from the paramagnetic state [77].The parameters required to generate the results in this work and the temperature dependencies of J E and J K are discussed in detail as part of the supplementary data.The supplementary material identifies the relationship between the model constants J E , J D , J K and their micromagnetic equivalents A (J m −1 ), D (J m −2 ) and K (J m −3 ).

Micromagnetic simulations
Micromagnetic simulations were performed using the MicroMagnum framework following the Landau-Lifschitz-Gilbert (LLG) equation.We added custom extensions for including a DMI tensor and temperature fluctuations via the effective field approach [89].The simulated system had dimensions of 1 × 1 × 0.02 µm 3 , with cell volumes of 1 × 1 × 20 nm 3 .The simulations were carried out for a range of M S and H K values, with the exchange and DMI constants set to A = 0.7 pJ m −1 and D = 0.14 mJ m −2 , respectively.The thermal fluctuations utilized an Euler-type solver with a spatially and temporally decorrelated thermal field and a fixed 10 fs time stepping [90].To model the experimentally realized ZFC procedure, the system was initialized with random magnetization at the micromagnetic T C (the temperature setting calibrated to where all magnetization becomes fully random).The temperature was subsequently reduced to 0.3T C within the next 3.5 ns and then held at that temperature for another 2.5 ns.The entire ZFC process thus took in total 7 ns.In addition to the thermal fluctuations, the saturation magnetization was scaled with the decreasing temperature following the experimentally observed behaviour (and thus also the anisotropy following K U = µ0 2 M S H K + µ0 2 M 2 S with the effective anisotropy field H K (in this definition H K = 0 mT corresponds to an out-of-plane to in-plane transition).Once at low temperatures, the thermal fluctuations were switched off, and the applied magnetic field was subsequently increased in steps of 10 mT, with the system relaxed at each step (initial state after cooling shown in supplementary data figure S18).No initial magnetization pattern or manipulation was used, such that the results present a close representation of a zero field cooling process through T C .
In these simulations, we utilized temperature fluctuations to closely model the introduction of randomness observed in the stripe domain state following the experimental ZFC process.However, the included thermal fluctuations are only an approximate model, due to the inability of the LLG to vary the magnitude of the magnetic moment, and because of the high temperature cutoff of spin waves from the finite cell size.Furthermore, although we present results obtained with a range of micromagnetic parameters, we note that the utilized parameters are not equal to the experimentally measured ones for each FGT flake.However, they do correspond to effective values that reproduce the relative sizes and morphology of the spin textures in the available simulation area.Simulating the full size of the FGT flake of several tens of µm would not be feasible.

Figure 1 .
Figure 1.Fe 3−x GeTe2 bulk and flake sample characterization.(a) Magnetization M measured as a function of temperature for four FGT bulk crystals with differing composition (denoted by the Fe deficiency value x), with a field of 30 mT applied along the c crystal axis.The Curie temperature TC of each sample is indicated by the vertical dashed lines.(b) Magnetization measured as a function of temperature for each composition with a field of 3 mT applied perpendicular to the c crystal axis.The TC and second transition temperature T * of each sample are labelled.(c) The extracted values of TC (purple), T * (orange) and the saturation magnetization MS (determined at 5 K) of each composition plotted as a function of x.Error bars on the T values indicate the error due to effects of broadening of the M curve.(d) Crystal structure of Fe 3−x GeTe2 (FGT) from the side and top viewpoints.The extent of the unit cell and the identity of the Fe I, Fe II, Ge and Te atom sites are labelled.The orientation of the crystal axes are indicated.(e)-(g) Optical microscopy images of the hBN-capped exfoliated FGT flakes on the Si3N4 membranes, for the x = 0.03, 0.27 and 0.37 compositions respectively.The region of interest (ROI) indicates the location of the scanning transmission x-ray microscopy (STXM) measurements in subsequent figures.(h)-(j) STXM images of skyrmions at the ROI created by field-cooling (FC) each flake from above TC to the indicated temperatures, under an applied field of 20 mT.The colour map indicates the out-of-plane component of the magnetization, mz (arbitrary units).

Figure 2 .
Figure 2. Scanning transmission microscopy measurements following the field-sweep procedure.(a)-(o) X-ray micrographs of the Fe 3−x GeTe2 (FGT) flake samples measured as a function of temperature and applied magnetic field for the x = 0.03 ((a)-(e)), 0.27 ((f)-(j)) and 0.37 ((k)-(o)) flakes respectively.The images were taken as a function of increasing out-of-plane applied magnetic field, starting in the saturated state at −250, as indicated by the orange arrows.The temperature as a fraction of the Curie temperature TC is labelled.The images were taken at the regions of interest as shown in figure 1.The colour map indicates the out-of-plane component of the magnetization, mz (arbitrary units).

Figure 3 .
Figure 3. Composition dependent magnetic phase diagrams.(a)-(f) Magnetic phase diagrams determined by x-ray microscopy of the three Fe 3−x GeTe2 (FGT) flakes as a function of temperature and applied field.Results following the field sweep (FS) procedure (a)-(c) and the field-cooling (FC) procedure (d)-(f) are shown for each composition x = 0.03, 0.27 and 0.37, as labelled.Red arrows in a and d indicate the measurement path of each field-temperature protocol.The extent of the skyrmion (Sk) stripe domain (SD) and uniformly magnetized monodomain (MD±) states is shown by the coloured red, blue and grey regions, respectively.Markers indicate the measured phase boundary points of the Sk (red squares) and SD (blue circle) states.The vertical dashed lines indicate the measured values of TC.For (a)-(c), additional data points exist at higher fields-here we have focused on the details in the low field regions.

Figure 4 .
Figure 4. Composition dependence of magnetization reversal and uniaxial anisotropy.(a)-(c) Measurements of the magnetization M versus applied field µ0H, at 5 K, of each Fe 3−x GeTe2 bulk crystal, with compositions x = 0.03, 0.27 and 0.37, respectively.Measurements were acquired with the magnetic field H applied both parallel (purple triangles) and perpendicular (orange circles) to the c crystalline axis.(d) Extracted values of the saturation magnetization MS as a function of temperature T for each FGT composition, x = 0.03, 0.10, 0.27 and 0.37.(e) Extracted values of the uniaxial anisotropy KU of each FGT sample as a function of T, taken from the H ∥ c data.We attempted subtraction of the shape anisotropy contribution from an estimate of the sample dimensions and the measured MS values (see methods).The inset plots the calculated value of the anisotropy field HK as a function of T for each composition.

Figure 5 .
Figure 5. Mean-field simulations modelling Fe 3−x GeTe2 flakes.Simulations 1-3 correspond to systems based on parameters for the x = 0.03, 0.27 and 0.37 experimental compositions, respectively.(a)-(c), Selected visualizations of simulated states acquired following a field sweep procedure starting from negative applied fields.The colour map indicates the out-of-plane component of the magnetization, mz (arbitrary units).(d)-(f) Simulated magnetic phase diagrams of the three systems.The extent of the skyrmion (Sk) stripe domain (SD) and uniformly magnetized monodomain (MD±) states is shown by the coloured red, blue and gray regions, respectively.Markers indicate the sampled positions in phase space.The simulated values of TC are indicated by the vertical lines.

Figure 6 .
Figure 6.Observation of skyrmionium states when zero-field cooling Fe 3−x GeTe2 flakes.(a)-(c) Exemplary x-ray micrographs of the x = 0.03, 0.27 and 0.37 FGT flakes, respectively.Images were recorded as a function of increasing applied magnetic field, after zero-field cooling (ZFC) to the target temperature from above TC.The formation of skyrmionium states is highlighted with the yellow dashed rings, and the edge of the flakes with green dashed lines.The colour map indicates the out-of-plane component of the magnetization mz (arbitrary units).In (b), the inset shows a zoomed view of one skyrmionium.(d) The total of number of skyrmioniums N SkM observed in each FGT flake during the ZFC field sweep performed at each temperature.The estimated TC of each FGT flake is indicated by the vertical dashed lines.(e)-(g) The observation of skyrmionium states as a function of temperature and applied magnetic field is shown for the x = 0.03, 0.27 and 0.37 FGT flakes, respectively.The average size of the skyrmionium states present at each temperature-field point is indicated by the size of the corresponding marker.

Figure 7 .
Figure 7. Micromagnetic simulations of composite skyrmion formation.(a)-(p) Selected field points of micromagnetic simulations modelling the zero-field cooling (ZFC) seeding procedure for a range of saturation magnetization MS and anisotropy field HK.After ZFC, the temperature is effectively 0 K. (q)-(t), Zooms into regions containing a skyrmion and skyrmionium (q), skyrmion bag (r), and objects with Bloch points (or chiral kinks) embedded within their domain walls (s), (t).The colour map indicates both the out-of-plane magnetization component, mz (black and white), and the in-plane components, mx and my (rainbow colour wheel).