Tunable gigahertz dynamics of low-temperature skyrmion lattice in a chiral magnet

Recently, it has been shown that the chiral magnetic insulator Cu$_2$OSeO$_3$ hosts skyrmions in two separated pockets in temperature and magnetic field phase space. It has also been shown that the predominant stabilization mechanism for the low-temperature skyrmion (LTS) phase is via the crystalline anisotropy, opposed to temperature fluctuations that stabilize the well-established high-temperature skyrmion (HTS) phase. Here, we report on a detailed study of LTS generation by field cycling, probed by GHz spin dynamics in Cu$_2$OSeO$_3$. LTSs are populated via a field cycling protocol with the static magnetic field applied parallel to the $\langle{100}\rangle$ crystalline direction of plate and cuboid-shaped bulk crystals. By analyzing temperature-dependent broadband spectroscopy data, clear evidence of low-temperature skyrmion excitations with clockwise (CW), counterclockwise (CCW), and breathing mode (BR) character at temperatures below $T$ = 40 K are shown. We find that the mode intensities can be tuned with the number of field-cycles below the saturation field. By tracking the resonance frequencies, we are able to map out the field-cycle-generated LTS phase diagram, from which we conclude that the LTS phase is distinctly separated from the high-temperature counterpart. We also study the mode hybridization between the dark CW and the BR modes as a function of temperature. By using two Cu$_2$OSeO$_3$ crystals with different shapes and therefore different demagnetization factors, together with numerical calculations, we unambiguously show that the magnetocrystalline anisotropy plays a central role for the mode hybridization.


INTRODUCTION
Chiral magnets with large spin-orbit interaction and significant Dzyaloshinskii-Moriya interaction (DMI) host different non-collinear spin configurations such as the helical, conical and skyrmion lattice configurations as well as the field polarized (ferrimagnetic) state in large applied magnetic fields. Many chiral magnets share a common magnetic phase diagram with a skyrmion lattice phase stabilized by thermal fluctuation located in a small pocket in temperature and magnetic field space in the vicinity of the Curie temperature [1][2][3][4][5][6][7]. Furthermore, their excitation dynamics in the GHz regime has been theoretically predicted and experimentally shown in prior studies [8][9][10][11].
The recent discovery of a novel low-temperature skyrmion (LTS) phase in various materials, including MnSi and Co 8 ZnMn 4 [12,13], has consequently sparked renewed interest in the particular chiral magnetic insulator Cu 2 OSeO 3 [14][15][16][17]. It has been shown that in this material, by applying appropriate magnetic field and temperature protocols, a low-temperature skyrmion phase can be nucleated, which is distinct from the high-temperature skyrmion (HTS) phase and which seems to be stabilized by a combination of DMI and the magnetocrystalline anisotropy [14,15,18]. As a consequence, it only appears if the externally applied magnetic field is oriented along a certain crystallographic direction. Furthermore, in contrast to the HTS phase, the nucleation process for the LTS is nontrivial: the metastable tilted conical phase seems to be a necessary prerequisite for the formation of the LTS lattice [14][15][16]. A recent study by Halder et al. [15] has shown that, to a leading order, a cubic magnetic anisotropy term is responsible for stabilizing the LTS lattice and that the characteristics of the phases can be described rather well using a Ginzburg-Landau formalism adapted for chiral magnets with cubic anisotropy [19,20]. The same formalism can be extended to compute the spin excitations in both LTS and HTS lattice phases. The description above suggests that one may influence the spatial configuration of the magnetization by controlling the strength of the magnetocrystalline anisotropy, e.g. by temperature.
From the above studies, it is reasonable to conclude that the stabilization mechanisms responsible for the LTS phase originate from the magnetic anisotropies when the magnetic field is applied parallel to 100 . Moreover, recently, the dynamic modes of this new LTS phase have also been thoroughly investigated using low-temperature broadband ferromagnetic resonance experiments, allowing in addition the identification of the low temperature elongated skyrmion phase [21].
In this paper, we address the dynamic modes in the LTS by applying broadband ferromagnetic resonance experiments as a function of temperature, T , and magnetic field, H. We report the influence of the number of field cycles on the intensity of the resonant spin excitations in the LTS phase and show that a large volume fraction of the investigated spin structure in bulk crystals can indeed be converted into the LTS phase by appropriate field-cycling [14,22]. Second, we observe the low-lying breathing and counterclockwise modes and find evidence of the clockwise gyrational modes at higher frequencies [8,9,11]. Using the dynamic modes as experimental evidence for the LTS, we map out the phase diagram as a function of H and T . Finally, we systematically trace the mode hybridization between a dark higher-order clockwise mode and the breathing mode [23] as a function of temperature for two differently shaped crystals. We observe that the demagnetization field contribution to the gap size is negligible, consequently revealing the linear dependence on the strength of the magnetic anisotropies. Our experimental results are accompanied by calculations based on a Landau-Ginzburg [14,24] energy functional, which includes the cubic magnetic anisotropy term that effectively reproduces the experimentally observed features. Discrepancies between the experimentally and theoretically observed size of the hybridization gap can be tentatively attributed to the influence of a gradient term of the exchange interaction.

EXPERIMENTAL PROCEDURES
Material and phase diagram. Copper-oxoselenite, Cu 2 OSeO 3 is a chiral magnetic insulator with small Gilbert damping which crystallizes in a non-centrosymmetric cubic structure hosting 16 copper Cu 2+ ions per unit cell. The magnetization properties of the crystal are owed by the formation of tetrahedral clusters by four Cu 2+ spins in a 3-up-1-down spin configuration, behaving as a "spin-triplet" with the total spin S = 1. For the detailed visualization of the magnetic structure, see Ref. [5]. The spin-spin interactions in the system are governed by the symmetric Heisenberg interaction, which favours collinear spin structures, and an anti-symmetric Dzyaloshinskii-Moriya To put our experimental results into context, we first show the generic magnetic phase diagram of Cu 2 OSeO 3 for static magnetic fields applied parallel to the 100 crystallographic direction, as illustrated in Fig. 1(a). Above the transition temperature T c ≈ 58 K, the system is paramagnetic. Below T c and without applied field, the system is described by a multidomain helical phase. When the magnetic field exceeds a critical field value H c1 , a phase transition into the conical phase occurs. By further increasing the magnetic field strength, the cone angle between the spins and the field axis decreases, resulting in a fully polarized state at H c2 . In addition to these states, a topological nontrivial HTS phase exists in a small temperature and field pocket, close to T c . This phase pocket can be extended considerably using various procedures, such as applying particular temperature-quenching techniques, applying pressure or using thin-film crystals where the thickness is of the order of several pitch lengths of the helical spiral (≈ 61 nm in Cu 2 OSeO 3 ) [12,13,[25][26][27][28][29][30]. Recently, it has been demonstrated that the second, low-temperature, skyrmion phase can be nucleated after traversing the so-called tilted conical phase and following a unique magnetic field-cycling protocol [21,22]. To enter the tilted conical phase, a state where the helix pitch tilts away from the field axis, the temperature needs to be decreased below T ≈ 30 K.
Then, as the applied field decreases further, the LTS phase emerges. Similarly, increasing the magnetic field from zero magnetic field, the tilted conical phase may form upon increasing magnetic field, again leading to the LTS generation. In both cases, the volume fraction of the LTS can be drastically increased by cycling the magnetic field within magnetic field boundaries where the LTS may be expected, as we will demonstrate below.
Experimental setup. We used two differently shaped bulk Cu 2 OSeO 3 crystals for our experiments to vary the demagnetization fields in the samples. One sample is a cuboid with dimensions is introduced into a closed-cycle cryostat, inserted into the gap of an electromagnet, as illustrated in Fig. 1(b). The magnetic field is aligned normal to the plane of the CPW, i.e. along the 100 direction of the Cu 2 OSeO 3 crystals. Since the samples are larger than the centre line (width 1.25 mm) of the CPW, magnetization dynamics are excited with both in-plane and out-of-plane oscillating magnetic fields, resulting in excitation of gyration and breathing modes respectively.
We have employed a microwave absorption spectroscopy to probe and identify these states, as demonstrated by Y. Onose et al. [9] The microwave magnetic fields are generated in the range of 1 -6 GHz in steps of ≈ 3 MHz by connecting the output port of a vector network analyzer (VNA) to the CPW. We measure |S 11 | or |S 21 |, the reflection and the transmission signal, respectively, of the microwave absorption spectra as a function of the applied field and subtract a background signal by re-measuring the spectra in a high magnetic field well above H c2 , i.e. in the field polarized state (see Supplementary Fig. S1 for details). The measurements were carried out in a temperature range of T = 5 − 60 K.
Field cycling protocol. After high-field cooling to the desired temperature from above the ordering temperature, the applied field is first reduced from a large magnetic field H init , within the field-polarized state to the upper boundary of the desired cycling window, H high where we subsequently start the cycling procedure. In this process, the magnetic field is repeatedly decreased and then increased between H high and H low in steps of ∆H. In this context, one complete cycle means that the magnetic field is ramped down to H low and subsequently increased again to its original value, i.e., H high → H low → H high . Upon completing n number of cycles, the magnetic field is either decreased to zero-field, H dec or increased to H init . at particular fields, separated by steps of 1 mT. For example, scans 0 and 1 correspond to single measurements at H init and H init + 1 mT. In the cycling region, the data was recorded every 10 th cycle, thus separated by 10 mT. Therefore, each lower vertex of the triangular-shaped signal corresponds to a completion of 10 cycles.
As the field decreases from H init to H high , an almost linear slope assigned to the Kittel resonance mode in the field-polarized regime is observed. Note that less prominent modes above and below the Kittel mode are standing spin waves across the whole Cu 2 OSeO 3 crystal, which can be observed due to the low damping of the material at low temperatures [31]. The gradient change around scan = 40 indicates the phase transition to the conical phase, exhibiting a frequency increase with decreasing magnetic field. At around scan = 60, the spectrum reveals a significant discontinuity, which we attribute to the onset of the tilted conical phase [21]. In the cycling regime, signatures of emergent resonant modes are observed at lower frequencies with their intensities increasing with the number of cycles. Completing the cycling protocol, both gyration and breathing LTS modes are enhanced and clearly observed at low frequencies when further decreasing the field towards zero magnetic field.

RESULTS AND DISCUSSION
Let us now discuss the magnetic excitations detected on the plate-shaped sample. Figs. 2(a -c) show the combined excitation spectra of H inc and H dec , measured for a various number of cycles n as a function of the magnetic field at T = 6 K. In the upper panel of Fig. 2, we show data with n = 0, i.e. without field-cycling. In reminiscence of the conical state, two modes may be distinguished, exhibiting an increase in absorption frequency for decreasing magnetic field. These resonance branches, which are expected to couple to the perpendicular excitation field with respect to the helix pitch vector, are therefore identified as ±Q modes [3,21]. Their extracted resonance positions are marked by yellow dots. Next, we investigate the evolution of the excitation spectrum with the number of field-cycles n = 100 ( Fig. 2 We notice that after their initial regime up to 500 cycles, both parameters steeply increase with cycle number up to around 4000 number of cycles and then start to saturate afterwards. Since the remnants of the tilted conical phase strongly contribute to the absorption spectra for small numbers of cycles (<500), the fitted amplitudes and resonance positions might be inaccurate, as indicated by the shaded region. In order to support these results theoretically, we employ the previously established phenomenological model for chiral magnets, which can be found in various publications [11,14,21].
The theoretical approach is based on the Ginzburg-Landau energy density, which consists of the following contributions: The first energy term interaction strength was already determined in prior reports and is set to τ ≈ 0.88 accordingly [11,14].
Since the focus of this work lies in the investigation of the LTS phase, we likewise limit our analysis to the eigenmodes and eigenresonances of the topologically nontrivial skyrmion states. In the following, the calculations are performed either for a spherical sample with demagnetization i.e. parallel to the 100 direction. The field boundaries were chosen to cover the range in which the skyrmion phase is energetically more favorable than the topological trivial states. Nevertheless, it was already shown that the LTS phase may remain even in smaller fields, merging into the elongated skyrmion phase, see Ref. [21]. The LTS phase is characterized by three distinct modes of the skyrmion lattice, the counterclockwise (CCW) gyrational mode, which decreases in frequency with decreasing magnetic field (violet), the breathing (BR) mode, which increases in frequency with decreasing magnetic field (red) and the clockwise (CW) gyrational mode again with decreasing frequency as a function of decreasing magnetic field (green). We find excellent agreement between the calculated excitation spectrum and the experimentally detected dynamic modes after field cycling.
A particular feature that becomes very clear from the theoretical analysis is the opening of a gap for the BR mode (red) which hints to a hybridization and which is also observed in the experiments. This hybridization feature is not observed in the HTS phase and will be discussed in the following. By introducing the cubic anisotropy term in our micromagnetic calculations, we observe a hybridization of the breathing, clockwise and counterclockwise mode with additional dark skyrmionic modes, i.e. modes that cannot be observed in the experiment due to low spectral weight. These anti-crossings differ significantly in magnetic field position, spectral weight and gap size and the dominant hybridization observed also in our experiments is the one of the BR mode. By means of real space images for different times t (see Supplementary Fig. S3 for details), those dark modes are identified as higher-order clockwise modes. Based on the number of nodes, we refer to these modes as Sextupole (6 nodes    Next, we are able to account for the size of the hybridization gap g quantitatively by extracting the minimal frequency difference between the breathing and octupole mode branches from the theoretical and experimental results. The calculated resonance frequencies, normalized by are given in dimensionless units, with internal critical field H int c2,0 determined without cubic anisotropy. For small K, we assume H int c2,0 ≈ H int c2 , in order to obtain an estimate for the physical units. The internal critical field is connected to the saturation magnetization M s by the internal susceptibility χ int con in the conical phase [11,32,33], In Fig. 5, the calculated field dependence of the excitation frequencies (black dots) is shown for K = 0, K = 0.0002 and K = 0.0004. While the colors of the dots indicate the excitation geometry with respect to the applied magnetic field direction (blue: in-plane, red: out-of-plane), the size of the dots reflects the spectral weight of the resonance modes. As already anticipated before, it can be seen, that including the cubic anisotropy term into the model results in hybridizations of different modes, which interaction strengths depend on the anisotropy value.
In order to extract the temperature dependence from our calculations quantitatively, it is sufficient to leave only the anisotropy value as a free parameter. The Ginzburg-Landau coefficient r 0 , the measure for the temperature, enters our model as a scaling factor proportional to the saturation magnetization M s . Therefore, by extracting the temperature-dependent parameters H int c2,0 and M s , the calculated spectra can be given in physical units at the corresponding temperature.
The anisotropy constantK is determined with the help of the results obtained in [14,15]. It was found that, exceedingK th c ≈ 0.0001, the LTS phase would stabilize in the theoretical model. Note that in order to avoid confusion, the dimensionless anisotropy constant is marked by a tilde, here. This value corresponds to the dimensionless value obtained from the experiments, with the anisotropy constant K σ,c given in units of energy density [14]. Considering now the observed anisotropy values in [15], the theoretical model can provide an estimate of the gap size, as a function of temperature.
The experimental data is evaluated the same way as described above. The resonance positions of the breathing and octupole mode in the vicinity of the hybridization were extracted from two Lorentzian peak fits. Calculating the frequency difference ∆f between both branches as a function of field, we define the minimum value as the gap size g. In Fig. 6(a) the fitting results are shown for the measurements performed at 6 K. At this temperature, g obtained from ∆f , given in the inset, reads 450 MHz.
For a conclusive temperature dependence of the gap size, we extend our analysis to the temperatures covered by our experiments. Our findings from both measurements and calculations are collated in Fig. 6(b). In both crystals studied, we observe a similar trend qualitatively, namely a decrease of g with increasing temperature. In order to serve as a guide to the eye, a grey line is added to the figure. Its intersection point is chosen to be consistent with the temperature-dependent anisotropy measurements from Halder et al. [15]. The theoretically estimated gap size, converted to physical units, is given by the blue line, reminiscent of the behaviour in the experimental data.
Its slope, however, differs by approximately a factor of 2.3 with respect to the experiment. A small deviation between experimental and theoretical data is induced by the conversion from dimensionless to physical units. However, this error is assumed to be small compared to the factor mentioned above, leaving the origin of this discrepancy open. In order to further investigate these findings, we included an additional exchange anisotropy term into the free energy functional. This anisotropy is allowed by the symmetry of the P 2 1 3 space group, and was already taken into consideration in [17,34]. Depending on the sign of C, the helix pitch aligns either along 111 (C < 0) or 100 (C > 0). Fixing all parameters apart from C, we calculate the resonance spectra and therefore also the gap size, in the same manner as before, as a function of the anisotropy strength. The extracted values are summarized in (Fig. 7). The blue line corresponds to the results already shown in Fig. 6 induced by the cubic anisotropy in this case, were already reported in [14,21]. Nevertheless, despite this additional energy contribution, the resonance spectra are left essentially unchanged, as shown in Fig. 7(b) -(c). The calculations were performed for Fig. 7(b) K = 0.0002, C = 0 and Fig. 7(c) K = 0.0002, C = 0.2, respectively, illustrating the significant impact on the hybridization process. These results suggest that further anisotropy terms [20,21], like the one discussed above, might play a non-negligible role and need to be considered for a comprehensive analysis of the temperature dependence of the hybridization.

SUMMARY
In conclusion, we report a comprehensive study of the nucleation of the low-temperature skyrmion lattice that is stabilized by the magnetocrystalline anisotropy of the sample via a fieldcycling protocol and systematically trace the characteristic hybridization in the breathing mode.
By employing two different sample geometries, we were able to remove contributions from the demagnetizing field and investigate the dynamics of low-temperature skyrmion modes for various temperatures. Our analysis of the hybridization gap highlights that the cubic magnetocrystalline anisotropy terms govern the origin of the hybridization following a particular symmetry selec- Note that a clear hybridization gap, g, is observed for both samples, decreasing with increasing temperature.  Fig. S1 (a) shows typical reflection spectra S 11 obtained directly from the VNA for the plateshaped sample at T = 6 K. For systematic measurement analysis and accurate determination of excitation dynamics, the data was processed over three distinctive stages. Firstly, we carefully remove the excessive noise in the signal utilizing an algorithm involving a numerical filter. In order to preserve all necessary information from the data obtained, a lossless analytical Savitzky-Golay filter [35] was employed, which the filter recalculates the signal from the convolution coefficients by fitting successive data points by a polynomial function using linear least squares. The filtered signal, S 11 is shown in Fig. S1(b). Secondly, the spectra in the region of interest were subtracted by a background signal (µ 0 H = 300 mT) obtained at a high magnetic field, as shown in Fig. S1(c).
Lastly, Fig. S1(d) shows the resulting signal after a full conversion into a linear scale, i.e., |S 11 | = 10 ∆S 11 /10 . The data were then fitted using multi-peak Lorentzian to extract the parameters of interest. We have used the same systematic protocol for all data presented in this paper.

S.2 ADDITIONAL CYCLING DATA
Here, we present additional spectra data obtained respectively, at 800 and 20000 field cycles as shown in Figs. S2(a -b) obtained at T = 6 K. A vertical line is drawn to indicate the approximate field (µ 0 H ≈ 75 mT) where hybridization is clearly seen for both number of cycles. In order to investigate the evolution of the dynamic modes, we further plot the line-cut data at this field for different numbers of field cycles in Fig. S2(c -g).
Our results highlight that the intensity of the hybridization becomes more visible with the number of accumulated field cycles which must be connected to the volume fraction of the LTS

S.3 THEORETICAL MODEL
The theoretical analysis presented in this report draws upon the well established Ginzburg-Landau theory for chiral magnets [11,14,21]. The energy functional F consists of the following energy terms, The parameters are the Ginzburg-Landau coefficients r 0 and U , the strength of the exchange J, Dzyaloshinsky-Moriya D, dipolar τ interactions, with demagnetization tensor N (tr(N ) = 1), magnetic field B and the cubic anisotropy constant K. It was already shown that the additional anisotropy term is sufficient to stabilize the LTS phase [14] and to quantitatively reproduce the changes in the resonance spectra, with respect to the high-temperature skyrmion phase [21].
The magnetization dynamics are described by the lossless Landau-Lifshitz equation of motion [36]: Here, γ = gµ B / denotes the gyromagnetic ratio and B eff = −δF/δM the local effective magnetic field, given by the derivative of the total energy functional with respect to the magnetization.
For further analysis, the free energy F and the magnetization M are divided into a static and dynamic component, Inserting this ansatz into the equation of motion 7 and only keeping terms linear in F (t) and δM(t), the resonance condition reads [32], with W = γM mf × χ −1 0 and χ −1 By performing the calculation of Eq 9 in momentum space, we are able to determine the normalized eigenvectors v α j (k) and eigenfrequencies ω α (k) of the given matrix W numerically, with internal indices j ∈ [1, 2, 3] and momentum index k. A more detailed derivation can be found in [11,14].
Based on the rather large dimensions of the CPW, we assume the driving field to be homogeneous in the performed calculations, therefore setting k = 0. The demagnetization factors of the samples under investigation read approximately N x = 0.17, N y = 0.12 and N z = 0.71 for the platelet and N i = 1/3 for a sphere, which we take as an approximation of the cuboid, where the external magnetic field is applied along the z-direction.
Next to the breathing, clockwise and counterclockwise mode, we focus our study on three additional modes, ranging in the same frequency regime. For the identification of the individual modes, real space images of the normalized dynamics magnetization δM z are calculated for different times t. The results are shown in Fig. S3. Based on the time evolution and the number of nodes, these modes are identified as higher-order clockwise excitations and are referred to as Sextupole (6 nodes), Octupole (8 nodes) and Dectupole(10) mode.
In order to verify the repulsive character of the hybridizing modes, we calculate the inner product of the interacting resonance branches for K = 0 and K = 0. For a homogeneous driving field, the internal product is given by [11], with reciprocal lattice L R . The results are presented in Table II of the main text, revealing the effect on the orthogonality of the eigenvectors of the hybridizing modes.
From our calculations, we observe that neither the shape of the sample nor a finite k vector changes the gap size, significantly (not shown here). In order to further investigate the discrepancy between the experimentally and theoretically obtained gap size we added an additional anisotropy into the free energy functional. Depending on the sign of the prefactor C, the helix pitch aligns either along the 111 (C < 0) or the 100 (C > 0) direction. In order to determine the effect of the exchange anisotropy on the hybridization gap size, we calculate the resonance spectra for various C values. For more details we refer the reader to the corresponding section in the main text.