Incommensurate--commensurate transitions in the mono-axial chiral helimagnet driven by the magnetic field

The zero temperature phase diagram of the mono-axial chiral helimagnet in the magnetic field plane formed by the components parallel and perpendicular to the helical axis is thoroughly analyzed. The nature of the transition to the commensurate state depends on the angle between the field and the helical axis. For field directions close to the directions parallel or perpendicular to the helical axis the transition is continuous, while for intermediate angles the transition is discontinuous and the incommensurate and commensurate states coexist on the transition line. The continuous and discontinuous transition lines are separated by two tricritical points with specific singular behaviour. The location of the continuous and discontinuous lines and of the tricritical points depend strongly on the easy-plane anisotropy, the effect of which is analyzed. For large anisotropy the conical approximation locates the transition line very accurately, although it does not predict the continuous transitions nor the tricitical behaviour. It is shown that for large anisotropy, as in CrNb3S6, the form of the transition line is universal, that is, independent of the sample, and obeys a simple equation. The position of the tricritical points, which is not universal, is theoretically estimated for a sample of CrNb3S6

The zero temperature phase diagram of the mono-axial chiral helimagnet in the magnetic field plane formed by the components parallel and perpendicular to the helical axis is thoroughly analyzed. The nature of the transition to the commensurate state depends on the angle between the field and the helical axis. For field directions close to the directions parallel or perpendicular to the helical axis the transition is continuous, while for intermediate angles the transition is discontinuous and the incommensurate and commensurate states coexist on the transition line. The continuous and discontinuous transition lines are separated by two tricritical points with specific singular behaviour. The location of the continuous and discontinuous lines and of the tricritical points depend strongly on the easy-plane anisotropy, the effect of which is analyzed. For large anisotropy the conical approximation locates the transition line very accurately, although it does not predict the continuous transitions nor the tricitical behaviour. It is shown that for large anisotropy, as in CrNb3S6, the form of the transition line is universal, that is, independent of the sample, and obeys a simple equation. The position of the tricritical points, which is not universal, is theoretically estimated for a sample of CrNb3S6.

I. INTRODUCTION
Chiral magnets are very promising ingredients for spintronic based devices since they support peculiar magnetic textures that affect the charge and spin transport properties in different ways 1 . As these magnetic textures can be deeply altered by magnetic fields, the transport properties can be magnetically controlled, and thus they might be used as functional components of magnetic devices 2,3 .
For example, in the mono-axial chiral helimagnet, an instance of which is CrNb 3 S 6 , the competition between the ferromagnetic (FM) and Dzyaloshinskii-Moriya (DM) interactions leads to an incommensurate magnetic helix propagating with period L 0 along a crystallographic axis, which will be called here the DM axis. The effect of the magnetic field on the helical magnetic ordering is an old topic and was discussed for instance in Ref. 4. In the case of mono-axial DM interactions, application of a magnetic field deforms the helix in a way which depends on the field direction. If the field is parallel to the DM axis, a conical helix is formed, with the spins tilted by a constant angle towards the DM axis, while revolving round it with a period L 0 independent of the magnetic field. If, on the other hand, the magnetic field is perpendicular to the DM axis, the spins remain perpendicular to the DM axis but tend to be aligned with the magnetic field. The spins rotate slowly about the DM axis in the regions where they are nearly parallel to the magnetic field, and the rotation becomes faster as the spin direction separates from the field direction. As the field increases, the regions where the spins ro-tate slowly become very wide and the regions where the spins rotate fast become very narrow, and a Chiral Soliton Lattice (CSL) is formed 1,[5][6][7] . This CSL, which is realised in CrNb 3 S 6 8 , supports dynamic modes like coherent sliding 9 and gives rise to exciting phenomena as spin motive forces 10 and tuneable magnetoresistence [11][12][13][14] .
If the applied field is strong enough, the incommensurate chiral helicoid is transformed in a commensurate FM state, called the forced ferromagnetic (FFM) state. Chiral and incommensurate (IC) to commensurate (C) transitions, which are universal and appear in many branches of physics, take place in these systems. The nature of the IC-C transition and its driving mechanism depend on the angle α between the field and the DM axis. If a field parallel to the DM axis increases, the cone on which the spins lie becomes narrower, but the helix period L 0 remains constant. At a certain critical field the cone closes completely and the FFM state is reached continuously. On the other hand, the effect of increasing a field perpendicular to the DM axis is to increase the period of the CSL, so that the FFM state is reached continuously at a critical field at which the period L of the CSL diverges 1 . Thus, in the two limiting cases of magnetic field parallel and perpendicular to the DM axis the FFM state is reached via two very different mechanisms. On the other hand, even if these systems have been extensively studied, the nature of the transition to the FFM has not been elucidated yet when α is continuously varied between the limiting cases 0 and π/2.
A theoretical attempt has been made recently based on an approximation which makes the model equations easily solvable 15 . As we shall show in section VI the approximation gives the transition line very accurately when the easy-plane anisotropy is large, as in the case of CrNb 3 S 6 . However, the nature of the transition, and in particular its continuity or discontinuity, is not well described within this approximation. In this work we address the important question of the phase diagram of the mono-axial chiral helimagnet in the magnetic field plane by solving the model without any uncontrolled approximation, and the results which come out are in some respects rather different from those of reference 15.
The paper is organized as follows. In section II the model is described and the method of solution is outlined. Sec. III is devoted to analyze the results for the phase diagram. In Sec. IV the nature of the singularities that appear on the transition line is studied. The effect of easy-plane anisotropy is discussed in Sec. V. A comparison with approximate solutions, namely the conical approximation of Ref. 15 and the decoupling approximation developed in this work, is presented in Sec. VI, where it is also shown that the form of the phase transition line is universal when the magnitude of the anisotropy is large. Sec. VII deals with the application of the results to CrNb 3 S 6 . And some conclusions are gathered in Sec. VIII. Finally, details on the numerical methods are given in the appendix.

II. MODEL AND METHOD OF SOLUTION
To settle the question raised in the introduction, we studied the zero temperature ground state of a classical spin system with FM isotropic exchange, mono-axial DM interactions, and single-ion easy-plane anisotropy, in the presence of an external magnetic field.
At zero temperature the spins lying on the planes perpendicular to the DM axis are FM arranged. The competition between the DM and FM interaction gives rise to a helicoidal spin structure along the DM axis and the model becomes effectively one dimensional. The classical effective 1D hamiltonian in the continuum limit reads: where z labels the coordinate along the DM axisẑ, the prime stands for the derivative with respect to z,n is a unit vector in the direction of the spin S, so that S = Sn, with S being the spin modulus, and Λ is the system length along the DM axis, which is assumed to be large. The constants entering Eq. (1) are related to the interaction strengths as follows: let J and D stand for the strengths of the FM and DM interactions, respectively, K for the strength of the single-ion anisotropy interaction, and a for the underlying lattice parameter; the energy H is measured in units of JS 2 /a, and therefore q 0 = D/Ja, γ = K/Ja 2 , and β = (gµ B /JSa 2 ) H, where H is the applied magnetic field. The parallel magnetic field is thus proportional to β z and, without any loss, we can take the perpendicular field proportional to β x and set β y = 0. The parameter q 0 is the helix pitch at zero magnetic field, so that its period is L 0 = 2π/q 0 . We set q 0 = 1 in the computations, what merely amounts to set the unit length equal to L 0 /2π. At zero temperature the classical ground state is given by the minimum of the hamiltonian (1), which is a solution of the corresponding Euler-Lagrange (EL) equations. The infinitely many solutions are characterized by the boundary conditions (BC). On physical grounds, we expect a periodic ground state for an infinite system. Then, if L is the period, the appropriate BC isn(0) =n(L). This however does not guarantee periodicity because, as the EL equations are second order, periodicity requires also the equality of the first derivatives, a condition which cannot be generally imposed since the associated boundary value problem (BVP) would be overdetermined. Hence, the value of the spin at the boundaries,n(0), has to be tuned in order to obtain periodicity with period L. This can in principle be done as we have the two degrees of freedom associated ton(0) and only two of the three Eqs. (2) are independent, due to the constraintn(z) ·n ′ (z) = 0. The periodic solutions are thus solely characterized by the period L. As the solution of a BVP is not necessarily unique, there might be several periodic solutions associated to a given L. We do not expect this on physical grounds, and indeed we never found more than one solution for the BVP in the course of our numerical computations. The equilibrium period is obtained by minimization of the energy density, E = H/Λ, which for a periodic state is given by where h(z) is the integrand entering Eq. (1). The approach proposed here is a generalisation of the well known procedure of expanding the magnetisation in harmonic modes and selecting the mode which minimises the free energy. The present problem, however, is highly non-linear and all modes contribute to the solution, which has to be found by numerical techniques. Dzyaloshinskii 5,6 solved analytically the case of a purely perpendicular field using a similar approach: he got the general solution of the differential equations and minimized the free energy in terms of the non-trivial integration constant, the jacobian elliptic modulus, which is directly related to the period of the IC structure.
The form of the EL equations depends on the coordinate system chosen to parametrizen. As the whole sphere cannot be smoothly parametrized with a single coordinate set, we used two different coordinate sets: a) coordinate set I, (ξ, ϕ), with ξ ∈ (−∞, ∞), and ϕ ∈ [0, 2π], so thatn = (cos ϕ, sin ϕ, ξ)/ρ, where ρ = 1 + ξ 2 ; and b) coordinate set II, (ϑ, φ), with ϑ ∈ [0, π] and φ ∈ [0, 2π], so thatn = (sin ϑ sin φ, cos ϑ, sin ϑ cos φ). Coordinate set I is closely related to the polar coordinate set with polar axis along the DM axis (ẑ) and the polar angle θ related to ξ by cot θ = ξ. This coordinate set is singular on ±ẑ, where ξ diverges. Coordinate set II is a polar coordinate system with polar axis alongŷ, and therefore is singular at ±ŷ, where sin ϑ = 0. The EL equations and the BC read in coordinate set I and in coordinate set II Asn contains two degrees of freedom and the differential equations are second order, the general solution contains four independent constants. The condition n(0) =n(L) removes two of them. We used translational invariance to eliminate one more, so that it remains either ξ 0 or φ 0 as tuneable parameter. The reason is that the configuration that minimizes the energy will have at least one point with the spin lying on the plane determined by the DM axis (ẑ) and the magnetic field (x). Using translational symmetry we can choose this point as z = 0, so that in coordinate set I we have ϕ(0) = 0 and in coordinate set II ϑ(0) = π/2. The ground state (or equilibrium) configurationn(z) with z ∈ [0, L] can be visualized as a closed path on the unit sphere, with the coordinate z acting as the parameter of the curve. Fig. 1 displays such paths, given by the red trajectories, in different situations. The green axis represents the DM axis (ẑ) and the FFM state at the transition point is represented by the blue axis. Coordinate set I and its associated BC is appropriate when the path does encircle the DM axis; coordinate set II is suitable when the paths are not too close toŷ, which lies on the equatorial plane. As already mentioned, notice that the BC (4c) forces the path to encircle the DM axis. When the values of the magnetic field are such that the equilibrium state is represented by one path which does not encircle the DM axis, as in the top panels of Fig. 1, the BVP (4) has no solution, and we are forced to use coordinate set II and solve (5). On the other hand, on a wide part of the phase diagram both coordinate sets can be used. In such cases, we solved both BVPs and got the same solutionn(z) within the tiny numerical uncertainties (see the appendix, where details on the numerical procedures are given).

III. PHASE DIAGRAM
To obtain a phase diagram we compare the energies of the IC state, E IC , which is computed numerically following the lines outlined in the previous section, and of the FFM state, E FFM , for which ϕ = 0 and the constant value of ξ is given by the solution of For small β, E FFM > E IC and the IC state is the ground state. The transition to the FFM state takes place when E FFM = E IC , and is continuous if at this point the IC state merges smoothly with the FFM. Otherwise the transition is discontinuous. The IC state can be visualized as a closed curve on the unit sphere while the FFM state is represented by a single point (Fig. 1). As discussed in the introduction, there are two mechanisms by which the IC state can be continuously transformed into the FFM state. The first possibility is displayed in the top left panel of Fig. 1: as the magnetic field is tuned to its critical value the IC curve reduces its size until it collapses onto the FFM state. In this case a helical conical state which, in the vicinity of the transition point, revolves around the direction of the FFM and not around the DM axis, becomes completely parallel at the transition point. The bttom right panel of Fig. 1 illustrates the second possibility. The length of the IC curve on the unit sphere remains finite as the transition point is approached. Near the transition, the vast majority of the spins, however, are concentrated on a narrow arc close to the FFM state and the number of spins lying on the remaining part of the curve becomes negligible, so that a CSL is formed. At the transition point the period of the soliton diverges and therefore its fundamental wavevector, q = 2π/L, tends to zero. Fig. 2 displays the phase diagram in the (β x , β z ) plane for γ = 2.584q 2 0 (a) and γ = 0 (b). We shall disccuss the γ = 0 case in section V and here the γ = 2.584q 2 0 case, which is relevant for CrNb 3 S 6 (see section VII). We see that the transition is continuous in the vicinities of the β x = 0 and the β z = 0 axes, while it is discontinuous in the intermediate regime. The behaviour of E FFM − E IC as a function of β z for fixed β x is illustrated in the top panels dz| is the magnetization in suitable units, as an "order parameter" (OP) which vanishes in the FFM phase. Its behaviour in each case is displayed in the bottom panels of Fig. 3.
For conciseness, we shall call the continuous transition lines which touch the β x = 0 and the β z = 0 axes respectively the continuos transition line 1 (CTL1) and 2 (CTL2), and we shall use DTL for discontinuous transition line. On the CTL1 the FFM state is reached continuously by the "closing cone" mechanism, as illustrated in Fig. 1 (top left). The period L remains finite and close to L 0 at the transition point. On the other hand, on the CTL2, the period of the IC diverges and the FFM state is continuously reached by the same mechanism as in the β z = 0 case. The bottom right panel of Fig. 1 is an instance of this mechanism.
On the DTL the IC and FFM state become degenerate and neither the "cone closes" nor the period diverges. These situations are illustrated in the top right and bottom left panels of Fig. 1. The OP is discontinuos on this line. Fig. 4 (left) shows the OP jump along the transition line, parametrised by β x . Along the DTL the IC state coexists with the FFM state as they are degenerate in energy. On each side of the line either the IC or the FFM state are metastable. This is illustrated in Fig. 4 (right), where the energy difference is displayed as a function of q = 2π/L in a case in which the IC state with q ≈ 0.84 q 0 is metastable and the FFM (q = 0) is stable. Hence, phenomena like phase coexistence, domain formation, and hysteresis are expected in some regions of the magnetic phase diagram. The metastability does not cause numerical problems since L is fixed in the numerical computations, and the BC (4c) prevent the switching between the IC and the FFM states. The DTL is separated from the CTL1 and the CTL2 by two tricritical points, called respectively TC1 and TC2. In Fig. 4 (left) these are the points where the OP jump ceases to vanish, and are signaled by the two green circles. For γ = 2.584q 2 0 TC1 is located at (β t1 x , β t1 z ) = (0.1593, 5.9552) and TC2 at (β t2 x , β t2 z ) = (0.498, 3.5301). Fig. 5 displays the fundamental wavevector q as a function of β z for different values of β x , including β t1 x = 0.159 and β t2 x = 0.498 (dashed lines with black symbols). Notice that q is almost constant (q ≈ q 0 ) for β x < β t1 x , while it decreases with β x and with β z for β x > β t1 x . The variation of q becomes very abrupt close to the tricritical point TC2.
The spatial dependence of the ground state as the transition is approached is shown in Fig. 6, where the cartesian components ofn are displayed as a function of z/L. Cartesian components are shown since they are regular for any β. They have been computed from the solutions obtained either with coordinate set I or II. The left panels correspond to β x = 0.15 and shown(z) for different values of β z which tend to a transition point on CTL1. We observe that the variation ofn(z) is distributed smoothly over the whole period for all β z , so that no soliton is formed. Also, the amplitude of the oscilations decreases smoothly to zero as the transition point is approached, attaining continuously the FFM. The right panels correspond β x = 0.55 and the values of β z tend to a transition point on CTL2. Observe that the variations ofn(z) are gradually concentrated on a narrow section of the period, and a soliton lattice is formed. It is a relative narrowness: the size of the region wheren(z) varies noticeably is approximately constant, but the period increases as the transition is approached. It is interesting that the soliton lattice can be created by increasing the component of the magnetic field parallel to the DM axis, keeping constant 0 1.0 the perpendicular field.

IV. SINGULAR BEHAVIOUR ON THE CONTINUOUS TRANSITION LINES
On the continuous transition lines CTL1 and CTL2 the IC and FFM states merge continuously but, as in any phase transition, physical observables may present peculiar singularities. Let us analyze the nature of these singularities. The OP vanishes linearly as the CTL1 (β x < β t1 x ) is attained from the IC phase (bottom left pannel of Fig. 3). Thus, no singularity appears on CTL1. As TC1 is approached, the slope of the OP becomes more abrupt and becomes singular at TC1. We assume that at TC1 the OP vanishes as a power law, OP ∼ A(β t1 z − β z ) n/m , with n and m small integers. The best fit, which is indeed very good, gives β t1 z ≈ 5.9552506 and n/m = 0.5002 ≈ 1/2.
As CTL2 is attained from the IC phase, the OP vanishes more abruptly than a power (bottom right pannel of Fig. 3). The OP scales with the period as 1/L, so that the OP singularity is directly related to the L singularity. For β z = 0 the following scaling law is derived from the analytically known solution: where β xc = (π 2 /16)q 2 0 is the critical field at β z = 0. Along CTL2, (β t2 x < β x ≤ β xc ), a generalization of the above scaling law, given by holds. In this expression β stands either for β x or for β z , and the other component of β has to be kept constant.
The parameters A and B vary smoothly along CTL2. The value of B depends on whether CTL2 is attained keeping β x or β z constant, but A, which characterises The scaling of L given by (8) is a general feature of the formation of the CSL, and it can be used to fit the experimental results near the transition point for oblique fields. The behaviour of the parameter A then can be used to locate the tricritical point TC2.

V. EFFECT OF THE EASY-PLANE SINGLE-ION ANIOSTROPY
The easy-plane anisotropy has an important effect on the phase diagram, as can be seen in Fig. 2b, where the γ = 0 case is displayed. We see that the CTL2 disappears, the DTL reaches the β z = 0 axis, and only one tricritical point (TC1) appears. For γ > 0 the phase diagram is qualitatively similar to the case discussed in the previous sections (γ = 2.584q 2 0 ). The transition lines in the (β x , β z ) plane for different values of γ are shown in Fig. 8. By increasing γ the lengths of CTL1 and CTL2 respectively decreases and increases. In a tridimensional parameter space (β x , β z , γ) there is a phase transition surface divided by two tricritical lines into three sectors, one in which the transitions are discontinuous and two in which they are continuous.
Not surprisingly, the easy-plane anisotropy provides stability to the CSL formed when the field is purely perpendicular, so that for large γ large parallel fields are necessary to modify the behaviour of the CSL and its transition to the FFM state. The effect of easy-plane anisotropy is important since low-anisotropy compounds described by the model studied in this work, as for instance thin films of MnSi 16 , will be characterized by low values of γ.
The case of easy-axis aniostropy (γ < 0), not analyzed in this work, is theoretically interesting -and perhaps also phenomenologically-as peculiar phenomena may take place in the crossover to Ising-like behaviour as |γ| increases 17 .

VI. COMPARISON WITH SIMPLIFYING APPROXIMATIONS
It is interesting to compare the results obtained here with those obtained by making approximations which drastically simplify the mathematical problem. We consider two of such approximations: the conical approximation of Ref. 15 and the decoupling approximation developed below.
In the conical approximation the spatial variation of the spin component parallel to the DM axis, n z , is neglected, and the ground state is obtained by minimizing the energy density over spin configurations with constant n z = cos θ 0 . The EL equations for the non-constant components give a sine-Gordon equation similar to that of the β z = 0 model with an effective perpendicular field β eff x = β x / sin θ 0 . Then, the ground state has the form 15 n x (z) = sin θ 0 cos ϕ(z), (9) n y (z) = sin θ 0 sin ϕ(z), (10) n z (z) = cos θ 0 , where ϕ(z) = −2 am β x /(κ 2 sin θ 0 )z, κ , (12) and am(x, k) is the Jacobi amplitude function with elliptic modulus κ. The energy density associated to the above solution depends on the two parameters cos θ 0 and κ, and reads 15 where K(κ) and E(κ) are the complete elliptic integrals of the first and second kind, respectively. The first three terms give the energy of the conical helix with wavevector q 0 , while the term proportional to β x is the difference between the conical helicoid (12) and the conical helix. We write the energy in this way for later convenience. Minimization of (13) with respect to θ 0 and κ gives the phase diagram. The IC-FFM transition is everywhere discontinuous except at the two end points of the transition line. Notice that the conical approximation is exact in the two limiting cases β x = 0 and β z = 0. The goodness of the approximation depends on the magnitude of the easy-plane anisotropy. Fig. 9 displays the phase transition lines obtained in this work by numerical integration of the EL equations (EL solution), and in the conical approximation, for γ = 0 (top) and γ = 2.584 (bottom). In the γ = 0 case the transtion line in the conical approximation deviates considerably from the line given by the EL solution. However, for γ = 2.584 the conical approximation predicts a transition line which is very close to the EL solution line. As the energy of the approximate conical helicoid is always higher than the energy of the ground state the approximate transition line lies always below the EL solution transition line, as can be appreciated in Fig. 9.
It is easy to understand why the conical approximation locates the transition line very accurately as the magnitude of the easy-plane anisotropy increases. The energy H given by Eq. (1) contains a term γn 2 z − β z n z = γ(n z − β z /2γ) 2 − β 2 z /4γ. If γ is very large the energy penalty caused by a large deviation of n z from the constant value β z /2γ cannot be compensated by a gain provided by the other terms in H. In the vicinity of the transition, however, the energy of the IC state is only slightly lower than the FFM energy and the conical approximation does not capture the nature of the transition. The differences between the IC and the FFM states are precisely due to the spatial variation of the spin. The small variation of n z allows the IC to reach the FFM state continuously. Consider for instance in the transition at constant β x depicted in the top left panel of Fig. 1. The transition takes place continuously when the spin path (red points) collapse onto the FMM (blue axis). This is not possible in the conical approximation, in which the spin configuration paths are forced to follow the parallels of the unit spheres of Fig. 1. As it gives discontinuous transitions along the whole transition line, the conical approximation does not predict the continuous transitions nor the tricritical behaviour.
Let us develop a further approximation, which we call the decoupling approximation, in which we assume also a conical ground state. We shall denote by β xc = (π 2 /16)q 2 0 and β zc = 2γ + q 2 0 the parallel and perpendicular critical fields, respectively. We notice that, for large γ, the value of cos θ 0 is determined basically by the three first terms of (13). The term proportional to β x , which, we recall, is the difference between the energies of the conical helicoid (12) and the conical helix with ϕ(z) = q 0 z, gives a minor contribution. Indeed, at low β x and β z the conical helicoid is a slightly distorted conical helix; at low β x and large β z the term proportional to β x is clearly a small correction; and for β x close to β xc the energy of the conical helicoid is very close to the energy of the FFM state, which happens to be close to the energy of the conical helix 18 . Hence, in the decoupling approximation we have cos θ 0 = β z /β zc . Minimization of the energy with respect to κ gives an equation similar to that of the well known β z = 0 case: E(κ)/κ = sin θ 0 β xc /β x . The phase transition takes place continuously when κ = 1, that is when β x /β xc = sin θ 0 , and, taking into account the value of sin θ 0 , we get the following expression for the transition line: The continuous green lines of Fig. 9 are the curves corresponding to the above equation. We see that they describe very accurately the transition line for large easyplane anisotropy, but it departs notably from the EL solution for low γ. Interestingly, equation (14) shows that the form of the transition line for high anisotropy is universal when expressed in terms of the dimensionless fields β x /β xc and β z /β zc . We see that the decoupling approximation is as accurate as the conical approximation, or even more, and it has the virtue of providing a simple analytic expression for the transition line. On the other hand, the decoupling approximation suffers from the same limitations as the conical approximation and does not predict correctly the nature of the transition, as it gives always continuous transitions.

VII. APPLICATION TO CrNb3S6
The general results obtained in this work can be readily applied to CrNb 3 S 6 . For this compound, a ≈ 1.2 nm. Other measurable parameters like the critical parallel and perpendicular fields, H xc and H zc respectively, and the the zero field helix period, L 0 , depend strongly on the sample. As an instance 8 , we take L 0 ≈ 48 nm and H xc ≈ 2300 Oe. A value γ = 2.584q 2 0 ensures the relation H zc ≈ 10H xc tipically observed experimentally. The large easy-plane anisotropy implies that the location transition line is accurately given by the results of the conical approximation of Ref. 15, and, morover, it is given by the Eq. (14) derived in the present work, so that The critical fields H xc and H zc depend strongly on the sample. Due to impurities, crystalline defects, etc., each sample will be theoretically described by a different set of the model effective parameters. The anisotropy, however, is expected to be large for all samples, so that Eq. (H t2 x , H t2 z ) ≈ (1860, 12490) Oe. A CSL can be formed by approaching the CTL2, increasing the perpendicular field H x with the parallel field H z hold constant below 12500 Oe, or increasing the parallel field with the perpendicular field hold constant above 1850 Oe, or increasing a field tilted from the DM by an angle α smaller than 81 degrees.
Furthermore, a line of discontinuous transitions appears in the phase diagram of CrNb 3 S 6 , and therefore phenomena typical of first order transitions, like phase coexistence and hysteresis in magnetisation and magnetoresistance, are expected.

VIII. CONCLUSIONS
In this work the important and long-standing problem of the determination of the phase diagram of the mono-axial chiral helimagnet at zero temperature in the magnetic filed plane has been addressed. The transition line from the incommensurate phase to the commensurate forced FM state has been thoroghly analysed as a function of the single-ion aniotropy and the nature transitions has been elucidated. We have shown that generically the transition line is composed by a line of discontinuous transition separated from two lines of continuous transitions by two tricritical points. The nature of the singularities at the transition line has also been thoroughly analysed and we showed that the transition from the chiral soliton lattice to the forced FM state along the continuous transtion line is characterized by logarithmic singularities analogous to those characteristic of the case of purely perpendicular magnetic field. This characteristic singular behaviour may be used to locate one of the tricritical points.
For large anisotropy the conical and the decoupling approximations give the location of the transition line very accurately, but they fail in characterising the nature of the transition, since the former predicts a discontinuous transition along the whole line and the later continuous transitions everywhere. Consequently, the tricritical behaviour is not predicted by none of these approximations. For large anisotropy we found that the form of the transition line is universal, given by Eqs. (14) or (15). The position of the tricritical points on the transtion line is not universal, however, and has to be computed for each sample.
Hence, unexpected predictions for the low temperature regime of CrNb 3 S 6 are given: discontinuous transitions, tricritical behaviour, and the universality (that is, the independence of the sample) of the form of the transition line due to the large anisotropy characteristic of this compound. These phenomena may have interesting applications in spintronics. It remains to study the very interesting question of how the phase diagram evolve by increasing the temperature. Work in this direction is in progress 19 . of the energy as a function of the BC for fixed L is the periodic solution which satisfies Eqs. (2).
To solve numerically the BVPs defined by eqs. (4) or (5) we used a relaxation method as described in Ref. 20 (see also Ref. 21). The BVP solver works as follows. First, the two second order differential equations are converted in a set of four first order differential equations in the usual way, by introducing two new variables ω and v and two new equations which relate them to the derivatives of ϕ and ξ: ω = dϕ/dz, v = dξ/dz. Then, this system of first order differential equations is converted in a system of finite difference equations by substituting the derivatives by forward finite differences. Therefore, a mesh with N + 1 points (z 0 = 0, z 1 , . . . , z N = L) is introduced. For simplicity, we used a regular mesh, so that z n − z n−1 = L/N for n = 1, . . . , N . The finite difference equations and the BC form a set of 4(N + 1) non-linear algebraic equations (4N of them given by the difference equations and the remaining 4 provided by the BC) with 4(N + 1) unknowns, which are the values of ξ(z), ϕ(z), v(z), and ω(z) at the mesh points. The algebraic equations are solved numerically with a Newton algorithm, which proceedes iteratively starting from one initial guess which has to be supplied. Subsequently, E is computed with a simple quadrature algorithm using the values h(z) calculated at the mesh points.
Periodicity is enforced by tuning ξ 0 until the conditions ω(L) = ω(0) and v(L) = v(0) are met. The two equations can be enforced by tuning only one variable, ξ 0 . The reason is that we exploited translational symmetry to eliminate the freedom in the BC for ϕ, which in general should read ϕ(0) = ϕ 0 and ϕ(L) = ϕ 0 + 2π. Hence, periodicity can be enforced by tuning ϕ 0 and ξ 0 . Translational invariance implies that we can set ϕ 0 = 0, taking into account that, on physical grounds, there has to be a point in which the spin lies in the plane formed by the DM axis and the magnetic field (i.e., with ϕ = 0). Hence, we enforce the condition v(L)−v(0) = 0 by tuning ξ 0 via a simple bracketing-bisection algorithm, and the other condition is automatically satisfied. Finally, the minimum of E as a function of L is obtained also with a simple bracketing-bisection algorithm.
To control the numerical errors we set very demanding values for the convergence tolerances of the different algorithms: 10 −14 for the BPS solver, 10 −12 for the computation of E, and 10 −11 for the bracketing-bisection algorithms. The effect of the discretization was taken into account when computing E via the quadrature algorithm. The mesh was refined iteratively, the number of points being doubled at each iteration, until the difference between the values of E computed in two succesive iterations was smaller than the tolerance (10 −12 ). In this way we were able to get the ground state with high accuracy.
The main issue is the convergence of the BVP solver, which works iteratively so that to get convergence the initial guess functions have to be close enough to the actual solution. We achieve this by varying β z by small steps and using the solution found at one β z as the seed for the next step. We start at β z = 0, taking advantage of the fact that the solution is analytically known there. This procedure is sound as we expect the IC structure to be continuous on β, as it is in the cases where the field is either perpendicular or parallel. The results confirm this point.