Skyrmion Tubes in achiral nematic liquid crystals

We analyze the interaction with uniform external fields of nematic liquid crystals within a recent generalized free-energy posited by Virga and falling in the class of quartic functionals in the spatial gradients of the nematic director. We review some known interesting solutions, i. e., uniform heliconical structures, which correspond to the so-called twist-bend nematic phase and we also study the transition between this phase and the standard uniform nematic one. Moreover, we find liquid crystal configurations, which closely resemble some novel, experimentally detected, structures called Skyrmion Tubes. Skyrmion Tubes are characterized by a localized cylindrically-symmetric pattern surrounded by either twist-bend or uniform nematic phase. We study the equilibrium differential equations and find numerical solutions and analytical approximations.


I. INTRODUCTION
The achiral nematic N phase is surely the most common state for thermotropic liquid crystals. Due to the uniaxial symmetry of the constituent molecules, it is possible to describe this liquid crystalline phase by means of a director field n ∈ RP 2 , prescribing point by point the average orientation of the molecular axes. In particular, in the N phase the ground state alignment of these axes is parallel to a fixed direction n ≡ n 0 . On the other hand, in the chiral nematic N * phase, formed by enantiomorphic molecules, the minimum energy configuration n = n c (r) is spontaneously twisted in a right-angle helix, whose pitch P usually lies in the range of micrometers.
However, new classes of nematics are expected when molecules are less symmetric, as for example, the biaxial nematic phase [1][2][3]. Moreover, for strongly bent mesogenic molecules, a new modulated nematic phase, now recognized as the twist-bend nematic N TB phase, has been recently observed and reported in several works, starting from the breakthroughs in [4][5][6]. It turned out that this phase is stabilized below the usual N phase and, although formed by achiral molecules, it exhibits doubly degenerate chirality, consisting of right and left Meyer's heliconical domains [7]. Thus, the appearance of the N TB phase represents a particularly interesting case of spontaneous breaking of the chiral symmetry. Since in these structures the director n is tilted by a fixed angle 0 < θ 0 < π/2, they may look similar to the smectic SmC * phases. However, at variance with them, the heliconical textures do not possess any layer periodicity. Moreover, the helical pitch is much smaller than the cholesteric one, i. e., of the order of 10 nm [4,5].
Several papers addressed the theoretical analysis of the N TB phase, both from the phenomenological and the static continuum theory points of view [8][9][10][11]. More specifically, in [9] a N -N TB phase transition was described by means of a generalized Maier-Saupe molecular field theory. In [10], a generalized Landau -de Gennes theory was applied to investigate the modulated nematic phases, possibly generated by achiral and intrinsically chiral bent mesogenic molecules. In [11], the N TB phase was studied as a mixture of two different ordinary N phases, both presenting heliconical structures with opposite helicities. A quadratic elastic theory, still featuring four Frank elastic moduli, was used for both helical phases. Similar models were proposed in [12][13][14], where also the effects of an external magnetic or electric bulk field were investigated. Moreover, authors in [15,16] presented coarse-grained elastic models which, similarly to the model for SmA * [17], make use of an extra scalar order parameter.
The prediction of the N TB phase dates back to the seminal paper [18] by Dozov, in which an elastic instability model was proposed with a bend constant K 33 turning negative. Higher derivative terms were added to the standard Frank-Oseen elastic energy in order to bound the energy from below. According to Dozov's model, depending on the ratio of K 11 and K 22 in the high-temperature non-modulated nematic regime, the low-temperature nematic phase can show either the twist-bend modulation when K 11 > 2K 22 or the splay-bend modulation otherwise. The same model predicted the existence of a nematic phase with spontaneous bend distortions [7,18]. However, unlike a twist, a pure bend distortion cannot fill the space without introducing frustration, possibly relieved by defects. This is certainly not the case for the mixed twist-bend distortion. Indeed, in [19] it was shown that in three space dimensions there exist only two families of director configurations which have uniform non-zero distortion characteristics at any point in space. It turned out that these latter configurations correspond to the right and left Meyer's heliconical domains, which form the N TB phase. Any other director field, apart from the constant nematic director, would be geometrically frustrated and become nonuniform if requested to fill the whole space.
The natural successive step was to see whether it is possible to build an elastic free-energy that penalizes the departures from one of these uniform director fields. Since in the uniform heliconical phases only one of the distortion characteristics vanishes, namely the splay one, Frank's quadratic theory is no longer sufficient. Thus, a higher-order elastic theory, in which the bend elastic constant may turn negative, was proposed in [19], allowing for fourthorder powers of ∇n in the free-energy. The author focused on an achiral scenario where the N TB phase has been experimentally identified and deliberately built his generalized elastic free-energy with the symmetry of the intended heliconical ground state, i. e., its double degeneration for right and left helicities. This choice makes the free-energy depend on only six elastic constants: three for the quadratic part and three for the quartic one. Then, for suitable choices of the elastic constants, it was shown how either the standard nematic or the heliconical phases minimize the proposed higher order free-energy. More specifically, the theory predicts the N TB phase arising from the standard nematic one for sufficiently negative values of the bend constant, passing through an intermediate pure bend state.
In [20] we reviewed the theory presented in [19] and found that, in the same region where N TB is preferred, localized excitations of the heliconical ground state are possible. In particular, we showed how axisymmetric structures, with a radial dependence of the conical angle and an additional twist around the heliconical axis, are admissible states of the generalized elastic theory, with energies falling in between those of the heliconical ground state and the nematic alignment. We found that our soliton configurations resemble interesting axisymmetric structures recently observed in chiral nematics and chiral ferromagnets [21][22][23][24], namely Skyrmion tubes. In contrast with these latter configurations, ours can be generated in an achiral framework without the need of external frustration.
Emergent topological defects in condensed matter systems are drawing much attention, particularly due to its potential technological applications, and Skyrmion tubes are not an exception. For instance, they have been recently proposed as magnonic waveguides channeling spin waves, based on the propagation of their breathing and rotational modes [25]. However, theoretical studies such as [23,24] have been focused on its realization in ferromagnets, although its experimental attainment seems to be easier in liquid crystals [23]. In this context, [20] paved the way for a better theoretical understanding of Skyrmion tubes in liquid crystals, where other localized configurations such as helicoids or Skyrmions are well known [26][27][28][29][30][31][32][33][34].
In the present paper, we generalize our previous work [20] with the addition of an external uniform magnetic field. We find Skyrmion-tube-like nematic textures that form when a uniform magnetic field is applied along the axis of a heliconical state. These Skyrmion tubes are surrounded by either a nematic uniform phase or by a uniform twist-bend phase. The paper is organized as follows. In Sec. II, we briefly revise previous material and set the model. Then, we study the interaction of an external magnetic field and find a uniform heliconical state. We also investigate the transition to the standard uniform nematic phase as the magnitude of the external field is increased from zero to a critical value. In Sec. III, we study nonuniform localized states under the action of the external field by imposing boundary conditions at the center of the heliconical state and at infinity. We find Skyrmion-tube configurations where the nematic texture is nonuniform in a localized radial region immersed in a uniform, either nematic or heliconical, state. Finally, in Sec. IV we draw our conclusions and outline future investigations.

II. TWIST-BEND PHASE UNDER EXTERNAL FIELDS
Nematic liquid crystals are usually modeled by Frank's elastic free-energy density. This is a general positive-definite quadratic form in the spatial gradients ∇n of a unit vector, the nematic director n, and it is written as where K 11 , K 22 , K 33 , and K 24 are the Frank elastic constants and they are such that known as Ericksen's inequalities [35]. The term K 24 is a null Lagrangian, it can be integrated over the domain B occupied by the nematic medium without producing any contribution to the total free-energy, provided that n is assigned over the boundary ∂B.
In [36] it was shown that Frank's elastic free-energy density can be written as a quadratic form in four quantities (S, T, b, D) as follows where S = divn is a scalar called splay, T = n · curln is a pseudo-scalar named twist, and B 2 = b · b, with the vector b = n × curln being the so-called bend . D is a symmetric traceless tensor such that Dn = 0. Accordingly, it can be given the form where q is the positive eigenvalue of D, named by Selinger [36] as biaxial splay, and n 1 and n 2 are the eigenvectors orthogonal to n. The tensor D can also be given the following form in terms of ∇n The quantities (S, T, b, D) are independent from one another and are called measures of distortion. Frank's energy (3) admits as global minimizer the state which corresponds to any constant field n ≡ n 0 . In [19] it was put forward a new energy functional with quartic powers of measures of distortion (S, T, b, D) as follows This represents the lowest order free-energy density that, for a suitable choice of the elastic constants, admits as global minimizer the so-called heliconical uniform distortion state [19,20], as opposed to the uniform state (6). By directly comparing (7) with (3) we get the following formal identification but as shown below k 3 can also assume negative values. From (8), it is clear that the number of independent Frank elastic constants is reduced from four to three as K 22 − K 24 = K 24 . This assumption is due to the choice of the same elastic constant in front of T 2 and 2trD 2 in the quadratic part of (7), and T 4 and 4(trD 2 ) 2 in the quartic part. This latter condition is related to the heliconical global minimizer of (7) (see below), which is such that T 2 = 2q 2 [19,20]. Hence, the free-energy density must be invariant under the transformation T 2 ↔ 2q 2 , implying that only the combinations T 2 + 2q 2 and T 4 + 4q 4 appear in the free-energy density. The above energy density turns out to be coercive provided that which is the condition of positive definiteness of the quartic part of (7). In terms of n and its gradients ∇n, (7) can be written as follows [20] Correspondingly, the free-energy stored in a region B occupied by the liquid crystal is given by the volume integral As mentioned above, (10) admits, as global minimizer, the uniform heliconical state [19,20]. This latter can be written as follows n h = sin θ 0 cos βze x + sin θ 0 sin βze y + cos θ 0 e z , where e x , e y , e z are the Cartesian unit basis vectors in R 3 , θ 0 is the conical angle and β a parameter that provides the pitch P = 2π/|β| of the twist. Here, we assume that β is a characteristic parameter that depends on the elastic constans k i only and needs to be optimized. The three-dimensional representation of such configurations is displayed in Fig. 1, where a set of (x, y)-plane cross sections showing how the configuration changes along z and a specific helix line are depicted. The nematic director n h rotates around e z making a fixed cone angle θ 0 with the rotation axis e z , which is called the helix axis. The structure (12) describes therefore the heliconical distortion predicted by Meyer [7] and it corresponds to the twist-bend liquid crystal phase N TB , experimentally detected in 2011 [37]. It is worth noticing that formula (12) also describes the nematic phase N when θ 0 = 0 and the chiral nematics when θ 0 = π 2 , implying that the twist-bend phase represents a structural link between these two extreme phases. Of course, as also observed in [19], the heliconical configurations cannot be minimizers of the standard Frank elastic energy, and a new elastic theory as (7) was needed to accommodate the heliconical phase as a ground state.
The corresponding free-energy density reads which depends on the pitch-related parameter β and the conical angle θ 0 , and is minimized by the values where ± signs label a counterclockwise or a clockwise heliconical configuration, and θ 0 = arcsin 2k 2 k 5 + k 3 k 6 (2k 2 k 5 + k 3 k 6 ) + 2(k 3 k 4 + k 2 k 6 ) . In [19,20] it was shown that, in order to have the heliconical states (12), the following constraints on the elastic constants must hold When an external magnetic field H is applied, the free-energy density (10) turns into where we added the term χ a being the magnetic susceptibility of the liquid crystal material. Otherwise, the free-energy can be also written as follows Correspondingly, the stored free-energy in a region B occupied by the liquid crystal is given by the volume integral The Euler-Lagrange equation associated with the above functional is given by where λ is a Lagrange multiplier for the unit director constraint. To obtain a pure equation it suffices crossing by n both sides of (21).
In the following, we will consider an external magnetic field along the z-axis, i. e., H = He z , and assume a nematic director field as in (12) n h = sin θ 0 cos βze x + sin θ 0 sin βze y + cos θ 0 e z , with the helix axis parallel to the magnetic field. As mentioned above [23], β is to be taken as fixed by the elastic constants k i only in accordance with (14). Thus, we assume that the external field just affects the nematic director by a torque imparted to the liquid crystal molecules. Accordingly, the conical angle θ 0 will change. The interaction with an external field H = He z is explicitely given in terms of θ 0 by the term Correspondingly, the reduced free-energy density takes the form where t = sin 2 θ 0 . We need now to minimize (25) by solving the stationary condition which is equivalent to (21) under parametrization (22). The latter equation becomes and the real solution to (27) is given by where with In order to have a unique real solution, the discriminant ∆ must be positive. This is the case when the external field vanishes [19,20]. When H = 0, it is clear from the definition of the quantities d and η that ∆ increases with respect to the zero-field value, thus keeping the positive sign and still yielding a unique real solution. Therefore, which generalizes (15) when an external field is present. In addition, we also study the transition to the standard uniform nematic phase corresponding to t = 0 where the director field lines up with the direction of the external field, i. e., along e z . For this to occur, the external field should solve equation (27) when t = 0, i. e., which leads to the critical field In the following, when dealing with an external field, we will give it in terms of this critical field, H cr . Thus, when H ≥ H cr , t = 0, that is to say, the phase is standard uniform nematic. In terms of the elastic constants, the critical field becomes which has been obtained by choosing for the parameter β the expression in (14). In Figs. 2 and 3 we represent the conical angle θ 0 as a function of the external field for different choices of the elastic constants. It is interesting to see in Fig. 2 how θ 0 goes to zero when increasing H. In particular, we have plotted (blue line) the curves θ 0 (H) for the values of the elastic constants k 1 = k 2 = k 4 = k 5 = k 6 = 1 and k 3 = −3 (from now on we will call this choice of values the standard set) together with a decreasing of k 3 to -5 (red dashed line), which corresponds to a critical field H 2 cr = 605/3χ a , and k 3 = −10 (green dotted line) with H 2 cr = 845/χ a . Interestingly, we see that in all cases the approach to 0, when we are close to the critical field, occurs in the same way. However, the behaviour for small values of the external field is quite different. When decreasing k 3 , θ 0 takes longer to become significantly smaller and a more abrupt reduction appears. Similarly, we can consider the case of an increasing elastic constant k 4 . In Fig. 3, besides the diminution of θ 0 for the standard set of parameters, the cases of only changing k 4 from 1 to 5 and 10 (red dashed line and green dotted line, respectively) are shown. In this case, the critical field when k 4 = 5 is H 2 cr = 841/3χ a and H 2 cr = 10443/19χ a for k 4 = 10. Unlike the previous case, for an increasing k 4 the conical angle shrinks in a more regular way. This might be also favoured because, even in the absence of an external field, θ 0 considerably decreases with k 4 . For small values of the ratio h = H Hcr , t gets the form where all quantities subindexed with 0 refer to the corresponding quantities when H = 0. In terms of the asymptotic angle one has In this formula the dependency on k 3 is quite involved, but algebraic, then analytic. Notice that in this model the response of the heliconical angle θ 0 to the external field is just quadratic as in a sort of Kerr effect [38]. Despite the complexity of the system under study, even when only considering the simple twist-bend configuration, simulations in 3 dimensions of this ground state have been successfully undertaken. One should note that within this theoretical setup, performing 3D simulations to minimize (20) via (21) is a challenging problem. On the one hand, the boundary conditions for the vector director are not constant at infinity, together with a free energy which does not vanish asymptotically. In addition, a high accuracy in the numerical calculation is needed in order to exactly match the values obtained from the analytical study, not only for the conical angle θ 0 but also for the pitch β. In the case of the twist-bend ground state, this may be done in a lattice of a reasonable size. In this way, it has been confirmed, for the standard set of parameters (with a lattice spacing of 0.02 and a gradient flow method), the dependence of θ 0 on the external field as in Fig. 2 (solid line), besides validating the assumption of constant pitch made before. Furthermore, in Fig. 4 one can see the energy per volume, f , as a function of H/H cr , where for H = H cr we arrive at the nematic phase.
As mentioned in the introduction, several studies found that the pitch P of the modulated nematic structure falls in the scale of a few nanometers [4,5]. The same studies report, under Freeze-Fracture Transmission Electron Microscopy (FFTEM), the presence of stripe-textured fracture planes which indicates fluid layers periodically arrayed in the bulk with a spacing of P . On the other hand, the authors in [4] found that this periodic structure is achieved with no detectable associated modulation of the electron density, and so it is not accompanied by a mass density wave, revealing a nematic rather than smectic molecular ordering. Thus, the "layers" in 3D found in the FFTEM are not images of molecular scale interfaces, but rather are 2D surfaces of constant azimuthal phase of the heliconical precession, sometimes called "pseudolayers" [39]. Due to the undulation of these surfaces, the direction of the heliconical axes periodically changes accordingly. In particular, in [39] a buckling of these "pseudolayers" under the action of an external magnetic field was observed until they flatten for sufficiently high fields. There, a Helfrich-Hurault model was proposed to theoretically describe this phenomenon and infere the value of the associated elastic constants from the experimental values of the critical magnetic fields. On the other hand, the N T B configurations presented here, and derived from the phenomenological elastic theory (17), uniformly fill the space with heliconical axes parallel to each other, meaning that the surfaces of constant azimuthal phase are flat independently of the presence of an external field. In order to allow for displacements from flat surfaces, one should add extra energy terms to (17) as reported in [40] (see eq.(61) therein) by considering a suitable parametrization of n in terms of the unit normal to the surface , based on the chiral Frank-Oseen framework. It turned out that both the pitch and the tilt angle increase as the field decreases from the critical value for the cholesteric-isotropic transition, originating a new configuration called oblique helicoid. Moreover, the CB7CB is known to be one of the recently discovered materials exhibiting the N T B phase. However, it is difficult to compare the results presented here with those obtained in [41] since quite different situations are examined. Indeed, the oblique N * phase, obtained through the addition of the chiral dopant, shows a period of some microns, while, if no chiral molecules are added from outside, the periodic structures typical of the N T B phase have periods of a few nanometers. Furthermore, because of the dopant additive, the oblique helicoidal phase has the same chirality everywhere, whilst the N T B phase is formed by both left-and right-handed domains. Finally, by fitting theoretical predictions to experimental results, the authors in [41] found the bend constant to be significantly smaller than others but still positive. We must stress that here we investigate an achiral scenario in which this latter constant turns negative, possibly allowing for the stabilization of N T B and other non-uniform configurations. Indeed, in [41] it is stated that the oblique helicoids appear in a specific range of intensity of the external field, which allows for the chiral twist to compete with the torque of the field. Above this range, the homeotropic alignment is favored; below this range, the right-angle helix typical of cholesterics is stabilized. The limiting values for the amplitudes of the external field are both determined by the ratio 0 < κ = K 33 /K 22 < 1: a negative value of K 33 would imply imaginary values for all the observables taken into consideration in the theoretical model proposed in [41]. Thus, it is not surprising that in our case the pitch keeps independent of the external field.

A. Skyrmion tube parameterization
At variance with the previous section, here we consider the case of nonuniform distortions leading to localized states. Bearing in mind that the uniform distortions are heliconical states, we slightly depart from this case by considering still heliconical structures, but with a nonuniform conical angle and an additional precession around the heliconical axis. These structures give rise to localized cylindrically-symmetric configurations, which can be referred to as Skyrmion tubes (SkT), of the general form n(r, z, ϕ; β) = sin(f (r)) cos(ϕ + βz)e x + sin(f (r)) sin(ϕ + βz)e y + cos(f (r))e z .
The ϕ−dependence prescribes a winding performed by the director around the heliconical axis e z for fixed z, f (r) is the profile function describing the conical angle and β has the same meaning as in the previous section.
One might also think about an ansatz without the angular dependence given by ϕ. However, as we already discussed in [20] and also checked from numerical calculations and within this setup under an external field, stable solutions of this kind do not exist. Hence, this ansatz with no ϕ−dependence gives us the uniform distortion as the ground state where now, the conical angle θ 0 will also depend on the value of the external field as showed in the previous section. The winding around the z axis given by the azimuthal variable is of key importance to give rise to these localized Skyrmion tubes and prevent them from directly decaying into the ground state.
In order to justify the above ansatz, one way to proceed is to resort to the so-called reduction by variational point symmetries [43], where these latter transform both independent and dependent variables, leaving unchanged the value of the functional (20) and the associated full Euler-Lagrange equations (21). A symmetry reduction procedure leads to an exact form of the solution with a less number of independent variables, as in (37), and it produces the corresponding equations in the remaining unknown functions. These latter obey to a restricted class of boundary conditions. Finding all symmetries admitted by the equations (21) may be challenging because of their high complexity. Still, one can exploit the constructive assumptions of translational and rotational invariance [19]. To this purpose, it is useful to parametrize the director field n in terms of two real stereographic variables ρ and Φ according to the following correspondence By expressing the infinitesimal point symmetries in terms of vector-fields v in the space of the independent and dependent variables, one can prove that any 1-dimensional subalgebra of the class i) leaves invariant the external magnetic field, ii) is a variational symmetry and iii) admits the following invariants: Thus, one can claim that the original variational problem has symmetry invariant solutions of the form where ρ and F are functions to be determined by a pair of symmetry-reduced partial differential equations in the independent variables r, ζ only. Now, supposing that the function F is smooth and non trivial in the angle ϕ, then it has to be independent of r in order to avoid discontinuities. Similarly, a non trivial dependence on r of ρ implies independence of ζ, otherwise it may lead to singularity and multi-valuedness in ϕ. Furthermore, by choosing in particular F = −ζ and identifying 2/α = β, we get the ansatz (37). Once the ansatz is justified, one can place it directly in the functional and find the corresponding reduced Euler-Lagrange equation for f (r), as detailed in the following.
In order to have localized configurations, we may impose the boundary conditions f (0) = 0 and f (r → ∞) = θ 0 , θ 0 being a suitable conical angle to be determined. We also consider the case f (0) = π and f (r → ∞) = θ 0 . Then, to study these configurations, we need to reduce the general free-energy in order to translate the ansatz into the equilibrium equations. The reduced free-energy integrated over the unit cell 0, 2π β × [0, 2π] and over r ∈ [0, ∞] will take the form We are interested in the reduced free-energy per unit cell 0, This latter can be rewritten as follows where G i are reported in the Appendix A and we have now defined Hence, we arrive at the following Euler-Lagrange associated equation As stated above, we will look for localized solutions of the form (37). The radial profile function f (r) solves the ODE (47). Here, we want to study the asymptotic behaviour as r → ∞. The asymptotic state will then be denoted as n ∞ = sin θ 0 cos(ϕ + βz)e x + sin θ 0 sin(ϕ + βz)e y + cos θ 0 e z , where θ 0 is the asymptotic conical angle, i. e., f (r) → θ 0 as r → ∞.
In order to determine θ 0 , as in the case of uniform distortions, we follow the route of free-energy minimization. Alternatively, we could study the asymptotic behaviour directly from (47). To find this value, we just need to consider the stationary condition of the free-energy with respect to f . In this way, we get an asymptotic angle depending only on the elastic constants, the external field, and the β parameter given in terms of the elastic constants only (14). The free-energy to be minimized is the asymptotic expression ofF H as in (45). As r → ∞, we get where the function G ∞ 0 is obtained from the function G 0 by dropping all the terms 1/r and 1/r 3 and keeping only linear terms in r, i. e., where entailing that We then need to minimize the function where we are now using the asymptotic value of f , i. e., f → θ 0 as r → ∞. Dropping the r in the above expression, the function to be minimized is Upon setting as above t = sin 2 θ 0 , we arrive at The corresponding stationary condition reads that is which reproduces the same equation as (27). Accordingly, we obtain for the asymptotic conical angle the same expression as in the uniform heliconical configuration (31) together with its dependence on the uniform external magnetic field. In particular, it follows that the asymptotic angle θ 0 vanishes when H ≥ H cr (34) (see also Figs. 2 and 3).
In the next section, we will look for localized solutions to (47) by numerically minimizing the free-energy with a gradient flow method.

B. Numerical results
As mentioned above, we performed numerical simulations to minimize (20) via the Euler-Lagrange equation (21) in order to find the uniform heliconical state. The same approach might be applied to find general configurations for a given direction of the external field and the boundary conditions. However, due to the issues highlighted in the previous section, the study of localized solutions with 3D simulations is a challinging project, outside the scope of the present paper. Here, we find the configurations corresponding to the Skyrmion tubes by minimizing the free-energy (45) within the ansatz (37). Hence, the profile f (r) can be numerically obtained by using a gradient flow method. For this purpose, we consider a 1-dimensional lattice L of 1000 points with a lattice spacing ∆r = 0.02, with spatial derivatives approximated by a fourth-order finite difference. Regarding the boundary conditions, we will consider two different cases concerning the value at the origin: f (0) = 0 and f (0) = π, where f (r → ∞) = θ 0 . Although we know that solutions taking the value π at the origin have higher energy [20], they are interesting when placed under an external field since they may survive the application of a magnetic field bigger than H cr (34). This should not be surprising since θ 0 is a function not only of β and the elastic constants but also of the external field H. Indeed, as an increasing H will decrease the conical angle θ 0 , the configuration with f (0) = 0 will converge to the nematic phase when H = H cr , whereas the solution with f (0) = π will remain, interpolating between π at the origin and zero at infinity. In other words, we can say that when increasing the external field, the first class of solutions start to dilute in the ground state, i. e., the difference between the profile values at the origin and at infinite decreases until vanishing, leading to the uniform ground state n = e z everywhere.
On the other hand, as for the second class, when H ≥ H cr the director assumes the two degenerate ground state configurations n = −e z and n = e z at r = 0 and as r → ∞, respectively. Thus, a non-uniform configuration around the center survives and a cylindrical domain wall connecting the two different ground states arises.
In Fig. 5, solutions with elastic constants k 1 = k 2 = k 4 = k 5 = k 6 = 1 and k 3 = −3 (the standard set) for different values of the external field are shown for these two different boundary conditions at the origin. For β we have chosen its value in terms of the elastic constants as in (14), i. e., β = 5 in this case. All this gives us a critical field H 2 cr = 75/χ a . One can easily see the different effect of increasing the external field in each class of configurations due to the diminution of the conical angle θ 0 with it. Moreover, Fig. 6 depicts the energy per pitch, P = 2π |β| , of both classes of solutions, once the energy of the ground state, i. e., is subtracted, i. e., ∆E/P = (E SkT − E GS )/P (we have introduced the notation E for the energy coming from the numerical calculation to make clearer it is calculated in the finite lattice L). For the configurations with f (0) = 0, the excitation energy decreases when approaching the critical field, as expected since for H = H cr they converge to the nematic phase. Otherwise, if f (0) = π the excitation energy remains about the same value for small fields before rapidly increasing. One can also go further and explore a little bit the parameter space. For instance, we can decrease the elastic constant k 3 . Doing this, the behaviour of the profile f (r) is similar to the one shown in Fig. 5, since an increasing external field always implies a diminution of the conical angle θ 0 (see Figs. 2 and 3). Finally, one can consider how the excitation energy changes both with the external field and the elastic constants when f (0) = π. This has been depicted in Fig. 7 for different values of the coupling constants k 3 and k 4 in a logarithm scale. In this way, it is manifest how an increasing k 4 has an important effect, raising the excitation energies in a considerable way, much more remarkable than when the value of the elastic constant k 3 is varied with respect to the standard set of values. This may be due to the fact that, even in the absence of external field, the twist-bend phase is characterized by a smaller conical angle than in the standard case, making the central value f (0) = π present a higher deviation from the ground state.
In order to obtain a better visualization of the nematic texture, in Figs. 8 and 9 we reported a three-dimensional representation of a Skyrmion tube with a profile function taking the value π at the center, i. e., f (0) = π. As it is clear, a tube is defined as an axially-symmetric region where the conical angle changes from π to the asymptotic value θ 0 . When the external field increases, the tube is surrounded by a uniform nematic phase as θ 0 vanishes when H ≥ H cr . Actually, in Fig. 10 we also report the profile function when H = H cr , 1.5H cr , 2H cr , where the shrinking of the Skyrmion tube with an increasing field is manifest.
As for stability, we have some numerical evidence of it by energy considerations, although there is not a complete proof based on the second variation of the free-energy functional. Nevertheless, the calculation and analysis of the second variation are not trivial, even for the quadratic Frank's functional [44]. This problem would deserve a separate further treatment and consideration, which are certainly beyond the scope of the present paper.
As a conclusion, the numerical analysis shows the existence of two families of solutions to equation (47). In the Appendix B, following an approach outlined in [20], one may find a global analytical approximation which can fit these numerical solutions.

IV. CONCLUSIONS AND PERSPECTIVES
In this paper we studied the interaction of external uniform magnetic fields with achiral liquid crystals according to a generalized fourth-order elasticity theory recently put forward in [19]. This theory is encoded in the free-energy density (10) which is parameterized by six elastic constants: k 1 , k 2 , k 3 associated with the quadratic terms in the Frank free-energy and k 4 , k 5 , k 6 related to fourth-order contributions. Under appropriate constraints on the six elastic constants, the proposed free-energy admits heliconical configurations as global minimizers. They are characterized by a director forming a constant conical angle θ 0 with respect to a fixed axis, say z, as shown in Fig. 1, continuously precessing when moving parallel to this axis and turning completely round over the length of a pitch P = 2π/β (12). These heliconical configurations have been recently identified experimentally in the ground state of twist-bend nematic phase N TB .
When an external magnetic field is applied, an interaction term Γ H (17) is added to the free-energy density. We studied the effect of a uniform field along the symmetry axis of the uniform heliconical state. The heliconical uniform state preserves its pattern and only the conical angle is affected by the external field as a consequence of the magnetic torque imparted to the nematic director. As the magnitude of the magnetic field is increased, the nematic undergoes a transition from the N TB phase to the uniform nematic phase N where the director lines up with the external field. The transition takes place at a critical value H cr (34). The twist-bend phase has been further reproduced by 3D simulations through a minimization of the energy functional (20) via Euler-Lagrange equations (21). According to this, we should stress that the pitch remains constant, thus consistently confirming our initial assumption.
Following our previous work [20], now in the presence of the external field, we generalized the heliconical configurations to nonuniform localized axially symmetric structures with a variable conical angle (37). Actually in [20], in the zero-field case, we had shown that there exists an axially symmetric state where the conical angle depends on the radial distance from the symmetry axis, the z-axis in our parameterization, going from 0 (or π) to a θ 0 at infinity in the radial direction, while the director winds uniformly once around the z-axis. The conical angle profile goes from 0 to its asymptotic value θ 0 in an exponentially fast way, thus singling out a central core. These localized structures are usually referred to as Skyrmion tubes [21][22][23][24]. The free-energy corresponding to the configuration starting from 0 at the origin of the radial axis has lower energy with respect to the one starting from π, although they both are excited states with respect to the uniform heliconical distortion. As shown in the present paper, for sufficiently low applied external fields, the Skyrmion tubes still keep their basic structure as localized configurations with a central core surrounded by a uniform heliconical distortion and a winding around the symmetry axis. Once a critical threshold H cr (the same as in the uniform heliconical phase) is reached, the Skyrmion tubes undergo a change in their patterns according to the value of the conical profile function at the origin of the symmetry axis. More precisely, when the conical angle is zero at the origin and H ≥ H cr , the conical profile function vanishes and the liquid crystal undergoes the transition to the uniform standard nematic phase, where the director lines up everywhere with the external field. Thus, in this case the central core of the Skyrmion tube tends to disappear. On the other hand, when the conical angle takes the value π at the origin, at sufficiently high external fields the central core still survives and it gets surrounded by a standard uniform nematic pattern where the nematic director lines up with the external field. Thus, in the central core the conical angle changes rapidly from the value at the center to zero.
As stated above, we would like to stress here that the configurations found in the present paper are of the same type as the so-called Skyrmion tubes found in [23,24] and there described numerically in ferromagnets and experimentally detected for chiral nematic liquid crystals under an applied external field. In contrast with these results, we found Skyrmion tube configurations in achiral nematics either with external fields or in their absence. When a sufficiently high external field is present, only a type of Skyrmion tube survives and it gives rise to an axially-symmetric localized configuration immersed in the standard nematic phase.
Thus, we reached a twofold target. On one hand this work presents a self-contained study about the formation and control of Skyrmion tubes under external fields and, more specifically, coaxial external fields. On the other hand it represents a first stone towards a general 3D study of these structures also including arbitrary orientations of the external field.
As a conclusion, the proposal of considering higher order free-energy expansions, as opposed to higher derivative ones [18], leads to interesting new perspectives in the liquid crystal science with many potential technological applications, where Skyrmion tubes might play an important role. On the theoretical side, according to our results it is clear that this kind of configurations emerge in a natural and straight way from the proposed energy and can be controlled by external fields.
As for future work, we plan to study the stability of the found solutions, the mutual interaction of Skyrmion tubes, the space arrangement of two or more of them and their lattice configurations. Furthermore, we aim at studying electro-optical effects and exploring other types of localized objects. Moreover, we also aim at exploring the effect of the compression of pseudolayers by using an appropriate compression energy and a representation of the director in terms of the geometric objects defining the layer surfaces as discussed at the end of Sec. II. Finally, we are also interested in studying liquid crystals confined within specific geometries and modeled by the quartic free-energy density along with its coupling with external fields. In this Appendix we collect the basic main functions and coefficients appearing in the equilibrium equations for Skyrmion tubes.
As for G 1 where As for G 2 G 2 = G 2 (r, f ) = g 21 + g 22 cos(2f ) + g 23 cos(4f ) + g 24 cos(6f ), where The function G 3 is given by where Finally, where Notice that in addition to (A7) other identities exist. This is due to the fact that there are seven free parameters: six independent elastic constants k i and a parameter β. Similarly to the previous case, we find that c π reads in order to fulfill the correct boundary conditions, i. e., f a (0) = π and f a (r → ∞) = θ 0 . Also in this case the best fitting procedure is successful, although for H = H cr the least square minimization must be performed for both the real and the imaginary part of f a simultaneously. More specifically, for every r both quantities Re[f (r) − f a (r)] and Im[f (r) − f a (r)] are taken into account when minimizing the sum of the squared differences. Indeed, the minimization of only the real part yields values of the best fitting parameters such that | − 1 + c π r 2 2 s(r)| > 1. As it can be noted also from the other cases, the ability of f a to fit f is weaker in the proximity of the bump through which the profile function reaches its asymptotic value. This is particularly true approaching the critical value of the external field. Thus, the values of a, c, d in this latter case are the best fitting ones when a reality condition is imposed to f a . The results are depicted in Fig. 12. The values of the optimal a, c, d for the cases taken into consideration are reported in Table II. In conclusion, our analysis leads to providing a good (within a few percentage) global approximated analytical expression of the Skyrmion tubes by using only 5 constants: θ 0 , a, c, d and β. The physical meaning of θ 0 and β is straightforward and they can be measured by suitable experiments. On the other hand, the remaining constants provide information about the shape of the Skyrmion tube which can be also experimentally observed.