Discontinuous to continuous transition changeover and magnetic helicity reversal in helimagnet nanodisks under torsion

A magnetic skyrmion is topologically protected because it possesses a non-zero topological charge. As a result, the creation or annihilation of a magnetic skyrmion is thought to be initiated by a sudden reversal of local magnetization, and thus cannot occur continuously. Here, we show that this viewpoint is only partially correct by studying the creation and annihilation of an isolated skyrmion in a nanodisk suffering coupled magnetic field and mechanical torsion. It was found that at a proper magnetic field, the torsion-induced chiral stress field can change the helicity of the magnetization structure and create or annihilate isolated skyrmion in a continuous way. In the torsion-magnetic-field phase diagram, there appears a critical point, where the type of the topological transition changes from a discontinuous to a continuous one. It was further proved that the critical condition for the continuous transition is the softness of magnetization at the center.

When isolated skyrmions are regarded as information memory carrying the information '1' , it is of great significance to understand the process of their creation and annihilation from or to the topologically trivial state carrying the information '0' .However, current models only consider the low-temperature case with constant magnetization modulus [41], or they ignore the influence of geometric boundaries [42].Whether the transition between the '0' and '1' states is discontinuous, as in [41], or continuous, as in [42], is to be explored, when taking into consideration the softness of the magnetization amplitude induced by temperature, as well as the geometric boundaries of the helimagnets.
Magnetic skyrmions possess internal structure and can be classified into Bloch skyrmion, Néel skyrmion, Bloch-Néel intermediate skyrmion, and antiskyrmion by two features: the vorticity (distinguishing skyrmion and antiskyrmion), and the helicity (distinguishing Bloch skyrmion, Néel skyrmion and intermediate skyrmion) [1,[43][44][45].This means that skyrmions have the potential to carry more information than a simple binary code of '1' [46].However, current research has primarily focused on controlling the number of skyrmions in nanostructures by applying various external fields [41,47,48], without considering the information that their internal structure may carry.This is due to the fact that in chiral magnets, the Dzyaloshinskii-Moriya interaction (DMI) has fixed the vorticity and helicity [43,49], making it difficult to alter them.We notice a recent experiment in which a so-called target skyrmion, composed of a whole isolated skyrmion at the center and the concentric helical stripes at the edge with opposite helicities, is stabilized in a helimagnet nanodisk [50].We may ask a question: is it possible to change the helicity of the edge part with an external field?
In this work, we apply the torsion which can induce a stress field, and study its effects on the magnetization structure in a B20 helimagnet nanodisk.This torsional stress field has two main features: first, it has either count-clockwise (CCW) helicity or clockwise (CW) helicity, similar to the field of magnetization; second, it is inhomogeneous, weak at the center but strong at the edge, thus, it may have obvious effects on the magnetization structure at the edge.By using the Monte Carlo simulations with natural boundary condition [51,52], we get the equilibrium state in the circular nanodisk at different magnetic fields and torsions.We find that a negative torsion stabilizes two skyrmion-related states, one has CCW helicity at the center and CW helicity at the edge, and the other one has CCW helicity, which is actually the fractional skyrmion [53].On the contrary, a positive torsion stabilizes two quasi-uniform states.One has CCW helicity at the center and CW helicity at the edge, and the other has CW helicity.We also plot the torsion-magnetic-field phase diagram and identify a critical point [54,55] where the transition type changes from discontinuous to continuous.Finally, by using the asymptotic analysis, we reveal the features of the continuous transition and give a qualitative explanation of the phase diagram.

Model
By symmetry analysis, we obtain the rescaled Gibbs free energy density g for a B20 helimagnet nanodisk under torsion as (see appendices A and B) Here, m = (m 1 , m 2 , m 3 ) represents the magnetization, x i (i = 1, 2, 3) are the spatial coordinates, b, t and τ are the rescaled magnetic field, temperature, and torsion, respectively.The last term in equation ( 1) represents the torsion-related magnetoelastic coupling.Similar to the torsional stress field, the strength of the magnetoelastic coupling increases with the distance from the center.In order that the nontrivial state in the nanodisk is the isolated skyrmion, we fix the radius of the nanodisk at an appropriate rescaled value R = 4.6.Since, if R is too small, skyrmion cannot be stabilized, and if R is too big, skyrmion cluster will appear [41].To simplify the study, we assume a fixed temperature t = 0.5, which is near the ordering temperature [56,57], and ignore its effects on the magnetization structure.At certain magnetic field b and torsion τ , we calculate the equilibrium state by applying the Monte Carlo simulation which finds the minimum of the Gibbs free energy of the system.Due to the Landau terms in equation (1), the spatial longitudinal modulations of the magnetization modulus should be taken into account in the simulation, i.e. m = |m| is not a constant and varies with spatial coordinates [57,58].The boundary condition for the simulation is the natural boundary condition [51,52], and is automatically satisfied if the Monte Carlo simulation finds a minimum of the free energy of the system.

Helicity reversal induced by torsion
We first consider a simple case without any torsion and calculate the magnetization structure in the nanodisk under different magnetic fields.The results are as expected: at a lower magnetic field, an isolated skyrmion is stabilized, while at a higher magnetic field, a trivial state without topological protection is stabilized (see figures 1(a) and (d)).The trivial state is characterized by positive out-of-plane magnetization components and a vortex-like, in-plane structure, and we call it the quasi-uniform state.At the edge of the nanodisk, the variational calculus results in the following natural boundary condition where n is the unit vector normal to the boundary of the nanodisk, and e z is a unit vector along the z-axis.Due to the circular shape of the nanodisk, the DMI-induced term (e z × n) × m of equation ( 2) is nonzero at the edge, causing dm dn ̸ = 0 and favoring the quasi-uniform state instead of the uniform state where all magnetization vectors point in the same direction.The value of temperature is fixed at t = 0.5.The in-plane and out-of-plane components of magnetization are represented by the arrows and the contour, respectively.The part inside the white dotted circle has CCW helicity, and the part outside the white dotted circle has CW helicity.
We then apply the torsion and investigate its influence on the magnetization structure in nanodisks.Our findings reveal that a positive torsion favors a quasi-uniform state.Specifically, At b = 0.3, introducing a positive τ = 0.1 alters the equilibrium state from an isolated skyrmion to a quasi-uniform state (figure 1(b)).We further increase the value of torsion to τ = 0.2, causing the disappearance of the white dotted circle (figure 1(c)), within which the in-plane magnetization swirling direction is CCW.This means that the helicity of the quasi-uniform state changes from CCW/CW to CW. Next, we examine the effects of negative torsions.Comparing figures 1(a) and (e), we observe that a negative torsion leads to an increase in the radius of the white dotted circle, resulting in a larger isolated skyrmion.
Given that the Gibbs free energy density terms and nanodisk geometry are all axially symmetric, the equilibrium magnetization structures of the nanodisk exhibit axial symmetry (see figure 1) and can be described via two functions m(ρ) and θ(ρ).In this context, m = |m| is the magnetization modulus, θ is the polar angle between the magnetization vector and z-axis, and ρ denotes the position.As a consequence, the magnetoelastic coupling in equation ( 1) can be reformulated as: From equation (3), we can discern that a positive torsion τ prefers a negative θ, while a negative torsion τ prefers a positive θ.This explains the observation in figure 1 that the region with θ > 0, representing CCW helicity, is enlarged (narrowed) by a negative (positive) torsion.
Inspired by the observation that negative torsion augments the size of an isolated skyrmion in a nanodisk, we investigate the effects of even stronger negative torsion on the behavior of skyrmions.We extract information about θ(ρ) and m(ρ) from the equilibrium magnetization structure obtained by Monte Carlo simulations, and plot θ(ρ) and m(ρ) in figure 2. Our results reveal noticeable changes in the magnetization structure near the edge of the nanodisk as the strength of negative torsion is increased from −0.2 to −0.8 at b = 0.3, while the changes near the center remain modest.These results are reasonable due to the fact that the magnetoelastic coupling is proportional to the distance from the center (see equation ( 3)).Notably, the polar angle θ(R) of the magnetization vector at the edge of the nanodisk increases from negative to zero and then to positive as the negative torsion strength increases.The inset of figure 2 shows k = (θ(0) − θ(R))/π as a function of the torsion.For τ > −0.5, we have k > 1, meaning that the magnetization structure is composed of a whole skyrmion with CCW helicity at the center and the concentric helical stripes with CW helicity at the edge [50].Conversely, for τ < −0.5, we have k < 1, indicating that the section with CW helicity disappears and the magnetization structure does not contain a whole skyrmion, and we call it a fractional skyrmion (see also figure 1(f)).
Under certain conditions of b and τ , radial profiles of θ(ρ) and m(ρ) can also be obtained by numerically solving the nonlinear Euler equations derived by the variational calculus (see appendix B).The orange curves in figure 2 show the profiles obtained using this numerical method, and they exhibit good agreement with the Monte Carlo simulation results, further validating the reliability of our simulation approach.

τ -b phase diagram and continuous transition
In our study, we carry out the Monte Carlo simulations with the initial state set as either the uniform state or the isolated skyrmion state.Interestingly, we observe that in certain cases, the final magnetization state depends on the initial state, while in other cases, it does not.This peculiar behavior motivates us to investigate the effects of torsion in greater detail.To begin, we fix the magnetic field at 0.35, set the uniform state as the initial state, and decrease the torsion from positive to negative.We find that the transition from quasi-uniform state to skyrmion occurs at about τ 1 = −0.016.Then, we set the isolated skyrmion as the initial state and increase the torsion from negative to positive.We find that the reverse transition from skyrmion to quasi-uniform state occurs at another value of torsion about τ 2 = 0.011.We choose the out-of-plane magnetization component at the center of the nanodisk m z0 (whose negative/positive value corresponds to the skyrmion/quasi-uniform state), and the average out-of-plane magnetization over the nanodisk mz , as the y coordinates, and plot the hysteresis loop in figure 3(a).Due to the abrupt change of m z0 and mz , the transition is referred to as discontinuous.The magnetization distributions at b = 0.35 and certain torsions are plotted in figure 3(d).Increasing τ from −0.03 to 0 results in the skyrmion state, while decreasing τ from 0.03 to 0 results in the quasi-uniform state.We repeat the process under another magnetic field.A similar hysteresis loop is obtained at b = 0.4 as shown in figure 3(b).At a higher magnetic field, the torsion region ∆τ = τ 2 − τ 1 , where skyrmion or quasi-uniform state can exist as a metastable state, is narrower.We expect that with a further increase in magnetic field, the region width will approach zero.This prediction is confirmed, as we have set b = 0.5 and found that m z0 and mz vary continuously with torsion (see figure 3(c)), indicating a continuous transition between the skyrmion and quasi-uniform state.Figure 3(e) shows the magnetization distributions at b = 0.5.At τ = −0.165,there appears a special state whose magnetization at the center is near zero.It is the existence of this special state that makes the continuous transition between the skyrmion and quasi-uniform state possible.We search for the transition points at different magnetic fields and plot the τ -b phase diagram in figure 3(f).We observe that as the magnetic field increases, the transition transforms from a discontinuous transition to a continuous transition.The discontinuous to continuous transition changeover indicates that the topological protection of skyrmion is disrupted by the softening of the magnetization modulus.
Next, we investigate the behavior of isolated skyrmions near the continuous transition.We select a relatively high magnetic field b = 0.5 and calculate, by Monte Carlo simulation, θ(ρ) and m(ρ) at various torsions ranging from −2 to −0.164, a point close to the continuous transition τ pt = −0.162.As shown in figure 4, the magnetization structure close to the center undergoes drastic changes with varying torsion, whereas the structure near the edge remains mostly unaffected.This is in contrast to the observations in figure 2. These findings suggest that near the continuous transition point, the central part of the nanodisk is sensitive to external fields.To further illustrate this, we introduce the concept of emergent elastic stiffness  coefficients, which characterize the ability of the magnetization structure to resist its own deformation (see appendix D).We hope that some of these emergent elastic stiffness coefficients near ρ = 0 approach zero when the continuous transition occurs.C e ρρ = 2 seems to be a candidate.At ρ = 0, the boundary conditions equations (B.4b) tell us m ′ (0) = 0, thus C e ρρ (0) = 2m 2 (0)θ ′2 (0).We know from figure 4(b) that lim τ →τpt m(0) = 0. Nevertheless, the inset of figure 4(a) indicates that the closer to the transition point, the faster the size of the skyrmion core, with π/2 < θ < π, decreases.In other words, lim τ →τpt θ ′ (0) = −∞.Consequently, the convergence of lim τ →τpt C e ρρ (0) cannot be determined directly.We, therefore, resort to the asymptotic analysis near ρ = 0 (see appendix C).It is found that only by setting θ(0) = π/2 in equation (C.1), we can obtain m(0) = 0. To fulfill the conditions lim τ →τpt m(0) = 0 and lim τ →τpt θ ′ (0) = −∞, it requires that θ jumps to a finite value π/2 and then decreases smoothly with ρ increasing from zero when the continuous transition occurs.Thus, we have lim During the continuous transition from skyrmion to quasi-uniform state, the magnetization structure at the center is softened.On the one hand, the magnetization modulus declines from a finite value to zero; on the other hand, it becomes more sensitive to external fields or even small perturbations.

Qualitative explanation of the phase diagram
Finally, we offer a qualitative explanation of the phase diagram.Let X = θ ′ (0) + 1 and Y = m(0), then P 1 = 0 in equation (C.3) can be rewritten as For fixed b, t and τ , the left-hand-side of equation ( 5) can be regarded as a cubic function of X without quadratic term, which is represented by the blue curves in figure 5; the right-hand-side of equation ( 5) can be regarded as a constant represented by the red dotted lines in figure 5.If b is large, the coefficient of the linear term tends to be positive, and the equation can have only one solution as shown in figure 5(b).In this case, increasing τ leads the red line to move upward, and a continuous transition occurs from a negative θ ′ (0) (corresponding to skyrmion) to a positive θ ′ (0) (corresponding to quasi-uniform state).It should be emphasized that when τ or b changes, the shape of the cubic function changes accordingly.If b is small, the coefficient of the linear term tends to be negative, the equation may allow multiple solutions as shown in figure 5(a).In this case, a discontinuous transition will happen if the red line passes the stationary points (with zero slopes) of the cubic function of X.

Summary
To summarize, we have investigated the effects of torsion on magnetization structures in helimagnet nanodisk, including the skyrmion state and quasi-uniform state.We find that at a certain temperature and magnetic field, the torsion can induce four states with different helicities.We also find that the transition from skyrmion to quasi-uniform state changes from discontinuous to continuous under certain conditions of magnetic field and torsion.Moreover, we prove that the continuous transition is achieved by softening the magnetization structure at the center.
Here, M is the magnetization vector and X = (X 1 , X 2 , X 3 ) is the position vector.The first three terms in equation (A.1) denote the Heisenberg exchange interaction with stiffness A, the DMI with coefficient D, and the Zeeman energy under magnetic field B, respectively.ω L = α(T − T 0 )M 2 + βM 4 contains the second and fourth Landau expansion terms with the expansion coefficients α and β, T represents the temperature, and T 0 is the ordering temperature of a ferromagnet.ω el is the elastic energy density, it has the following form for B20 Materials where ε ij are the strains, C ij are the elastic coefficients.ω el is the magnetoelastic energy, it can be written as a most general form including fourth-order and fifth-order tensors where λ ijkl and γ ijklm are the coefficients of the fourth-order and fifth-order tensors, respectively; and the Einstein summation convention is applied.According to the transformation of tensors, we have, where R is the transformation matrix, Λ ijkl and Γ ijklm are the tensor coefficients after transformation.If R corresponds to a symmetric operation of the system, i.e.R is a element of the point group of the system, we have Considering the constrains of equation (A.5) and of Lifshitz invariants, which require γ ijklm = −γ ijlkm , we can finally obtain independent coefficients λ ijkl and γ ijklm and derive the expression of ω el for B20 materials as where λ n and γ n are the independent magnetoelastic coefficients, In this work, the stress boundary conditions are chosen; therefore, the appropriate thermodynamic potential is the Gibbs free energy density.Apply the Legendre transformation, we get the Gibbs free energy density as New J. Phys.26 (2024) 023009 and S = C −1 the flexibility matrix.Here, in equation (A.9), we have ignored the terms − 1 2 F T SF (due to its smallness) and −ω el (σ).With the following rescaling parameters the Gibbs free energy density can be expressed as We consider the torsion-induced shear stresses as the boundary conditions and ignore the terms that are not axially symmetric for simplicity.In this case, we have σ 11 = σ 22 = σ 33 = σ 12 = 0, and the magnetoelastic coupling term can be simplified as The torsional shear stresses σ 13 and σ 23 are proportional to the distance from the center, where R is the rescaled radius of the nanodisk, ρ is the distance from the center, ϕ is the polar angle, and T is the applied torsion.Equation (A.13) can be finally derived as

Appendix B. Numerical method of solving nonlinear Euler equations
The system we study possesses axial symmetry, and the circular cell approximation can be used to describe the magnetization structure of the nanodisk.We introduce the cylindrical coordinates for the position vector x = (ρ cos ϕ, ρ sin ϕ, z) and the spherical coordinates for the magnetization vector m = m(sin θ cos ψ, sin θ sin ψ, cos θ).Substitute x and m into equation (A.15), we have  To solve the Euler equations, we use the iteration method [61].As a first step, we linearize equations (B.3) near the approximation solutions m a and θ a , which are chosen appropriately, and get the linear differential equations of the corrections m c and θ where where X represents the corrections θ c and m c , A is the coefficient matrix, and B is the constant term.When an eigenvalue of A approaches zero, X = A −1 B will be large, leading to the divergence of iteration.In other words, the transition between different magnetic states occurs when A has a zero eigenvalue.We should also emphasize that some dynamic information can be extracted from the matrix A. The dynamic Landau-Lifshitz equation is given by where B eff , the effective magnetic field, is defined as the functional derivative B eff = δg δm .To analyze the stability of a given state of the nonlinear system controlled by the Landau-Lifshitz equation, we should linearize and then discretize equation (B.9) [62].We can finally obtain a similar coefficient matrix related to A, whose eigenvalues determine the stability of the a given magnetic state.Both the convergence of the iteration method and the stability of magnetic states of the dynamic system dependent on the eigenvalues of A. We can also conclude that a soft mode will appear during the magnetic transition, i.e. the frequency of a collective mode will be zero.

Appendix C. Asymptotic analysis
To get more information on the magnetization near the center of the nanodisk, we approximate θ and m by their Maclaurin series: where θ(0) = π for skyrmion (or θ(0) = 0 for quasi-uniform state).Substitute the Maclaurin series (C.1) into the Euler equation (B.3), we can get the following polynomial equations: which require the coefficients P i (i = 0, 1, 2, • • • ) and Q i (i = 0, 1, 2, • • • ) being zero.We then have (D.5) Obviously, they are non-negative at any point.When the eigenvalues are large (small), it is relatively hard (easy) to deform the magnetization structure.For this reason, we call them 'emergent elastic stiffness coefficients' .

Figure 2 .
Figure 2. (a) The polar angle θ between the magnetization vector and the z-axis, and (b) the modulus m of the magnetization vector, as functions of the position ρ.The dots and orange curves are obtained by Monte Carlo simulations and by solving the Euler equation (B.3), respectively.The results are calculated at b = 0.3, t = 0.5.The inset shows the torsion dependence of k in the nanodisk.

Figure 3 .
Figure 3. mz0, the out-of-plane component of magnetization vector at the center of the nanodisk, and mz, the average mz over the nanodisk, as functions of torsion τ at different magnetic fields (a) b = 0.35, (b) b = 0.4, and (c) b = 0.5.Magnetization distributions at certain conditions are plotted in (d) and (e).(f) τ -b phase diagram.The red and blue curves represent the discontinuous and continuous transitions, respectively.
Then, we solve equation (B.6) numerically and let m a + m c and θ a + θ c as the new approximation solutions.Repeating the steps above, we can finally get the numerical solution m and θ of equation (B.3) at certain b, t 12 equation (D.1), ε 11 , ε 22 , ε12and ω represent the deformation and rotation of the magnetization structure, respectively, and we call them the 'emergent strains' .According to equations (D.1) and (D.2),We have Then, substitute equations (D.1)-(D.3)into the Gibbs free energy density, one obtains a function of emergent strains g = g (ε e 11 , ε e 22 , ε e 12 , ω e ) .(D.4)If the solution of the magnetization is stable, the Hessian matrix H(g) of g(ε 11 , ε 22 , ε 12 , ω) should be positive definite, i.e. the eigenvalues of H(g) are non-negative.In cylindrical coordinate system, the four eigenvalues of H(g) can be derived as