Cu₂OSeO₃ represents a unique example in the family of B20 cubic helimagnets [1, 2] exhibiting two well-defined skyrmion pockets at the temperature-magnetic field phase diagram [3, 4] (Figure 1A). The “high-temperature” (HT) pocket is commonly referred to as the A-phase: it is located at the boundary with the paramagnetic state and is also intrinsic to other cubic helimagnets such as MnSi [8, 9] and/or FeGe [10, 11]. The “low-temperature” (LT) pocket, however, arises around zero temperature and only in this particular chiral magnet. It forms exclusively for magnetic fields applied along the easy crystallographic 100 directions [6, 12].

Within a small pocket of A-phase, skyrmion lattice (SkL) appears spontaneously and its stability is commonly attributed to the thermal fluctuations, which “work” at relatively high temperatures [9, 13]. Moreover, the temperature versus field regime of the HT-SkL is almost isotropic for different field directions with respect to the crystal lattice. The small angle neutron pattern has a sixfold intensity for a random orientation of a sample in the A-pocket [9]. At the same time, the boundaries of the A-phase can be drastically changed by applying pressure [14], electric fields [15–18], chemical doping [19] or uniaxial strains [20, 21].

The principle of HT-SkL stabilization rests on a specific field-driven evolution of the energetic difference to SkL’s main competitor, the conical spiral [22]. The difference between the energies of the hexagonal skyrmion lattice W _SkL and the conical spiral W _cone , ΔW _min = W _SkL − W _cone , has minima along a curve ξ(T) exactly for those magnetic fields that stabilize the A-phase [22–24] (Figure 1A). In the following, the conical spiral is defined as a solution of the isotropic Dzyaloshinskii model (1) with the propagation direction along the magnetic field in which the magnetization rotation retains single-harmonic character, but becomes additionally distorted by different anisotropic contributions (4) (Figure 1C).

Remarkably, the same stabilization principle holds for the LT-SkL [3]. The stabilization mechanism, however, is represented by a well-known cubic anisotropy with easy 100 axes. Neutron diffraction measurements in Ref. [3] confirmed this mechanism and revealed the formation of skyrmion states in large areas of the magnetic phase diagram, from the lowest temperatures up to the A-phase along the line ξ(T) (Figure 1A). LT-skyrmionic scattering, however, shows a clear hysteretic behavior: it concentrates around the H _c2 line when the magnetic field increases from zero to above H _c2, and around the H _c1 line when the magnetic field subsequently decreases from above H _c2 to zero (as indicated by small green arrows at the sketch of the experimental phase diagram, Figure 1A). Moreover, the LT-skyrmion correlations underlie a ring-shaped scattering, which gains the strongest intensity in the LT-SkL region and the reduced radius as compared with the spiral peaks.

The role of cubic anisotropy in the LT-region is not limited only to the SkL stabilization, but is also manifested by many other remarkable phenomena:

(I) In particular, cubic anisotropy fixes the propagation directions of spiral states along easy anisotropy directions below the critical field H _c1 and underlies multidomain spiral states: for MnSi [8, 9] spirals propagate along the easy 111 directions, for Cu₂OSeO₃ [1]—along 100 . If the propagation vector of a spiral state is perpendicular to an applied magnetic field, we will call such states transversal spirals (Figure 1E). If the wave vector composes some angle with respect to the field direction below the H _c1 line, such spirals will be called oblique spirals.

The multidomain spiral state of Cu₂OSeO₃ below H _c1 for H‖[001] consists of spiral domains with mutually perpendicular wave vectors: one domain represents the most energetically favorable conical state (Figure 1C), and the others—metastable transversal spirals (Figure 1E), which eventually flip along the field direction. According to the theoretical predictions, the domains of transversal spirals do not reappear with the lowering field, which is also indicated by the experimental results in Ref. [12]. Thus, the field H _c1 is not shown at the theoretical phase diagram (Figure 1B) for this field direction [22, 25].

Multidomain spiral state of MnSi for the magnetic field applied along the easy anisotropy direction, H‖[111], consists of energetically degenerate oblique spirals leaning towards the field and thus deflecting from their zero-field propagation directions along the 111 axes. Since such spirals represent global minima of the system, the multidomain state will be repopulated even when lowering the magnetic field.

To discern the spiral behavior at the critical fields, we will use the notion H c 1 * in the former case and H _c1 in the latter case.

(II) Interestingly, the spiral flips were theoretically predicted to occur also for higher fields (near H _c2) and relatively strong values of cubic anisotropy [22], but they have not been yet identified experimentally. Although the conical phase co-aligned with the field is considered to be the most favorable state above H _c1 that benefits from the Zeeman interaction, it may jump away from the field to be able to embrace a larger number of easy anisotropy axes.

For example, for H ‖ 111 or H ‖ 110 and easy cubic axes 111 , the wave vector jumps closer to the hard axes 100 (see for details Refs. [22, 25]). Theoretically, if the spiral jumps from a conical into a tilted spiral (we will use this definition to distinguish such spirals from the mentioned oblique spiral states at low fields), the canting angle with respect to the field remains almost fixed by the strong cubic anisotropy [22, 25].

No tilted spiral states near H _c2, however, were identified theoretically for the easy cubic axes 100 and the field along the high-symmetry directions 111 , 110 , and 100 .

(III) For moderate and small anisotropy values, cubic anisotropy defines the character of the phase transition between the conical and the homogeneous states at the critical field H _c2: first-order phase transition (FOPT) occurs when the magnetic field is applied along the easy anisotropy axes (in the case of Cu₂OSeO₃, for H ‖ 100 ) whereas it retains its second-order character for the field along the hard anisotropy axes 111 . From the ascent of the critical line H _c2 from zero temperature and onwards to T _C , as indicated at the sketch of the experimental phase diagram (Figure 1A), we can deduce that in Cu₂OSeO₃ the cubic anisotropy exhibits a pronounced temperature dependence and gradually weakens as shown by the dashed blue line at the theoretical phase diagram (Figure 1B).

(IV) Additional deviation from the generic phase diagram of cubic chiral magnets [9] is manifested in Cu₂OSeO₃ by the tilted spiral state (TS), which appears for the same direction of the magnetic field H ‖ 100 and coexists with the stability region of the LT-SkL [5]. As stated before in II, such a TS cannot be stabilized solely by the cubic anisotropy. In Ref. [5] it was found that the TS is presumably stabilized by the competing effect between cubic and exchange anisotropies. The exchange anisotropy in this case has easy 111 axes and is sufficiently strong to deflect the wave vector of the conical state towards these crystallographic directions. With increasing temperature, apparently, the exchange anisotropy also fades away to keep the ratio of anisotropy constants below the threshold limit when spirals at zero field would flip along the easy 111 directions of the exchange anisotropy. Besides, the tilted spiral state is hailed as a nucleation source of LT-skyrmions [6], which can form within domain boundaries between different TS-domains and thus leads to the aforementioned hysteretic behavior of the LT-SkL.

Such a pronounced temperature dependence of anisotropy coefficients alongside with their complex interplay constitutes the main obstacle on the way to construct the complete quantitative model of a bulk cubic helimagnet Cu₂OSeO₃. In previous reports, only separate assessments of anisotropy values have been undertaken without any regard to their temperature dependence and interplay [26].

In the present manuscript, we underline distinctive features of tilted spiral states that stem from the effect of the interplaying cubic and exchange anisotropies. We discuss these results in the frame of the phenomenological theory introduced by Dzyaloshinskii in an attempt to establish a quantitative comparison between the model and the experimental results for a bulk cubic helimagnet Cu₂OSeO₃. In particular, we find a crossover behavior of a tilted spiral state driven by the varying ratio of anisotropy constants: i) in the case of strong cubic and exchange anisotropies, the TS directly jumps into the homogeneous state and exhibits almost no rotation; ii) for weaker anisotropies, however, the TS can smoothly rotate towards the easy axes of exchange anisotropy and even shows a reverse rotation back to the easy axes of cubic anisotropy. We argue that experimentally these two regimes are defined with respect to the critical temperature T = 18K (point A at the phase diagram in Figure 1A), which introduces an important threshold. By a systematic study of the stability and ordering of the low temperature magnetic states in Cu₂OSeO₃ done in Refs. [3, 5, 7] we conclude that the former case of the TS behavior is inherent for low temperatures (T = 2K) whereas the latter case allows to explain the experimental observations in the vicinity of the point A.

Additionally, we show that skyrmion states may also become oblique with respect to an applied magnetic field, which, however, occurs for the field directions 111 and 110 making skyrmions metastable states. Their field-driven jump-like reorientation is accompanied by a gradual rotation perpendicular to the field and may come into a conflict with versatile oblique spiral states, which makes the experimental phase diagram rather complex.

Furthermore, by comparing the experimental and theoretical phase diagrams, we are able to delineate the ranges of anisotropy constants, which may subsequently serve as a guide for future experiments. Although such theoretical undertakings require some additional fine tuning, our approach and the reasoning behind it provide a viable strategy for an in-depth and quantitative understanding of chiral magnets in view of tailoring their properties for future applications.

Within the phenomenological theory introduced by Dzyaloshinskii [27] the magnetic energy density of a bulk non-centrosymmetric ferromagnet with spatially dependent magnetization M can be written as W 0 M = A ∑ i , j ∂ m j ∂ x i 2 + D w D M − M ⋅ H (1) where A and D are coefficients of exchange and Dzyaloshinskii-Moriya interactions (DMI); H is an applied magnetic field; x _i are the Cartesian components of the spatial variable. w _D is composed of Lifshitz invariants L i , j k = M i ∂ M j / ∂ x k − M j ∂ M i / ∂ x k (2) that are energy terms involving first derivatives of the magnetization with respect to the spatial coordinates. In the following, all calculations will be done for cubic helimagnets with w D = m ⋅ ∇ × m (3) although the results may be applied for magnets with other symmetry classes [28] including different combinations of Lifshitz invariants. W ( M ) includes only basic interactions essential to stabilize skyrmion and spiral states and specifies their most general features attributed to all chiral ferromagnets.

In the forthcoming calculations, we use reduced values of the spatial variable, x = r/L _D , where L _D = A/D is the periodicity of the modulated states. m is the unity vector along the magnetization vector M = m M , and h = H / H D with μ ₀ H _D = D ²/(AM) being the reduced value of the applied magnetic field.

For the cubic helimagnet Cu₂OSeO₃, we supplement the isotropic energy density (1) by the exchange and cubic anisotropic contributions [26, 29], Φ a = b e a T ∑ i ∂ m i / ∂ x i 2 + k c T ∑ i m i 4 , (4) where b _ea = B _ea /A and k _c = K _c A/D ² are reduced anisotropy constants, which are in general temperature-dependent. The constants k _c and b _ea are typically one order of magnitude smaller than the exchange stiffness A and are conventionally considered as a third level of hierarchy of energy scales following the exchange and DM interactions.

The theoretical explanation of SkL stability (Figure 1B) by the anisotropies (4) is mainly based on the effect imposed on one-dimensional spiral states. In fact, the ideal magnetization rotation in the conical state can be impaired by the easy and hard anisotropy axes for specific directions of the magnetic field. Through this mechanism, skyrmions, which are more resilient to anisotropy-induced deformations, due to their two-dimensional nature, gain stability [22].

In the following simulations, we restrict ourselves to k _c > 0 with easy 100 axes and b _ea < 0 with easy 111 axes, which is the presumable layout for the aforementioned Cu₂OSeO₃ [3]. The emphasis will be given to the results with H‖[001]. We also neglect the influence of dipole-dipole interactions due to the volume charges presumably formed within different oblique states. We assume that the DM interactions suppress demagnetization effects and are the main driving force leading to the magnetization rotation and to the equilibrium periodicity. The influence of dipole-dipole interactions on the effects found in the present manuscript will be considered elsewhere.

The Euler-Lagrange equations derived from the energy functional (1) are non-linear partial differential equations. These equations have been solved by numerical energy minimization procedure using finite-difference discretization on grids with adjustable grid spacings and periodic boundary conditions. Components of the magnetization vector m have been evaluated in the knots of the grid, and for the calculation of the energy density (1) we used finite-difference approximation of derivatives with different precision up to eight points as neighbors. To check the stability of the numerical routines we additionally refined and coarsened the grids. For axial fields, we used grid spacings Δ _y ≈ Δ _x so that grids are approximately square in the xy plane in order to reduce any artificial anisotropy incurred by the discretization. The final equilibrium structure for the modulated states was obtained according to the iterative procedure of the energy minimization using simulated annealing and a single-step Monte-Carlo dynamics with the Metropolis algorithm. In detail, numerical methods on energy minimization procedure are described in, e.g., Ref. [30] and hence will be omitted here.

To avoid an impediment introduced by the periodic boundary conditions, which would arise due to tilted/oblique spiral states in case of using a three-dimensional numerical grid, we perform two-dimensional simulations. For these, we write the energy density in a coordinate system ( x ̃ , y ̃ , z ̃ ) connected with the wave vector of a TS and a corresponding plane of rotation (Figure 2A). First, we rotate the coordinate system (CS) around z by an angle ψ. Then, we perform a subsequent CS rotation around y ̃ in such a way that the wave vector of a spiral state, which points along x ̃ , makes an angle α with the z-axis. For a TS, we assign ψ = π/4 and thus can address the TS structure in the plane ( 1 1 ̄ 0 ) . The energy density is then minimized with respect to α and a tilted spiral with an optimal α is found. In a general case, we will also minimize the energy density of a spiral state with respect to ψ.

For the present setup of easy 100 anisotropy axes and h‖[001], the tilted spiral state was found by both theoretical and experimental means. Its behavior and occurrence at the phase diagram are pronouncedly temperature-dependent and exhibit well-distinguished “regimes”.

First of all, to stabilize TS in the vicinity of the critical field H _c2, the cubic anisotropy k _c must overcome some threshold value (the exact value will be computed elsewhere). A relatively weak cubic anisotropy k _c = 0.05 does not lead to any intermediate tilted spiral, although, according to the theoretical phase diagram in Refs. [3, 25], it is still able to stabilize a LT-SkL (Figure 1B). An increasing exchange anisotropy in this case would only switch the conical phase with α = 0 immediately into a 111 oblique spiral even for zero field, which then aligns with the field at the H _c1 value, a situation, which was never observed experimentally in Cu₂OSeO₃. Higher values of cubic anisotropy, on the contrary, are found to be responsible not only for the stability of the tilted spiral state, but also for the crossover of its behavior.

In the following, as instructive and representative examples, we consider two cases with k _c = 0.1 (as becomes evident later, this value corresponds to higher temperatures in Cu₂OSeO₃) and k _c = 0.15 (lower temperature range, see Sect. V. for details). Further increase of k _c leads to a high-field jump from the conical into a transversal spiral, as featured by the theoretical phase diagram (see point C in Figure 1B and the pink-shaded region of the transversal spiral). Since such a limiting case was also not observed experimentally in bulk Cu₂OSeO₃, we may restrict ourselves to a finite range of anisotropy coefficients, which can be subsequently narrowed down for a better fitting with the experimental results.

In the same way as it was implemented for tilted spirals, the spatial orientation of SkLs is specified by the competing effect of easy anisotropy axes and an applied magnetic field. For example in zero magnetic field and k _c > 0 (in the following, k _c = 0.1), SkLs are metastable states and occupy the crystallographic {001} planes with easy axes ⟨001⟩ of the cubic anisotropy (Figures 6, 7). For h‖[001], the equilibrium position of SkL was found to be codirectional with the applied magnetic field. Consequently SkL gains stability in the vast area of the theoretical phase diagram (Figure 1B). Two energetically close SkL minima were found to occur in the transversal (001) plane. The first minimum corresponds to the easy ⟨100⟩ cubic axes pointing along the diagonals of the hexagonal SkL and along their apothems. In the second minimum, the hexagonal lattice is rotated by the angle π/4. The cores of skyrmions in both states become square shaped with the tendency either to elongate or to shorten along particular directions (see for details Ref. [25]). In the experiment [3], one can discern 12 peaks originating from these two SkL domains within the ring of scattering. For larger anisotropy values, the easy cubic axes along the diagonals may induce an elliptical instability of the SkL similar to that of isolated skyrmions [31, 32] and trigger the phase transition into the helical state. Thus only the second SkL minimum is preserved. At the phase diagram of Figure 1B, the line a − b separates the two skyrmion lattice phases. For larger values of b _ea , SkLs may occupy {110} planes, which include the easy axes of the exchange anisotropy ⟨111⟩. However, even in a small magnetic field, SkL planes align perpendicular to the field.

For other directions of the field and k _c > 0, the SkL still remains a metastable solution. In an applied magnetic field h‖[110] (Figure 6A), the SkL rotates and orients perpendicular to the field, i.e., occupies the (110) plane. Whereas for b _ea = 0 such a rotation is smooth as shown by the energy density curves (Figure 6C), it becomes abrupt under the influence of exchange anisotropy. It is obvious that such SkLs’ jumps into the (110) plane are dictated by the easy axes ⟨111⟩ of the exchange anisotropy alongside with the Zeeman interaction. Figure 6D shows the angle ϕ between the field and the normal n to the SkL plane for different values of b _ea . The equilibrium angle ϕ is obtained by the energy minimization procedure, in which we smoothly rotate n in the plane (001). For b _ea = −0.15, the SkL is located in the plane (110) even for zero magnetic field, i.e., ϕ = 0 in the whole field range. Dashed lines indicate that the tilted SkLs exist as local energy minima even above the transition fields into the state perpendicular to the field.

For the magnetic field applied along [111] direction, the reorientation processes of LT-SkLs are more involved. Figures 7A,B show the details of SkL jumps for the “pure” effect of cubic anisotropy. SkLs, located in {001} planes, smoothly rotate towards {111} planes with increasing magnetic field. Interestingly, a SkL in the plane {110} represents a local energy minimum as shown by the first panel of Figure 7B. Two energy minima merge when the SkL aligns perpendicular to the field (third panel). For a dominating role of the exchange anisotropy (Figures 7C,D), SkLs occupy {110} planes and then gradually rotate towards {111} planes perpendicular to the field. For the chosen set of the anisotropy constants, the planes {100} do not even represent a local energy minimum (first panel in Figure 7D). For intermediate values of exchange anisotropy, however, SkLs undergo the following reorientation processes: first, they jump into a position close to {110} planes and then rotate towards {111} planes (Figure 7F). Such jumps become apparent at the field dependencies of ϕ(h) for fixed values of b _ea (Figure 7E).

To establish a quantitative phenomenological model of the bulk helimagnet Cu₂OSeO₃, one should determine all material constants entering Eqs. 1, 4. In this sense, hierarchically the most important Heisenberg exchange and Dzyaloshinskii-Moriya interactions may be evaluated at the atomic level, using, e.g., ab initio density functional theory (DFT) calculations. A multi-scale approach pursuing this goal is described in detail in Ref. [26] and Table 4 in the aforementioned paper gives the material constants A, D and K _c . Alternatively, the strength of the DMI together with the exchange interaction is known to govern the pitch of the spin spiral, the diameter of a single skyrmion, as well as the saturation field (H _c2) required to fully align the magnetic moments along the external field, and thus can be roughly evaluated. The obstacle to finalize the quantitative model, however, lies in the fact that anisotropic contributions with non-dimensional anisotropy coefficients k _c and b _ea usually appear together in the magnetization processes and moreover exhibit pronounced temperature dependence. In the following, we discuss qualitatively a strategy that could be utilized to complete a quantitative model for Cu₂OSeO₃. We will use some behavioral patterns drawn from previous theoretical results.

First of all, we notice that the point A at the experimental phase diagram—an intersection point between the curve H _c2 and an upper boundary of a LT-SkL pocket–demarcates two temperature regimes: for lower temperatures, SkL exists above H _c2 and its nucleation is believed to be directly related to the stability region of a tilted spiral [3, 4]; for higher temperatures, SkL dissolves before the cone saturation field H _c2; at the same time, identification of the tilted spiral in this temperature interval bumps into the accuracy of experimental procedures. According to the theoretical phase diagram constructed in Ref. [25], the point A corresponds to the value k c ( A ) ≈ 0.06 (Figure 1B). The minimal anisotropy value k _c needed to stabilize skyrmions is computed theoretically, k c ( B ) ≈ 0.04 , which would correspond to the farthest temperature point of the LT-SkL pocket in Figure 1A. Such a cubic anisotropy, however, will not lead to a tilted spiral state, which sets in only for k _c > 0.05. At the same time, we argue that the value of k _c must be smaller than 0.15, since in the experiments no tilted spiral state was observed for H‖[110] (although it is predicted in Figure 5). The value of the exchange anisotropy b _ea can be then finely tuned by comparing different characteristics of the tilted spiral from the experiment and the theoretical modeling. An exhaustive analysis will be done elsewhere. Here, we just claim that the two regimes of the tilted spirals considered before roughly correspond to the experimental results for 2K and 18K, respectively.

A systematic study of the stability and ordering of the low temperature magnetic states was done in Refs. [3, 5, 7]. Analysis of experimental data for low temperatures (e.g., T = 2K) prompts the following conclusions: i) an experimental value of α _max = 0.61 cannot be explained by the theoretical data for k _c = 0.1, since theoretically only a smaller angle value can be achieved; on the other hand, for k _c = 0.15 such an angle value is readily reached on one of the plateau-like dependencies α(h) (Figure 3D) for b _ea < − 0.15; ii) comparing the experimental field interval of tilted spiral existence ΔH ≈ 0.59 with Figure 3F one may estimate b _ea ≈ − 0.18. In any case, theoretical results for k _c = 0.15 allow to qualitatively address the following experimentally observed phenomena in Ref; [7]: i) first-order phase transition with respect to the homogeneous state, which also becomes apparent in the coexistence of conical and tilted spiral states; ii) rather abrupt cone evolution into a tilted spiral state, which is identified as a blind zone, for which a tilted spiral angle does not exist (as seen in Figure 2E); iii) comparable angles and field intervals of spiral canting followed from the theory and the experiment.

Analysis of experimental data for temperatures close to the critical point A (e.g., T = 18K) can be summarized as follows: the first-order phase transition between a tilted and homogeneous states manifests itself in a rather unpronounced way; indeed, the spirals smoothly slant starting directly from the conical state and almost return back, which may be addressed by the anisotropy value k _c = 0.1; an experimental ratio of two fields H _c2 for the considered two temperatures is 1.0575, which theoretically results in k _c = 0.12 if one assumes k _c = 0.15 for T = 2.

The value k _c = 0.05 is located on the other side from the critical point A and thus corresponds to experimentally higher temperatures. The value of the exchange anisotropy must have an upper boundary, b _ea > − 0.08, since no oblique spirals in zero field were observed. One may remark that in this case other nucleation mechanisms should be considered to explain SkL stability [6]. Experimentally, the tilted spiral is found to persist up to above T = 35 K, which is not supported by the theoretical results with k _c < 0.04. Still, due to the “flat” energy minimum for the conical state, one observes experimentally some broadening of the conical peaks without a preferable tilt direction, an effect that spans the whole conical phase and persists up to above 35 K [7]. The subsequent endeavor to construct a quantitative model for a bulk helimagnet Cu₂OSeO₃ must take into account experimental data for other field directions. Moreover, non-hysteretic magnetization processes with high experimental accuracy must be given a fitting preference.

In order to narrow down the range of anisotropy coefficients in the phenomenological model of a bulk helimagnet Cu₂OSeO₃ and to build its complete quantitative model, we examined the low-temperature behavior of tilted skyrmion and spiral states. By comparing the results of numerical simulations with the plentiful experimental data from our previous publications [3, 5–7] we were able to ascribe the range of the cubic anisotropy, k _c = 0.15–0.10, to the temperature range T = 2–18K. The course of reasoning leading to a complete quantitative model for a bulk helimagnet Cu₂OSeO₃ is based on simultaneous examination of stable skyrmion and tilted-spiral states. One could not rely on the behavior of the tilted spiral alone, since it would result in some alternative parameter sets inconsistent with the behavior of LT-skyrmions. The complex reorientation processes of metastable skyrmions (Figure 6; Figure 7) are barely reflected at the experimental phase diagrams and are often hidden by the reorientation of spiral states. However, an experimental indication of tilted skyrmion states is presumably manifested by their abrupt disappearance above the H c 1 * line for H‖[110] (see Figure 5B in Ref. [3]).

In the present manuscript, we used definitions of different modulated states (for example, transversal and conical spirals) conventional in the theoretical modeling. In the experiments, however, one would probably not specify different spiral states below the critical field H _c1 and would describe it as a multidomain spiral state. Moreover, the line H _c1 is not shown at the theoretical phase diagram in Figure 1B, since it is related to the reorientation processes of metastable transversal spirals, which, according to the theory, must not be reproduced with the decreasing magnetic field. Therefore, besides the further step to finalize the quantitative model of a bulk helimagnet Cu₂OSeO₃, we pursued the goal to make the terminology used in the theoretical and experimental aproaches consistent.