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

Abstract

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¹¹, magnetic sensors or high-density magnetic storage devices¹² 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¹³, 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¹⁸, sol-gel¹⁹, electrochemical synthesis²⁰, Kirkendall effect²¹, self-assembly of nanoparticles (NP)²² and others²³, 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²⁴, magnetite nanotubes on MgO wires²⁵, nanoparticles deposited on the surface of nanotubes²⁶, or carbon nanotubes either filled with magnetic NP, NP attached to their surface^27,28,29 or networks of molecular magnets^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³². 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⁴¹, spin-polarized scanning tunneling spectroscopy⁴², X-ray microscopy⁴³ and MFM⁴⁴.

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⁵⁴.

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⁵⁷.

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 simulations 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.

Results

Model and Computational details

Our model system to simulate single-wall nanotubes consists of Heisenberg classical unitary spins $\overrightarrow{S}=({S}_{x},\,{S}_{y},\,{S}_{z})$ placed on the surface of a cylinder of radius R and finite length L⁵⁶. 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=\sqrt{2}a\mathrm{/2}$ and radius $R=\frac{N\sqrt{2}}{2\pi }a$. Therefore, the distance between two consecutive spins in the ring is given by ${d}_{1}=R\sqrt{\mathrm{2(1}-\,\cos \,{\phi }_{N})}=2R\,\sin ({\phi }_{N}\mathrm{/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⁵³. 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.

Figure 1

Dipolar energy landscape for a (8, 15) tube (φ₁ = π/4) with spins aligned along the surface tangent to the tube: (a) as a function of ϕ, the angle between nn in consecutive rings, and interlayer separation l, (b) as a function of θ (their common angle with respect to the tube axis) and ϕ, when the interlayer separation is fixed at $l=\sqrt{2}a\mathrm{/2}$.

The spin interaction Hamiltonian can be written as

$$ {\mathcal H} ={E}_{ex}+{E}_{dip},$$

(1)

where E_ex stands for an isotropic short-range exchange coupling between nearest neighbors (nn) spins:

$${E}_{ex}=-\,\sum _{\langle i,j\rangle }{J}_{ij}{\overrightarrow{S}}_{i}\cdot {\overrightarrow{S}}_{j},$$

(2)

being J_ij = J the positive exchange constant accounting for ferromagnetic interactions. E_dip stands for the long-range magnetic dipolar interaction given by

$${E}_{dip}=D\sum _{i < j}(\frac{{\overrightarrow{S}}_{i}\cdot {\overrightarrow{S}}_{j}-\mathrm{3(}{\overrightarrow{S}}_{i}\cdot {\hat{r}}_{ij})({\overrightarrow{S}}_{j}\cdot {\hat{r}}_{ij})}{|{\overrightarrow{r}}_{ij}{|}^{3}}),$$

(3)

where summation expands over all the set of spin pairs, $D=\frac{{\mu }_{0}}{4\pi }{\mu }^{2}/{a}^{3}$ is the dipolar coupling constant, with μ₀ and μ the magnetic permeability of vacuum and the magnetic moment per spin respectively, and ${\overrightarrow{r}}_{ij}$ 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⁵⁸. The total energy per spin can then be written as

$${\varepsilon }_{total}/J={\varepsilon }_{ex}+\gamma {\varepsilon }_{dip}\mathrm{.}$$

(4)

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₀α^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₀ = 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³, with 2 × 10³ 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 $\overrightarrow{{v}_{i}}$ inside a sphere of radius S_max which is added to the current spin vector $\overrightarrow{{S}_{i}}$ and then the resulting vector $\overrightarrow{{S^{\prime} }_{i}}$ is normalized, thus $\overrightarrow{{S^{\prime} }_{i}}=\overrightarrow{{v}_{i}}+\overrightarrow{{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.

Figure 2

Thermal dependence of energy (a), specific heat (b,c), magnetization (d) and susceptibility (e) for different values of γ for a (8, 15) tube. Error bars correspond to averages over 5 independent runs.

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(d), 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:

$${\overrightarrow{\rho }}_{i}=\overrightarrow{\nabla }\times {\overrightarrow{S}}_{i}.$$

(5)

Correspondingly, the average z component of the vorticity of the whole tube can be written as a discrete sum over lattice points:

$$\langle {\rho }^{z}\rangle =\frac{1}{N{N}_{z}}\sum _{k,l}(\frac{({S}_{k+\mathrm{1,}l}^{y}-{S}_{k,l}^{y})}{{\rm{\Delta }}{x}_{k}}-\frac{({S}_{k,l+1}^{x}-{S}_{k,l}^{x})}{{\rm{\Delta }}{y}_{l}}),$$

(6)

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:

$${\chi }_{eff}=\frac{\langle {({\rho }^{z})}^{2}\rangle -{\langle {\rho }^{z}\rangle }^{2}}{{k}_{B}TN}.$$

(7)

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⁵⁹, since a γ increase can also be related to a size increase.

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.

Figure 4

Transition temperature as function of γ for a (8, 15) tube. Region I stands for FM-PM transitions whereas region III stands for V-PM transitions. In region II, H-PM transitions are present and the system is highly metastable. Error bars correspond to averages over 5 independent runs.

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⁵⁶ 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.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.

Figure 5

Energy dependence with both the number of MC steps and the magnetization per site for an intermediate γ = 0.035 value where a H-PM transition is expected. Lines referred as 3D in the legend correspond to the cooling protocol (1) mentioned in the main text, whereas those referred as Cone correspond to protocol (2). Results, revealing a high metastability, were obtained by using several initial random seeds (different colors) and low temperature T = 0.1, J/k_B.

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.

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.

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 vectors along a local reference frame centered on the spin positions ${\overrightarrow{r}}_{i}$ = (Rcosφ_i, Rsinφ_i, z_i) (see inset in Fig. 7) the spin vector components can be written as:

$$\begin{array}{c}{S}_{\rho i}=\,\sin \,{\theta }_{i}\,\cos \,{\phi }_{i}\\ {S}_{\phi i}=\,\sin \,{\theta }_{i}\,\sin \,{\phi }_{i}\\ {S}_{zi}={S}_{zi}=\,\cos \,{\theta }_{i},\end{array}$$

(8)

so the vortex and helical states formed at intermediate γ values can be better described. The magnetization components of the configurations obtained at the lowest temperature after simulated annealing are plotted in Fig. 7 for different values of γ. For γ > 0.04, the configurations present vanishing radial and longitudinal magnetizations (i.e. m_ρ ≈ 0 and m_z ≈ 0) and a tangential component m_ϕ ≈ 1, characteristic of a V state perpendicular to the z axis. In contrast, for γ < 0.03 the obtained configurations have vanishing radial and tangential magnetizations (i.e. m_ρ ≈ 0 and m_ϕ ≈ 0) and a longitudinal component m_z ≈ 1, characteristic of FM alignment along the tube axis. However, in the region 0.03 < γ < 0.04 [see inset in Fig. 7(a)], configurations have only vanishing radial magnetization, while both m_z and m_ϕ have non-vanishing values that decrease (increase) with increasing γ. 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.

Figure 7

(a) Cylindrical components of the total magnetization per site of the low temperature magnetic configurations obtained by simulated annealing as a function of the γ parameter. (b) Zoom of the intermediate region where transition from ferromagnetic alignment to vortex states takes place. The inset in (a) shows the coordinates used to describe the magnetic configurations.

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.

Figure 8

Height profile of the average polar angle 〈θ〉 for a (8, 15) tube. Results for intermediate values of γ are given in the left panel and while in the right panel results are shown for the highest γ values considered in the simulations. In the right panel, 〈θ〉 ≈ 30° stands for a FM alignment whereas values of 〈θ〉 close to 90° correspond to pure vortex states.

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 $\gamma \lesssim 0.3$, spins are still mostly aligned along some common 〈θ〉 angle which deviates from 0 (FM states) with increasing dipolar interactions, while for $\gamma \lesssim 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.

Figure 9

Magnetization (a) and z– component of the vorticity (b) of the equilibrium configurations as a function of the number of layers N_z of a tube with N = 8 for γ ∈ [0.015, 0.040] as indicated in the legend.

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. Snapshots of the spin configurations are reported in Supplementary Fig. S5. As it can be observed, for the thinner N = 8 tubes, increasing the length leads to an increase of the central part with FM aligned spins at the expenses of the decrease of the vortices at the tube ends and also of the width of the domain walls connecting the vortices to the central FM region.

Figure 10

Typical longitudinal profile of the average polar angle 〈θ〉 for γ = 0.05 nanotubes with different radii and lengths as indicated in the legend. FM alignment in the central part of the tube and a trend to form vortices at the edges is observed as the tube length is increased.

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₀ + A₂exp(N/B₂) with A₂ and B₂ 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 γ.

Figure 11

Configuration phase diagrams for tubes with different aspect ratios (N/Nz) and three values of γ = 0.01, 0.015, 0, 02. Red circles identify FM-H transitions and green stars V-H transitions. Continuous lines are fits of the transition curves to the laws described in the text.

Table 1 Results of the fits of the curves separating the V-H and FM-H configurations shown in Fig. 11.

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):

$${\varepsilon }_{ex}=-\frac{z}{2}(1-\frac{1}{{N}_{z}})[1-{\sin }^{2}\theta \mathrm{(1}-\,\cos \,({\phi }_{N}\mathrm{/2))}],$$

(9)

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}\theta$ 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.

Figure 12

Dependence of exchange, dipolar and total energy per spin on the polar angle θ for a configuration with spins tangential to the tube surface all forming an angle θ with respect to the tube axis. Results are shown for a (8, 15) nanotube. Panel (a) is for γ = 0.005 and panel (b) is for γ = 0.095. Panel (c) displays a contour plot of the dependence of the total energy on θ for γ values close to the formation of H states, where isoenergetic curves are plotted with black lines.

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 stars), ε_dip is minimized for θ = π/2 and increases when the magnetic moments progressively align along the tube axis. 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

$${\gamma }^{\ast }=\frac{{\varepsilon }_{ex}(\pi \mathrm{/2)}-{\varepsilon }_{ex}\mathrm{(0)}}{{\varepsilon }_{dip}\mathrm{(0)}-{\varepsilon }_{dip}(\pi \mathrm{/2)}}.$$

(10)

Inserting the values for the (8, 15) tube (see Supplementary Information), this gives ${\gamma }^{\ast }=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 $\gamma < {\gamma }^{\ast }$ simulation results have exactly the same energy as the FM states whereas, for $\gamma > {\gamma }^{\ast }$, 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.

Figure 13

Energies of ground state configurations of a (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.

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, ${\gamma }^{\ast }$ 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 θ₀ with respect to the tube axis¹⁵. 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 the energies numerically for a selection of values of θ₀ ∈ [0, π/2] and λ ∈ [0, L/2] and located the minima of the E(θ₀, λ) surface for a range of γ values. The results of this minimization procedure are shown in Fig. 14, where the γ dependence of the θ₀ and λ values that minimize the energy [panels (b) and (c)], and the corresponding energies E_min [panel (a)] are shown for tubes with N = 8 and different lengths. Similar results have been obtained for tubes with radi up to N = 16 (not shown).

Figure 14

Results of energy numerical minimization of analytical mixed states characterized by a FM central region of length λ and angle θ₀ at the tube ends. Dependence of (a) E_min, (b) θ₀ and (c) λ, and on γ is shown for tubes with a fixed radius N = 8 and different number of layers N_z = 8, 15, 20, 30, 45. The inset in panel (c) displays the energy surface for the (8, 15) tube with $\gamma ={\gamma }^{\ast }=0.034$.

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 θ₀, λ 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 $\gamma < {\gamma }^{\ast }$, θ₀ ≈ 0,λ ≈ L and, therefore, FM configurations are obtained independently of the tube length. For $\gamma > {\gamma }^{\ast }$, θ₀(γ) 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 $\gamma > {\gamma }^{\ast }$ 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.¹⁵, 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,62,63,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 θ₀ 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⁻² meV for atomic spins to 40 meV for nanoparticles while values of J vary in between 1 and 10 ths of meV for most FM materials. Therefore, materials may be easily found with γ’s falling within the interval 10⁻² −10⁻¹ 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.