Change in the magnetic configurations of tubular nanostructures by tuning dipolar interactions

We have investigated the equilibrium states of ferromagnetic single wall nanotubes by means of atomistic Monte Carlo simulations of a zig-zag lattice of Heisenberg spins on the surface of a cylinder. The main focus of our study is to determine how the competition between short-range exchange (J) and long-range dipolar (D) interactions influences the low temperature magnetic order of the nanotubes as well as the thermal-driven transitions involved. Apart from the uniform and vortex states occurring for dominant J or D, we find that helical states become stable for a range of intermediate values of γ = D/J that depends on the radius and length of the nanotube. Introducing a vorticity order parameter to better characterize helical and vortex states, we find the pseudo-critical temperatures for the transitions between these states and we establish the magnetic phase diagrams of their stability regions as a function of the nanotube aspect ratio. Comparison of the energy of the states obtained by simulation with those of simpler theoretical structures that interpolate continuously between them, reveals a high degree of metastability of the helical structures that might be relevant for their reversal modes.


Introduction
When planar structures are extended to 3D by changing their curvature, the magnetic behavior of the corresponding surface can change significantly, leading for example to the appearance of new magnetic configurations that are not energetically favorable in the planar case 1, 2 and also to the induction of anisotropic and chiral effects that merely originate from the curvature of the surface [3][4][5][6][7] . A nice example of these, are tubular magnetic nanostructures that have gained recent interest because the increased functionality arising from their inner and outer surfaces may contribute to improve existing applications 1, 8 based on other proposals. This fact makes them very attractive for diverse applications in several fields involving biomedicine 9, 10 , spintronics 11 , magnetic sensors or high-density magnetic storage devices 12 among others. The hollow nature of nanotubes favors the formation of closed flux, core free configurations. These kind of states becomes particularly relevant for magnetic storage purposes because they do not produce stray fields (flux-closure configuration) avoiding the consequent leakage of magnetic flux spreading outward from the tube. Additionally, the core-free aspect offers advantages over nanowires, as it allows more controllable reversal processes, increased stability 13 , and faster domain wall propagation 14,15 .
The rapid advance in fabrication techniques available nowadays allows to synthesize tubular nanostrutures by following different routes like chemical or electro-deposition 16,17 , hydrogen reduction 18 , sol-gel 19 , electrochemical synthesis 20 , Kirkendall effect 21 , self-assembly of NPs 22 and others 23 with a high degree of control on their geometrical parameters. In particular, different realizations of magnetic nanotubular structures have been achieved in the form of rolled up magnetic nanomembranes 24 , magnetite nanotubes on MgO wires 25 , nanoparticles deposited on the surface of nanotubes 26 , or carbon nanotubes either filled with magnetic nanoparticles (NP) [27][28][29] or with NP and networks of molecular molecules or nanoparticles attached to their surface 30,31 , only to mention some examples.
Vortex states having the magnetic moments circulating around the symmetry axis are expected to be formed in structures having cylindrical symmetry 32 . As predicted by theoretical studies, vortex and helical states have been found giving rise to a rich magnetic phase diagram where both temperature and geometry (aspect ratio) play an important role 33,34 . However, since their magnetization is close to zero, standard magnetic characterization precluded their experimental evidence until the advent of more sensitive techniques allowed their observation in individual nanocylinders, nanotubes and dot-shaped elements of different compositions by means of electron holography 35,36 , cantilever-based torque magnetometry [37][38][39][40] , nanoscale torque-SQUID magnetometry 41 , spin-polarized scanning tunneling spectroscopy 42 , X-ray microscopy 43 and MFM 44 . Different theoretical and simulation approaches have been used in order to study the thermodynamic properties of magnetic hollow nanocylinders, including micromagnetic studies 15,45 , scaling techniques 46,47 , many-body Green's functions [48][49][50] , ab-initio calculations [51][52][53] , and Monte Carlo (MC) simulations 54 .
Most of these works, have focused on the study of magnetic configurations as a function of the geometric parameters of the nanotubes for given values of material parameters. However, the occurrence of such phases and their thermal stability becomes also determined by the competition between the short-range exchange and long-range magnetic dipolar interactions 55,56 , quantified by γ. The value of γ can be tuned not only by changing parameters of the material (J and µ) , but also by changing the lattice spacing between the magnetic entities (magnetic ions or nanoparticles) and their spatial arrangement, since both contribute to the dipolar interaction energy. Real systems in which this can be achieved are, for example, carbon nanotubes used as structural templates or guides to attach or decorate their surface with high-moment magnetic nanoparticles, clusters or single-molecule magnets 57 .
Therefore, the present work will focus on studying the changes in the equilibrium configurations with γ considering first a fixed tube geometry and then extending the calculations to a range of tube radii and lengths. To address the problem, we use standard Metropolis MC to study the finite-temperature properties, and we revert to analytical or numerical calculations when comparing the properties of ground-state configurations obtained from simulations to theoretical ones. It turns out that the different equilibrium states with well-defined chiral order have magnetization close to zero and, therefore, this magnitude is no longer a valid order parameter to characterize the transitions between the mentioned states. We will show that these states and the transitions in between can be correctly identified and characterized by introducing the rotational of the magnetization as a measure of vorticity or circularity of the moments around the tube. This, combined with a detailed inspection of the equilibrium spin patterns, will allow us to determine the magnetic phase diagrams of their stability regions as a function of the nanotube aspect ratio.
The outline of the manuscript is a follows. We start by presenting the model for the nanotubes and computational details of the simulations. Then, the results of the equilibrium properties are presented for a (8,15) tube, followed by a detailed analysis of the ground state configurations. In the next subsection, we extend the calculations to tubes with different radii and lengths and present the corresponding phase diagrams as a function of the parameter γ. Finally, we analyze the metastability of the helical states by comparing the simulation results to analytical calculations based on two kind of states that can be theoretically parametrized. We finish with a discussion and final conclusions.

Model and Computational details
Our model system to simulate single-wall nanotubes consists of Heisenberg classical unitary spins S = (S x , S y , S z ) placed on the surface of a cylinder of radius R and finite length L 56 . The structure is formed by rolling a planar squared spin lattice with lattice parameter a along the (1,1) direction onto a cylindrical geometry in order to get a zig-zag ended tube that can be characterized by N z layers of rings stacked along the z axis with N spins per layer. For a given N, the tube is characterized by the angle between consecutive spins in the ring ϕ N = 2π/N , interlayer separation l = √ 2a/2 and radius R = N √ 2 2π a. Therefore, the distance between two consecutive spins in the ring is given by d 1 = R 2(1 − cos ϕ N ) = 2R sin(ϕ N /2) (see Supplementary  Fig. S1), and that between nn in adjacent rings is a. In what follows, we will designate the tubes by the notation (N, N z ).
Of course, there are other ways to arrange spins on a regular lattice on the surface of a cylinder (for example, varying the interlayer separation and the kind of stacking in between layers), but not all of them will be stable from the point of view of dipolar energy alone 53 . In order to justify the election of zig-zag tubes for our study, we have computed the dipolar energy of tubes with fixed (N, N z ) and a spin configuration forming a perfect vortex. In Fig. 1 (a), we show the energy landscape as a function of l and φ , the relative angle between consecutive layers. The dipolar energy presents minima at φ = ϕ N /2 and maxima at φ = 0, π/4, π/2 independently of l. This observation indicates that zig-zag tubes with spins in a vortex configuration are minimum energy states, whereas the formation of tubes with AA stacking (consecutive rings staked onto each other, with spins forming columns) is not energetically favored by dipolar interactions. Next, we fix l to the minimum possible separation between rings and we allow the polar angle θ of the spin orientation to vary from parallel [ferromagnetic (FM) configuration, θ = 0] to perpendicular [vortex (V) configuration, θ = π/2] to the tube axis. The computed energy landscape shown in Fig.  1(b) demonstrates that, depending on the degree of vertical misalignment between spins in consecutive rings, different spatial arrangements of the spins are favored by the dipolar interactions. Whereas for AA stacked tubes, alignment along the tube axis is preferred, for zig-zag tubes, a vortex configuration (θ = π/2) has the minimum possible dipolar energy among all considered kinds of ordering.
The spin interaction Hamiltonian can be written as where E ex stands for an isotropic short-range exchange coupling between nearest neighbors (nn) spins: being J i j = J the positive exchange constant accounting for ferromagnetic interactions. E dip stands for the long-range magnetic dipolar interaction given by where summation expands over all the set of spin pairs, D = µ 0 4π µ 2 /a 3 is the dipolar coupling constant, with µ 0 and µ the magnetic permeability of vacuum and the magnetic moment per spin respectively, and r i j is the relative position vector between i and j spins. We define the parameter γ = D/J in order to quantify the degree of competition between the short-range and long-range interactions. Temperature (T ) is expressed in reduced J/k B units where k B is the Boltzmann constant. The Hamiltonian does not contain a magnetocrystalline anisotropic energy term as so far we are interested in soft magnetic materials where the associated energy can be neglected compared to the exchange and magnetostatic counterparts 58 . The total energy per spin can then be written as We have conducted standard Metropolis MC simulations to analyze the temperature dependence of thermodynamic observables and low temperature spin configurations arising from the competition between dipolar and exchange interactions. In order to be able to reach low temperature configurations close to the vicinity of the thermodynamic equilibrium for a given γ, we have employed the following protocols for the thermal dependence and spin updates: (1) linear decrease of temperature with random homogeneous spin trials; (2) linear decrease of temperature with spin trials restricted to lie inside a cone with aperture varying with temperature; and (3) simulated annealing with a temperature power law decay of the form T = T 0 α k where spin updates are also restricted to lie inside a cone. For the first two, initial and final temperatures are 50 J/k B , 0.2 J/k B respectively, with a temperature step of 0.2 J/k B . For protocol (3) simulations are started at T 0 = 50 J/k B and ended at 0.113 J/k B after 200 temperatures with α = 0.97. The number of MC steps used at each temperature is 5 × 10 3 , with 2 × 10 3 used for averaging. Error bars in subsequent results have been obtained after averaging over 5 independent simulations starting with different seeds.
In order to generate spin trials inside a cone, we generated a random vector v i inside a sphere of radius S max which is added to the current spin vector S i and then the resulting vector S i is normalized, thus S i = v i + S i . In order to keep the acceptance  rates high, we use an adaptive scheme for the spin updates by reducing the aperture of the cone progressively as T decreases. We initially set S max = 2 at high T so that initially the trial spin can reach any state within the Bloch sphere and then, at every temperature, we replace it by the maximum value accepted at trials during the previous temperature.

Thermal dependence of near-equilibrium properties
Simulation results for the thermal dependence of energy, specific heat, magnetic susceptibility and magnetization per site of a (8,15) nanotube are shown in Fig. 2 for a range of γ values between 0 and 0.082. Energy curves shown in panel (a) display a concavity change pointing to a well-defined thermal-driven phase transition for all the values of γ considered. For a given temperature, the energy decreases as γ is increased independently of the temperature. In particular, for γ ≤ 0.03, the magnetization curves in Fig. 2(c) are smooth and well-behaved accounting for a well-defined ferromagnetic (FM) to paramagnetic (PM) transition above some pseudo-critical temperature T c , below which the spins tend to be aligned along the tube axis. This is also signaled by the peak in the specific heat curves of Fig. 2(b) that appear rounded due to finite-size effects. Similar features can be found in the thermal dependence of the susceptibility as observed in Fig. 2 where we observe that the pseudo-critical temperature for the PM to FM transition (signaled by the position of the peak in the susceptibility curve) is shifted to higher values when dipolar interactions are turned on.
However, for γ ≥ 0.04, the magnetization tends to zero when decreasing T and it is no longer a valid order parameter since, as we will show latter, the dominant dipolar interactions favor the appearance of closed flux states in the form of a vortex circulating around the tube axis. In order to better characterize the onset of these vortex states, we have found convenient to define a new order parameter that we will call vorticity that characterizes the degree of circularity of a given magnetic structure, which is defined locally as the vector: Correspondingly, the average z component of the vorticity of the whole tube can be written as a discrete sum over lattice points: 6 Figure 3. Thermal dependence of the vorticity (a) and the associated susceptibility (b) for different values of γ for a (8,15) tube. ρ has been normalized to the maximum attainable value. Error bars correspond to averages over 5 independent runs.
where ∆x k and ∆y l are components of interparticle vector connecting nn spins and the sum is over indexes that characterize lattice coordinates. It can be demonstrated that ρ is maximum for a perfect vortex configuration with a value that depends on the tube radius and length. In analogy to the susceptibility associated to magnetization, we also define an effective susceptibility to characterize transitions associated with the appearance of states with finite vorticity: As shown in Fig. 3(a), in the regime γ ≥ 0.04, the normalized vorticity behaves as expected for an order parameter. It attains a finite value at low temperature that saturates to a value that corresponds to a perfect vortex (V) state having the spins pointing tangent to the tube surface and perpendicular to the tube axis. The transition temperatures T c for FM-PM and V-PM transitions were extracted by averaging the locations of the maxima of the corresponding susceptibilities [see Fig. 3(b)] and specific heat and their dependencies on γ are shown in Fig. 4. Results indicate that T c increases with γ for both transitions as the dipolar interaction becomes stronger. However, for the V-PM transition, T c seems to saturate when the dipolar energy dominates over exchange, in agreement also with finite size scaling theory 59 since a γ increase can also be related to a size increase.
For intermediate values 0.03 < γ < 0.04, the magnetization begins to exhibit a noisy behavior with big error bars [see Fig.  2(d)], although characterized by non-vanishing values, despite of the several configurational averages performed. Concomitantly, the vorticity [see Fig. 3(a)] behaves less noisy and acquires finite values at low temperatures, which means that certain degree of circulation of the spins around the tube surface emerges. Thus, the low temperature configurations are characterized by helical (H) states 56 having spins pointing in the plane tangent to the tube surface as for the V states, but now with some misalignment relative to the tube axis.
For this intermediate range of γ, it was not possible to infer precisely a H-PM transition temperature due to the high degree of metastability. The reason for this can be clearly noticed in Fig. 5, where several long runs using the previously mentioned protocols 1 and 2 led to different low temperature states through different paths along the energy landscape, a situation that does not happen for lower or higher γ values. Therefore, in order to study the low temperature configurations, we ran additional simulations using simulated annealing protocol 3 varying the temperature according to a T = T 0 0.95 k dependence and sampling trial spin movements in a cone with progressive reduction of its aperture in order to ensure convergence to the lowest energy states.
By using the vorticity as an order parameter, we were able to obtain the phase diagram shown in Fig. 6 for two different tube sizes. The phase diagram evidences an increase in the vorticity of the system as the strength of the long range dipolar interaction increases. For both sizes, pure vortex states are stable at the highest considered values of γ, and this is the reason why γ values greater than 0.06 were not explored. Finally, we point out that a finite-size effect can be observed on the phase diagram. With a slight variation of the tube radius and length from (8,15) to (7,14), the onset of stable vortex states takes place at smaller γ. We will further elaborate on this point in a subsequent subsection.

Ground state characterization
In order to better characterize magnetic ground state configurations, we have expressed the total magnetization of the tubes in cylindrical coordinates (azimuthal m φ tangential to the tube surface, radial m ρ and longitudinal m z ). By projecting the spin    Figure 6. Dependence of the normalized vorticity ρ z on γ calculated for the low temperature configurations obtained using simulated annealing for tubes of size (8,15) and (7,14). The continuous lines serve to define a phase diagram separating the regions corresponding to to FM, Helical and Vortex states, whose corresponding magnetic configurations are shown in the insets for a (7,14) tube. Blue spheres correspond to lattice nodes and red arrows to spin directions.
vectors along a local reference frame centered on the spin positions r i = (R cos ϕ i , R sin ϕ i , z i ) (see inset in Fig. 7) the spin vector components can be written as: This indicates that H states can be seen as a combination between V and FM alignment that originates from the competition between the exchange interaction trying to align spins along a single direction and dipolar interactions that favor circularity with spins tending to point tangent to the surface of the tube to reduce magnetostatic energy. These helical states can be better characterized by plotting the angle θ averaged over all the θ i of spins in every layer, i.e. for each z value along the length L of the tube, as shown in Fig. 8. An alternative representation is given in Supplementary Fig.  S2, where all the spins of the tube are represented by dots on a unit sphere according to their orientation and are given a color that depends on the value of γ of the considered configuration. As it can be seen in Fig. 8, the angle θ is not uniform along the tube length and tends to increase towards the tube ends for all the computed γ values. Only for γ = 0 the vorticity is strictly zero, and θ is constant along the tube. For γ 0.3, spins are still mostly aligned along some common θ angle which deviates from 0 (FM states) with increasing dipolar interactions , while for γ 0.4, all configurations average around θ ≈ π/2. In the unit sphere representation of Supplementary Fig. S2, FM states appear as a groups of dots localized around the poles, while V states are revealed in the form of groups distinguishable along the equator of the sphere (note that the amount of groups equals twice the number of spins per layer due to geometrical considerations and the discreteness of the lattice).
However, for intermediate values in the region 0.03 < γ < 0.04 where H states arise, profiles are not longer flat evidencing different spin configurations at the edges and the inner part of the tube. More concretely, θ is lower in the middle part of the tube and greater at the edges. This means that the degree of vorticity is higher at the ends than at the central part of the tube. This behavior can be more clearly observed for longer tubes, as we will show in the next section. Due to the competition between exchange and dipolar energies, helical states transform onto vortex states as γ increases, and this change occurs first on the rings near the tube edges.

Size dependence and phase diagram
In the previous subsections, we focused on determining the low energy configurations as a function of γ for a fixed tube radius and length, exemplified by a (8,15) tube. Now, we carry out a systematic study of the low temperature configurations for different values of N z , N at different fixed interaction strengths γ. For this purpose, we start first by varying the tube length while keeping the number of spins per layer at N = 8. Using simulated annealing, we have obtained the low temperature magnetization and vorticity of the low temperature configurations, as shown in Fig. 9 for a number of layers N z in the range from 4 to 25 and 0.015 ≤ γ ≤ 0.045.
Length driven transitions from FM to V states are identified by the N z for which m ≈ 1 and ρ z ≈ 0 and transitions from FM or H to V states by m ≈ 0 and ρ z ≈ 1 , while intermediate values of these parameters characterize regions where H states are favored. As observed in Fig. 9, transitions are sharp for γ < 0.03, reflecting the fact that, in this range, exchange interactions dominate and states with almost uniform alignment along the tube axis are favored. Only for tubes with low aspect ratio, formation of V states is observed. Starting at values γ > 0.03, transitions between ρ z = 0 and 1 become smoother and the region in between indicates the interval of N z for which H states appear. With increasing values of γ ≥ 0.04, V states prevail up to higher tube lengths, until a point where formation of FM states is not favored even for the longest simulated tubes.
As an example of what happens for longer and wider nanotubes, we show representative magnetization height profiles for γ = 0.05 in Fig. 10. For wider tubes (N = 16), a perfect V configuration is obtained for short tubes but this transforms into a H state when the length increases. However, even for a tube with aspect ratio 1 : 4, only the central layer of spins is oriented along the tube axis. Much longer tubes have to be considered in this case for the observation of an extended FM central region, due to the strong in-plane dipolar fields generated by the vortices formed at the tube ends.
We have extended the calculations to additional tube radii N=10, . . . , 16, using the procedure described previously in order to identify the N, N z at which FM-H and H-V transitions are observed. We have obtained the phase diagrams in the (N, N z ) plane shown in Fig. 11 for three values of γ. We realize that the V-H critical curve can be fitted to a linear law N z = A + BN (dotted line) with A and B values given in Table 1. As observed, when increasing γ, the slope of this line increases and the N axis intercept tends to zero, meaning that, even though if exchange is dominant, V states can be stabilized for low aspect ratio tubes and big enough radii. However, the line H-FM curve exhibits a more pronounced change with the radius of the tube. We have found that this line can be fitted to an exponential dependence of the form N z = N 0 + A 2 exp(N/B 2 ) with A 2 and B 2 values given in Table 1. Comparing the results for different γ, we see that there is always a narrower region for formation of H states for tubes having length equal to radius (∼ 12, 10, 8 for γ = 0.01, 0.015, 0.02) and that is shifted to lower N with increasing γ.  Fig. 11. Values for A and B are the fitted parameters to a linear law for the V-H curve and A 2 and B 2 are the fitted parameter to an exponential law for the FM-H curve.

Comparison to analytical calculations
Given the narrow region of γ values for which H configurations with complex magnetic order are stabilized, next we will study their metastability by comparing their total energy to that of simpler configurations that interpolate continuously between FM and V states.

Quasi-uniform states
As a first case, we have considered configurations with the spins pointing along the local plane tangent to the tube surface and all oriented along a direction characterized by a polar angle θ with respect to the tube axis. We start computing the exchange energy per spin in J units (i.e. total exchange energy divided by N tot = NN z ) , which can be shown to be given by (see Supplementary Information for details): where z is the coordination number (z = 4 in our case), ϕ N = 2π/N is the angle subtending two consecutive spins on the same layer containing N atoms, whereas ϕ N /2 is the basal projection angle between two n.n. interacting spins belonging to consecutive layers according to the zig-zag geometry. Thus, exact values can be obtained for the FM or V configurations by setting θ = 0, π/2 in Eq. 9, while intermediate values of θ would correspond to H states. Notice that the energy depends on the tube radius (or number of spins per ring N) through ϕ N and also on N z due to the finite tube length. In the long tube limit N z → ∞, it becomes ε ex = −(z/2) for the FM state and ε ex = −(z/2) cos(ϕ N /2) for the vortex state, as expected. Due to the sin 2 θ dependence, ε ex increases progressively as the system evolves from FM to V states passing through intermediate H states. The variation of this energy with θ is shown in Fig. 12 in blue circles for a (8, 15) tube.
The dipolar energy can also be computed using analytical expressions, although their derivation is more complex due to the long-range character of the interactions and the fact that contributions from spins in the same ring and from different (odd or even) rings have to be considered separately (see Supplementary Information for the details). As shown in Fig. 12 (green symbols), ε dip is minimized for θ = π/2 and increases when the magnetic moments progressively align along the tube axis.  (8,15) tube obtained from the simulations following an annealing (green squares) or thermalization (blue dots) procedures are compared to those of configurations having spins oriented at θ = 0 (FM state, black circles) or θ = π/2 (V state, red squares) tangentially to the tube surface. Yellow triangles correspond to analytical results for the mixed states with minimum energy as described in the text. Inset: zoom around the range 0.03 < γ < 0.04.

14/21
Therefore, as a consequence of the competition between exchange and dipolar energies, the state with minimum energy will depend on the value of γ. This can be clearly observed by looking at the total energy plotted in Fig. 12(a), (b) for two different values of γ: for small γ, the FM state with θ = 0 is favored, while for large γ the V state has the lowest energy.
By looking at the evolution of the total energy when going from low to high γ, it is clear that there must be an intermediate γ for which ε tot does not depend on the polar angle θ . This value can be found by equating the total energies per spin for θ = 0, π/2, which results in Inserting the values for the (8,15) tube (see Supplementary Information), this gives γ = 0.034 with a value for the energy ε tot /J = −2.31 which lies precisely at the center of the region of γ where H states are obtained in the simulations. We notice that, as can be seen in Fig. 12(c), for γ close to this critical value, configurations with different θ have energies that differ in less than 1% and, as a consequence, they are quasi-degenerated. Therefore, within this range of γ's, it is very difficult to ascertain if the configurations obtained from simulation are the true ground states, even if the annealing protocol is used. In fact, comparing the energies of the states obtained from the simulated annealing process with those of the states with uniform θ (see Fig. 13), we find that for γ < γ simulation results have exactly the same energy as the FM states whereas, for γ > γ , configurations obtained through simulated annealing have slightly higher energies than V states. Notice also that, in the same figure, we can see that configurations obtained after annealing (green squares) have indeed lower energies than those obtained after the usual thermalization process (blue circles) used to compute thermodynamic averages, as anticipated in the previous section.
Finally, we have computed the dependence of γ on N z as shown in Supplementary Fig. S3 for different tube radii. We observe that the region of γ for which quasi-degenerated configurations are expected, is displaced to higher values of γ as the tube length increases and that, for a fixed tube length, γ decreases when increasing the tube radius. This serves also as an indication for the expected location of the region favoring H states and it is in agreement with what is found on the phase diagrams presented in Fig. 11.

Mixed states
The low temperature H configurations that we have obtained by simulations present some non-uniformities at the tube ends that depart from the states with uniform θ considered in the previous subsection. Therefore, in an effort to find a better approximate analytical description for the equilibrium magnetization patterns, we will consider now the energy of configurations for which a central region (of length l − 2λ ) of the nanotube with a uniform FM state is connected by a pair of incomplete vortex walls that extend towards the two ends of the tube, where the spins form an angle θ 0 with respect to the tube axis 15 . An example of such a configuration is displayed in Supplementary Fig. S4.
For this kind of configurations, no analytic expressions can be obtained for the energies. Therefore, for a given tube with length L and radius R, we have calculated numerically energies for a selection of values of θ 0 ∈ [0, π/2] and λ ∈ [0, L/2] and located the minima of the E(θ 0 , λ ) surface for a range of γ values. The results of this minimization procedure are shown in First, let us stress the perfect linear dependence of E min on γ that indicates a linear scaling of the energy on this parameter when considering the values of θ 0 , λ that minimize the energy. Second, regarding the angle at the ends of the tubes and the central region with uniform magnetization, we observe that, for γ < γ , θ 0 ≈ 0, λ ≈ L and, therefore, FM configurations are obtained independently of the tube length. For γ > γ , θ 0 (γ) follows the same trend independently of the tube length, increasing progressively towards 90 • since for higher γ the dipolar interaction favors flux closure at the ends. For the shortest considered tubes (N z ≤ 15), the vortex region extends to the central part of the tubes, whereas for long tubes [see for example the case N z = 45 in Fig. 14 (c) ] this region progressively expands as γ is increased.
Finally, we have compared the energies of these optimized mixed states to those of the configurations from our MC simulations. As indicated by the yellow triangles in Fig. 13, these states have lower energy than FM ones for γ > γ but clearly higher energies than those obtained by both simulation procedures, which are closer to the energies of perfect V states. This result is in contrast with the work of Landeros et al. 15 , who obtained good agreement between simulations and optimized mixed states. However, they considered tubes with AA stacking of the spins instead of zig-zag stacking as in our case, and this may make a difference concerning magnetic properties. To confirm this, we have checked (not shown) that for AA tubes of the same size as ours, the energies of the optimized mixed states are indeed lower than perfect V states. This means that the spatial distribution of the spins on the surface of the nanotube has an influence on their magnetic order due to the anisotropic and the long-range character of the dipolar interaction.

Discussion and Conclusions
Using atomistic MC simulations, we have studied the ground state configurations that arise as a result of the competence between exchange and dipolar interactions in magnetic nananotubes of finite size, establishing radius-length phase diagrams for their stability for different γ's. We have shown that, for a given radius and length, ferromagnetic (FM), vortex (V) or helical (H) states can be stabilized, depending on the value of γ. Our results are in agreement with experiments that have recently reported the observation of these states in real samples. We would like to emphasize that, although experimental studies have been conducted on nanotubes which have considerably bigger dimensions than those considered here, the application of scaling techniques 46,47,60 , would bring the range of γ for the observation of the different states towards the order of magnitude of real material sizes .
As mentioned previously, in our model Hamiltonian we did not include magnetocrystalline anisotropy as we were aiming to simulated magnetically soft materials. However, the results of previous studies 61-64 by means of micromagnetic simulation of nanotubes with finite width considering uniaxial or cubic anisotropy, can give a hint on the effect of its inclusion in our results. Therefore, we speculate that inclusion of uniaxial anisotropy along the tube axis in our model would lead to an increase of the length of the central FM region and a reduction of the region λ where vortices are formed at the tube ends. Probably, this would lead also to a reduction of the θ 0 angle. Regarding the phase diagrams of Fig. 11, this would translate in displacement of the FM and V stability regions towards higher N values. However, it is difficult to guess in what way the transition lines between the three regions would be affected by the anisotropy. We are planning to address these issues in a forthcoming publication.
Furthermore, we would like to stress that our results are applicable to a wide variety of magnetic systems independently of the physical realization of the spins that reside on the nanotube surface. They could represent magnetic moments of atomic spins or of a magnetic cluster or nanoparticle as this would only change the values of D and J. As shown in Supplementary  Information, typical values of D are in the range from 10 −2 meV for atomic spins to 40 meV for nanoparticles while J vary in between 1 and 10ths of meV for most FM materials. Therefore, materials may be easily found with γ's falling within the interval 10 −2 − 10 −1 where we have found H configurations as ground states.
We have shown that a modified simulated annealing protocol has to be employed in order to be certain that the obtained configurations are indeed the true ground states, specially the ones showing helicoidal order. In order to determine the interval of γ values for which H states are the quasi-equilibrium configurations for a given radius and length of the nanotube, we have found crucial to introduce the so-called vorticity order parameter which clearly identifies the formation of states with zero magnetization but a well defined circularity or chirality. This has allowed us to establish the stability regions for their formation in radius-length phase diagrams, showing that, with increasing γ, the region of stability of H states becomes wider for higher aspect ratio tubes.
In an effort to check the robustness of the helicoidal states, we have compared the energies obtained by simulation to those of similar configurations that can be expressed analytically, such as quasi-uniform states with spins having a common angle with respect to the tube axis or mixed states formed by a FM ordered regions connected to vortex states at the tube ends. We have found that the states resulting from MC simulated annealing have indeed lower energies than any of the considered analytical ones for the range of γ's where the H states appear. However, the energy differences in this range are really small, indicating the high degree of metastability of the H states. This fact can be relevant for the nanotube reversal modes in hysteresis loops, which we will study in a future publication.