Stability of the skyrmion lattice near the critical temperature in cubic helimagnets

The phase diagram of cubic helimagnets near the critical temperature is obtained from a Landau-Ginzburg model, including fluctuations to gaussian level. The free energy is evaluated via a saddle point expansion around the local minima of the Landau-Ginzburg functional. The local minima are computed by solving the Euler-Lagrange equations with appropriate boundary conditions, preserving manifestly the full nonlinearity that is characteristic of skyrmion states. It is shown that the fluctuations stabilize the skyrmion lattice in a region of the phase diagram close to the critical temperature, where it becomes the equilibrium state. A comparison of this approach with previous computations performed with a different approach (truncated Fourier expansion of magnetic states) is given.

The phase diagram of cubic helimagnets near the critical temperature is obtained from a Landau-Ginzburg model, including fluctuations to gaussian level. The free energy is evaluated via a saddle point expansion around the local minima of the Landau-Ginzburg functional. The local minima are computed by solving the Euler-Lagrange equations with appropriate boundary conditions, preserving manifestly the full nonlinearity that is characteristic of skyrmion states. It is shown that the fluctuations stabilize the skyrmion lattice in a region of the phase diagram close to the critical temperature, where it becomes the equilibrium state. A comparison of this approach with previous computations performed with a different approach (truncated Fourier expansion of magnetic states) is given.

I. INTRODUCTION
Cubic helimagnets support chiral magnetic textures called skyrmions that are very interesting from the fundamental and applied physics points of view. From the point of view of fundamental physics, skyrmions are solitonic states that emerge as topologically nontrivial solutions of nonlinear field equations 1 . From the point of view of applied physics, they are very promising ingredients for spintronic devices 2,3 , since they are textures modulated at the nanoscale, very robust due to the protection provided by the topology, and appear spontaneously under certain conditions that can be externally controlled via the temperature, the magnetic field, or the electric current.
Theoretically, skyrmions can appear in bulk materials as isolated excitations of the forced ferromagnetic phase (FFM) 4 . However, mean field computations indicate that they cannot condense spontaneously into equilibrium states such as a skyrmion lattice (SKL), since the competing conical state (CS) has lower free energy 1,5 . To stabilize the SKL other ingredients, like high magnetic anisotropy 6 or modifications of the magnetic stiffness 1 are necessary. In two dimensional space, however, the SKL become the equilibrium state in some regions of the phase diagram 7 . Experimentally, SKLs have been found in thin films samples of some materials 8 .
Despite the theoretical considerations, the A-phase that appears in bulk cubic helimagnets has been identified with a SKL [9][10][11] . It has been proposed by Mühlbauer et al. that it is stabilized by fluctuations that modify the mean field computations 9 . This idea is supported by Monte Carlo simulations 12 . Mühlbauer et al. addressed the problem by minimizing the appropriate free energy using a truncated basis of Fourier modes for the magnetic configurations. Adding the contribution of fluctuations, they found as equilibrium solution a magnetic state different from the CS and compatible with a SKL. However, it is unclear how this equilibrium state is related to the magnetic skyrmions of Bogdanov and collaborators 4,13,14 , which have a highly nonlinear structure, formed by a central core that has a magnetic moment pointing opposite to the applied field, surronded by a FFM background in which the magnetic moments point in the diretion of the field. The Monte Carlo simulations, although valuable since they provide some insight, are far from being conclusive, as they are plaged of technical problems: first, the physics of the problem involve two very different scales, the lattice parameters of the underlying crystal and of the SKL, that cannot be accomodated in a Monte Carlo simulation; second, topology divides the configuration space into separated sectors which cause the lost of ergodicity of the Monte Carlo algorithms.
In this paper we explore the idea that fluctuations stabilize the SKL using a different approach, working in configuration space with the fully nonlinear skyrmions of Bogdanov. In this way we clarify the nature of the SKL and get insight into the physics of the problem, showing that the equilibrium skyrmion lattice is a crystal of solitons. This means that the strong chiral fluctuations that characterize the precursor phase of the A-phase [15][16][17] should be described by the nonlinear skyrmions as elementary excitations, rather than by a bunch of weakly coupled Fourier modes.
To close the introduction, it is worth mentioning that very recently a SKL has been found in a cubic helimagnet at low temperature for magnetic field in the range close to the transition to the FFM 18 . The existence of a stable SKL in this part of the phase diagram, stabilized by fluctuations, was predicted in Ref. 19.

II. THEORETICAL FRAMEWORK
Consider a cubic helimagnet whose equilibrium properties are described by the partition function where m is the local order parameter, proportional to the local magnetic moment, and W = q 0 d 3 r W , a Landau-Ginzburg (LG) functional, with (2) The parameter q 0 has the dimensions of inverse length and sets the scale for the spatial modulation of the magnetic configurations. The index i runs over {x, y, z}, summation over repeated indices is understood and ∂ i stands for ∂/∂x i . The first and second terms in (2) correspond to the ferromagnetic and Dzyaloshinkii-Moriya exchange interactions, respectively, and the last term within the brakets is the Zeeeman energy, with h proportional to the applied magnetic field, which we take along the zaxis: h = hẑ. The quadratic and quartic terms in m are standard in LG models: a controls the temperature variation, and we fix the coefficient of the quartic term to 1/4. This can always be done by a suitable rescaling of m and a redefinition of h and β. Alternatively, by a rescaling of m we could set β = 1. In this case a different coefficient of the quartic term appears and we would have a more standard form of the LG theory. To study the fluctuations, however, we prefer the equivalent form presented here.
It is convenient to work with the dimensionless free energy density where V is the volume. To evaluate the partition function we assume that it is dominated by the local minima of W, which is a good approximation if β is large. The local minima are solutions of the Euler-Lagrange (EL) equations, δW/δ m = 0, that explicitely read (4) Let us denote a generic local minimum by m 0 and expand W in powers of ξ = m − m 0 up to quadratic order: where the indices α and β run over {x, y, z} and is a symmetric differential operator that is define positive since m 0 is a local minimum of W. The symbol ǫ αβγ stands for the totally antisymmetric tensor. The linear term in (5) vanishes on account of the EL equations. The partition function gets a contribution from each local minimum, obtained as an integral over ξ that can be readily evaluated since the integrand is gaussian. The free energy density associated to m 0 has the form f = f 0 + f 1 , where is sometimes called the mean field or classical contribution, and is the contribution of the fluctuations to gaussian order. The operator K 0αβ = −∇ 2 δ αβ , which does not depend on m 0 , is introduced merely as a convenient form to normalize the term coming from fluctuations. In the continuum limit f 1 is divergent and a short distance cut-off, Λ, has to be introduced. The crystal lattice, whose latttice parameter is denoted by a L , provides a natural cut-off, Λ = π/a L . To compute the total free energy one has to solve the spectral problem K αβ ξ β = λξ α . However, away from criticality, f 1 is dominated by the short distance fluctuations. That is, ln det(KK −1 0 ) is dominated by the largest eigenvalues of K, provided it does not have infrarred divergences, as it is the case if we are outside the critical region. The largest eigenvalues of K correspond to the short distance modes, and for these K ∼ K 0 , so that we have The last approximation in the above expression relays on the dominance of short distance modes, for which 0 ] using the basis of eigenstates of K 0 given by plane waves ξ α = e i k rû α , with the three polarization vectors satisfyingû α ·û β = δ αβ . We obtain where Expression (10) is general for any local minimum, m 0 , of the LG functional. Notice that, in spite of its appearance, f 1 is not merely a renormalization of the LG m 2 term that can be absorbed in a redefinition of a. The local minimum m 0 has to be computed with the LG functional, with the parameter a multiplying the m 2 term. Then, the gaussian fluctuation around this local minimum give the f 1 contribution to its free energy, which depends on a through m 0 .
The saddle point expansion is reliable if the coeffcient c 0 is small, that is, if β is large. For the cubic helimagnet MnSi q 0 ≈ 0.03Å −1 and a L ≈ 4.5Å, so that Λ/q 0 ≈ 20. Hence we may take c 0 = 2.5/β. The value of β is unknown, but it has to be large to justify the saddle point expansion.
The validity of the saddle point expansion and of the short distance approximation is limited by the behavior of the soft modes. Usually these become important only in a neighborhood of the critical region, where they produce infrared divergences that invalidate the saddle point expansion and, of course, the short distance approximation. The fact that the transitions are of first order 16,17 means that it is possible that the soft modes play little role in the chiral states, so that that most of the phase diagram obtained in the gaussian approximation is qualitatively valid. In any case, within the present approach it is not possible to delimite the region of validity of the computations, given that we do not consider the soft modes. In addition, our treament of the CS, presented in Sec. III, is not valid for very low values of the magnetic field. In this case the CS is nearly degenerate, as its wave vector can be rotated away from the magnetic field direction with a small energy cost. This situation requires a modification of mean field theory known as Brazovskii theory 21,22 .

III. CONICAL STATE
The CS is a solution of the EL equations of the form m 0 = (m t cos qz, m t sin qz, m z ), where m t , m z , and q are constants. The EL equations are satisfied if and only if where ∆ = (q − q 0 )/q 0 . Notice that ∆ 2 is limited by the inequality m 2 t ≥ 0, what implies that ∆ 2 ≤ ∆ 2 max , where ∆ 2 max satisfy the equation It is not difficult to see that ∆ 2 max ≤ 1 − a and that ∆ max = 0 for a = 1 − h 2 . For a > 1 − h 2 Eq. (13) gives m t < 0 for any ∆ 2 . Thus, the CS exists only for a < 1 − h 2 .
The free energy density including gaussian fluctuations is The last term, proportional to c 0 , is the contribution of the fluctuations in the short distance approximation, given by Eq. (10). The equilibrium value of ∆ is obtained by minimizing the free energy. In absence of fluctuations, c 0 = 0, the minimum of the free energy is attained at ∆ = 0, that is, q = q 0 , for any value of a and h. A second order phase transition to the FFM state takes place when a = 1−h 2 . This is the mean field critical line. For small c 0 the free energy minimum is still attained at ∆ = 0. This minimum disappears on the critical line a = 1 − h 2 − 2c 0 . We see that, as expected, the fluctuations lower the critical temperature.
Using the parametrization of the local order parameter in polar coordinates m = m (sin θ cos ψ, sin θ sin ψ, cos θ) , and cylindric coordinates (r, ϕ, z) for the spatial points, the skyrmion is a solution in which ψ = ϕ + π/2 and m and θ are functions of the radial coordinate, r, only, so that it is invariant under rotations around the magnetic field axis,ẑ. The EL can be cast to the form where the prime stands for the derivative with respect to r. The boundary conditions are where the non negative constant m s is the value of the order parameter in the FFM homogeneous state, since as r → ∞ Eq. (17) approaches the EL equation for the homogeneous state. The condition m ′ (0) = 0 is necessary to have a finite solution at r = 0. Thus, the skyrmion has no degree of freedom. Since θ and m tend exponentially to 0 and m s as r → ∞, the skyrmion is a localized structure that consists of a central core surronded by a FM background 23 . The core size is controlled by the external parameters a and h.
It is believed that A-phase that appears in several cubic helimagnet is formed by the condensation of skyrmions in a lattice with hexagonal cells that have a skyrmion core in its center. The lattice breaks the skyrmion rotational symmetry. However, if the lattice cell is larger than the skyrmion core, we can approximate the lattice cell by a central axisymmetric skyrmion core surronded by a deformed FM background. Then, the energy of the skyrmion lattice can be computed in the so called circular cell approximation (CCA) 14 , in which the skyrmion profile is a solution of Eqs. (17) and (18)  The skyrmion lattice (SKL) has thus two free parameters, R and m s , that are fixed by minimizing the free energy. Notice that in this case m s is not forced to be equal to the value of the order parameter in the homogeneous state.
At mean field level, the SKL has higher free energy than the CS for any value of a and h; thus it is at most metastable. To obtain a stable SKL one has to either modify the model 1,6 or go beyond the mean field approximation and include the effect of fluctuations at gaussian level 9 . Let us analyze the last possibility. Fig. 1 displays the difference of f 1 /c 0 , given by Eq. (10), between the SKL and the CS, computed with the mean field equilibrium solutions. As this difference is negative except in a very narrow interval of a in the vicinity of the transition point, we see that the free energy of fluctuations contributes to lower the free energy of the SKL with respect to the CS. Thus, thermal fluctuations are a potential mechanism to stabilize the SKL.
To explore this possibility, we solve numerically the boundary value problem defined by Eqs. (17) and (18) and the boundary conditions (20). Then, we determine the equilibrium values of m s and R by minimizing the free energy that includes the fluctuations, f = f 0 + f 1 . We fix c 0 = 0.01. With such small value of c 0 , the equilibrium values of m s and R are only slightly shifted from their mean field values. A typical result is shown in Fig. 2, which displays the components of the free energy as a function of q 0 R. In these plots the value of m s has been fixed by minimizing the total free energy, keeping R fixed. The blue line is the difference of total free energies between the SKL and the CS. The red and green lines repre- sent respectively the mean field (f 0 ) and fluctuation (f 1 ) components of the this free energy difference. The SKL is more stable than the CS since its total free energy is lower in a neighbourhood of the minimum. Thus, in this case (h = 0.15 and a = 0.9), the SKL is the equilibrium state. The equilibrim skyrmion profile in the unit cell is displayed in Fig. (3). ∼ h < ∼ 0.167, it becomes the equilibrium state in an interval of a; below and above this interval it is metastable or unstable. c) For 0.155 < ∼ h < ∼ 0.162 the SKL is the equilibrium state in two disjoint intervals of a: by increasing a from negative values, where the SKL is metastable, it becomes the equilibrium state at a certain a; at a higher a the SKL loses stability and becomes again metastable; and it regains stability at a still higher value of a, and remains the equilibrium state until it disappears. d) For h < ∼ 0.155, the SKL becomes the equilibrium state at high enough a and remains so until its disappearance. The phase diagram is reconstructed from this results. Let us discuss it in the following paragraphs.

V. PHASE DIAGRAM AND DISCUSSION
At mean field level (c 0 = 0), the phase diagram is displayed in Fig. 5. A second order phase transition takes place on the red line, which separates the homogeneous FFM phase (green region) from the conical phase (red region). The SKL is metastable within the conical phase, in the region with blue stripes. The local minimum of the free energy that defines the SKL disappears on the dashed blue line. In contrast with what happens in models where the modulus of the order parameter is fixed, the SKL does not disappear through a nucleation process, in which the lattice size R diverges and the homogeneous state is attained smoothly 6,14 . In the present case, with a soft modulus, the size of the SKL remains bounded and the metastable state disappears because the minimum depth gets gradually shallower until the minimum becomes a inflection point. Figure 6 displays the phase diagram including fluctuations, with c 0 = 0.01. The effect of the fluctuations is quantitative, shifting the line boundary between the homogeneous and CS to lower values of a (lower temperatures), and qualitative, stabilizing the SKL in the blue region. On the phase boundary (blue line), a first order transition takes place, since the SKL and the CS cannot be smoothly conected. With higher values of c 0 the phase diagram is similar, with the SKL stable phase extended to lower values of a. For instance, with c 0 = 0.1 the region where the SKL is stable has a similar form to the blue region of Fig. 6 (right), but it covers the interval from a = −0.4 to a = 0.75 and from h = 0 to h = 0.5. In this case, the transition from the FFM to the CS takes place a = 0.8 for h = 0.
The results may be invalid near the phase boundary in the high a region, where soft modes may become important, causing the failure of both the saddle point expansion and the short distance approximation. Also, the high (approximate) degeneracy of the CS for low field values invalidates the computation in the low field region, where the CS has to be treated with a mean field theory of Brazovskii type.
The phase diagram is very similar to that obtained by Mühlbauer et al. 9 using a a truncated Fourier basis to obtain the mean field local minima and to minimize the total free energy. These authors identify the new equilib-rium state stabilized by the fluctuations with a SKL. Although this identification is reasonable, it is unclear how the spin texture determined through the summation of a finite number of Fourier modes is related to the highly nonlinear skyrmion configurations of Bogdanov. We have shown here that results similar to those of Mühlbauer et al. can be obtained by working directly with the fully nonlinear skyrmion texture in configuration space: a lattice formed by the condensation of skyrmions is the equilibrium state at low enough field near the critical temperature. It can be identified with the A-phase of cubic helimagnets. These results shed light on the internal structure of the SKL and thus provide more insight to the physical aspects of the problem. Of special importance is the fact that the strong chiral fluctuations that characterize the precursor phase of the A-phase 16 should be described by nonlinear skyrmion tubes as elementary excitations, rather than by the dynamics of a bunch of weakly coupled Fourier modes.