1 Introduction
Spontaneous self-organization of magnetized plasma lies at the basis of many fascinating phenomena in both fusion reactor operation and astrophysical plasma observations. In such situations, magnetic helicity in the plasma is of crucial importance in determining the evolution of the system. Magnetic helicity is an integral quantity calculated by
$H_{m}=\int \boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}^{3}x$
where
$\boldsymbol{A}$
is the vector potential and
$\boldsymbol{B}=\unicode[STIX]{x1D735}\times \boldsymbol{A}$
is the magnetic field. Helicity was given its name by Moffatt (Reference Moffatt1969) who recognized its topological interpretation; a measure of the self- and interlinking of magnetic field lines. This notion was extended to non-closing and ergodic field lines by Arnol’d (Reference Arnol’d1986). In a perfectly conducting plasma, the magnetic field can be seen as being ‘frozen in’ and is advected with the fluid motion (Alfvén Reference Alfvén1943; Batchelor Reference Batchelor1950). Since the fluid motions can then only distort and reshape the magnetic field lines, but cannot break or cause unlinking, the conservation of helicity is easily understood from a topological perspective.
In perfectly conducting plasma helicity is exactly conserved (Berger & Field Reference Berger and Field1984), whilst in resistive plasma the rate of energy decay is strongly constrained by its presence (Del Sordo, Candelaresi & Brandenburg Reference Del Sordo, Candelaresi and Brandenburg2010). The most famous example where helicity determines a self-organizing process is the Taylor conjecture (Taylor Reference Taylor1974, Reference Taylor1986) which states that the magnetic field in a toroidally bounded plasma relaxes to a linear force-free state (shown by Woltjer (Reference Woltjer1958) to be the lowest-energy state) with exactly the same helicity as it started with.
But what state is achieved when a helical plasma relaxes in an environment without a boundary? If the plasma’s fluid pressure is high compared to the magnetic pressure (high plasma
$\unicode[STIX]{x1D6FD}=p/(B^{2}/2)$
), the linking in the initial field gives rise to a self-organized magnetohydrodynamic (MHD) equilibrium where field lines lie on nested toroidal magnetic surfaces (Smiet et al.
Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015). The approximately axisymmetric field is in a Grad–Shafranov equilibrium (Shafranov Reference Shafranov1966), and the Lorentz force is balanced by a gradient in pressure. The pressure is lowest on the magnetic axis, which is the field line lying at the centre of the nested toroidal magnetic surfaces. In this structure the rotational transform is nearly constant from surface to surface. As a consequence the magnetic field line structure is topologically similar to the mathematical structure of the Hopf fibration (Hopf Reference Hopf1931) or its generalization to torus knots (Arrayás & Trueba Reference Arrayás and Trueba2014; Smiet et al.
Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015). It should be noted that the Hopf structure has previously also been used in a beautiful paper by Finkelstein and Weil (Finkelstein & Weil Reference Finkelstein and Weil1978) to generate linked and knotted magnetic fields for astrophysics. The relaxation of the Hopf field to the Grad–Shafranov equilibrium has been shown using topology conserving relaxation in our recent paper (Smiet, Candelaresi & Bouwmeester Reference Smiet, Candelaresi and Bouwmeester2017). It is remarkable that this structure is obtained for a wide class of initially helical fields; it emerges from trefoils (Smiet et al.
Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015), twisted rings and even Borromean linked flux tubes (Candelaresi & Brandenburg Reference Candelaresi and Brandenburg2011).
The robust generation of this ordered magnetic structure makes it natural to assume that a similar process emerges after the twisted magnetic field of a coronal loop is ejected into the high-
$\unicode[STIX]{x1D6FD}$
interplanetary plasma of the solar wind. The pressure in the solar wind at 1 A.U. is of the order of
$p=1.4\times 10^{-11}~\text{N}~\text{m}^{-2}$
and the field strength
$B=6\times 10^{-9}T$
such that the plasma
$\unicode[STIX]{x1D6FD}=2\unicode[STIX]{x1D707}_{0}p/B^{2}\simeq 1$
(Goedbloed & Poedts Reference Goedbloed and Poedts2004). Such events are called coronal mass ejections (CMEs), and CMEs are correlated with the observation of magnetic clouds (Burlaga Reference Burlaga1991). A magnetic cloud is a localized magnetic structure in the interplanetary plasma with increased magnetic field strength, and where the direction of the field varies by a large angle (Burlaga Reference Burlaga1991). These magnetic signatures are observed by interplanetary satellites with increasingly accurate magnetic instruments (Raghav & Kule Reference Raghav and Kule2018). Unfortunately due to the low density of probes in the interplanetary medium, high-resolution measurement of the complete magnetic structure in these clouds is challenging.
There are several models for the magnetic structure of these clouds published in the literature (Burlaga Reference Burlaga1991). Some models assume the field in the cloud is still magnetically connected to the surface of the Sun, but there are several models that assume a localized magnetic field which is created by internal currents and balanced by the external plasma pressure (Ivanov & Harshiladze Reference Ivanov and Harshiladze1985; Burlaga Reference Burlaga1991; Vandas et al. Reference Vandas, Fischer, Pelant and Geranios1992; Garren & Chen Reference Garren and Chen1994; Kumar & Rust Reference Kumar and Rust1996). In this paper we posit the self-organized state identified in Smiet et al. (Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015) as a new model. This model gives different predictions for the structure from the models above.
In this paper, we study the evolution of the self-organizing equilibrium starting from a twisted flux tube. We vary both the resistivity of the simulation as well as the amount of twist in the initial flux tube. The evolution of the structure is governed by two processes. First, the lowering of the magnetic field strength changes the equilibrium condition, such that the depression in pressure becomes smaller. Second, finite resistance breaks the frozen-in condition, allowing the plasma fluid to slip perpendicular to the magnetic field lines of the configuration. The net effect of this is a fluid flow directed towards the region of lowest pressure. The combined effect of these two processes is a structure that grows on a resistive time scale. The magnetic topology, characterized by the rotational transform of the magnetic surfaces, quickly reaches a nearly flat profile with a slight positive curvature. The rotational transform decays according to a power law, the characteristic exponent of which depends on the aspect ratio of the structure.
Coordinate system used for the construction of the initial magnetic field. The surfaces through which the toroidal flux
$\unicode[STIX]{x1D713}_{t}$
and poloidal flux
$\unicode[STIX]{x1D713}_{p}$
are defined are shown by the blue and green circles respectively;
$\unicode[STIX]{x1D753}$
is the coordinate pointing in the poloidal direction of the torus, and
$\unicode[STIX]{x1D703}$
is the coordinate pointing in the toroidal direction. The magnetic axis is given by the red circle located at
$R^{\ast }=1$
.
2 Initial field
In our previous paper (Smiet et al. Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015), we showed how a self-organizing equilibrium is generated through the chaotic reconnection of linked magnetic flux rings. The violent reconnection and high spatial gradients generated when the flux tubes meet, necessitated a high value of resistivity and low magnetic field strength to numerically resolve.
In this paper we use an initial condition which reorganizes in a much more ordered fashion. The initial field consists of a twisted flux tube with field lines lying on nested toroidal flux surfaces with varying rotational transform. We use cylindrical coordinates
$R,z,\unicode[STIX]{x1D703}$
and flux functions
$\unicode[STIX]{x1D713}_{p}(R,z)$
and
$I(\unicode[STIX]{x1D713}_{p})$
. The coordinate system is described in figure 1. Physically
$\unicode[STIX]{x1D713}_{p}(R,z)$
represents the poloidal magnetic flux, passing through a circular surface of radius
$R$
around the
$z$
-axis (green circle).
$I$
physically represents the total poloidal current, the current through the green, shaded circle. The magnetic field is calculated from the flux functions by the standard methods for axisymmetric fields (Goedbloed, Keppens & Poedts Reference Goedbloed, Keppens and Poedts2010):
$$\begin{eqnarray}B_{R}=-\frac{1}{R}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{p}}{\unicode[STIX]{x2202}z},\quad B_{z}=\frac{1}{R}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{p}}{\unicode[STIX]{x2202}R},\quad B_{\unicode[STIX]{x1D703}}=\frac{1}{R}I.\end{eqnarray}$$
Using this construction guarantees that the magnetic field is divergence free, and that field lines lie on magnetic surfaces of constant
$\unicode[STIX]{x1D713}_{p}$
.
We choose the following flux function,
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{p}=\left\{\begin{array}{@{}ll@{}}B_{0}\cos ^{4}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}{\displaystyle \frac{a}{a_{0}}}\right),\quad & a\leqslant a_{0}\\ 0,\quad & a>a_{0},\end{array}\right.\end{eqnarray}$$
where
$0<a_{0}<1$
and with
$B_{0}$
a scaling parameter that sets the magnetic field strength and
$$\begin{eqnarray}a=\sqrt{(1-R)^{2}+z^{2}},\end{eqnarray}$$
denoting the distance from the unit circle (the circle
$R=1,z=0$
, indicated in red in figure 1), which is the magnetic axis of this initial condition.
With this choice we can see that
$\unicode[STIX]{x1D713}_{p}$
is constant and zero for
$a\geqslant a_{0}$
such that, according to (2.1), the poloidal field vanishes. The magnetic surfaces in the region where
$a<a_{0}$
form concentric tori with circular cross-section that enclose the magnetic axis.
We can choose any function of
$\unicode[STIX]{x1D713}_{p}$
for the toroidal current function
$I$
, which here we define as:
$$\begin{eqnarray}I=\frac{\unicode[STIX]{x1D713}_{p}\unicode[STIX]{x03C0}^{2}}{\imath _{0}^{\ast }a_{0}^{2}},\end{eqnarray}$$
where a scaling parameter
$\imath _{0}^{\ast }$
is introduced, named such because it sets the rotational transform on axis, as we will show in § 2.1. With this choice for
$I$
the toroidal magnetic field also vanishes for
$a\geqslant a_{0}$
, and thus the field describes an axisymmetric flux tube with major radius
$1$
and minor radius
$a_{0}$
with
$\boldsymbol{B}=0$
outside of the tube. Note: we use the convention that a subscript zero denotes a value at time
$t=0$
, and a superscript asterisk denotes that the quantity is measured on the magnetic axis.
2.1 Rotational transform profile
The winding of field lines in a toroidal magnetic structure is quantified by the rotational transform
$\imath$
or its inverse, the safety factor
$q$
. The rotational transform geometrically represents the ratio of the number of times a field line wraps around the poloidal direction of a torus to the number of times it winds around the toroidal direction. The safety factor can be calculated using the well-known formula (Wesson & Campbell Reference Wesson and Campbell2011):
$$\begin{eqnarray}q=\frac{1}{2\unicode[STIX]{x03C0}}\oint \frac{1}{R}\frac{B_{\unicode[STIX]{x1D703}}}{B_{p}}\,\text{d}l,\end{eqnarray}$$
where
$B_{p}$
is the magnitude of the poloidal magnetic field
$B_{R}\hat{R}+B_{z}\hat{z}$
, and the integration is carried out over a constant
$\unicode[STIX]{x1D703}$
cross-section of a magnetic surface
$(\unicode[STIX]{x1D713}_{p}=\text{const.})$
. This integration path is indicated by the blue circle in figure 1.
The poloidal magnetic flux enclosed between concentric magnetic surfaces of radii
$a$
and
$a+\text{d}a$
is
$\text{d}a2\unicode[STIX]{x03C0}RB_{p}$
. Hence conservation of poloidal flux implies that
$RB_{p}$
is constant on each surface. Evaluating the poloidal field at
$z=0$
,
$R\geqslant 1$
, where
$\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{p}/\unicode[STIX]{x2202}z=0$
and
$\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{p}/\unicode[STIX]{x2202}R=\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{p}/\unicode[STIX]{x2202}a$
its value in the direction of the poloidal vector
$\unicode[STIX]{x1D753}$
is equal to
$$\begin{eqnarray}RB_{p}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{p}}{\unicode[STIX]{x2202}a}=B_{0}\frac{2\unicode[STIX]{x03C0}}{a_{0}}\cos ^{3}\left(\frac{a\unicode[STIX]{x03C0}}{2a_{0}}\right)\sin \left(\frac{a\unicode[STIX]{x03C0}}{2a_{0}}\right).\end{eqnarray}$$
The toroidal magnetic field is:
$$\begin{eqnarray}B_{\unicode[STIX]{x1D703}}=B_{0}\frac{\unicode[STIX]{x03C0}^{2}\cos ^{4}\left({\displaystyle \frac{\unicode[STIX]{x03C0}a}{2a_{0}}}\right)}{R\imath _{0}^{\ast }a_{0}^{2}}.\end{eqnarray}$$
On the magnetic surfaces the parameter
$a$
is a constant, so filling this in (2.5) becomes:
$$\begin{eqnarray}q(a)=\frac{\cos \left({\displaystyle \frac{a\unicode[STIX]{x03C0}}{2a_{0}}}\right)}{4\imath _{0}^{\ast }a_{0}\sin \left({\displaystyle \frac{a\unicode[STIX]{x03C0}}{2a_{0}}}\right)}\oint \frac{1}{R}\,\text{d}l.\end{eqnarray}$$
Using
$\unicode[STIX]{x1D753}$
to parametrize the integral over the surface at
$a$
, and the identities
$R=1+a\cos (\unicode[STIX]{x1D753})$
and
$\text{d}l=a\,\text{d}\unicode[STIX]{x1D753}$
, we get:
$$\begin{eqnarray}\oint \frac{1}{R}\,\text{d}l=\mathop{\int }_{0}^{2\unicode[STIX]{x03C0}}\frac{a}{1+a\cos (\unicode[STIX]{x1D753})}\,\text{d}\unicode[STIX]{x1D753}=\frac{2\unicode[STIX]{x03C0}a}{\sqrt{1-a^{2}}}.\end{eqnarray}$$
This gives us the safety factor
$$\begin{eqnarray}q(a)=\frac{{\displaystyle \frac{a\unicode[STIX]{x03C0}}{2a_{0}}}\cot \left({\displaystyle \frac{a\unicode[STIX]{x03C0}}{2a_{0}}}\right)}{\imath _{0}^{\ast }\sqrt{1-a^{2}}},\end{eqnarray}$$
and hence a rotational transform of:
$$\begin{eqnarray}\imath (a)=\imath _{0}^{\ast }\sqrt{1-a^{2}}\frac{\tan \left({\displaystyle \frac{a\unicode[STIX]{x03C0}}{2a_{0}}}\right)}{{\displaystyle \frac{a\unicode[STIX]{x03C0}}{2a_{0}}}}.\end{eqnarray}$$
The rotational transform profile is flat near the magnetic axis, and increases to infinity when
$a\rightarrow a_{0}$
. At the magnetic axis the rotational transform is given by
$$\begin{eqnarray}\lim _{a\rightarrow 0}\imath (a)=\imath _{0}^{\ast }.\end{eqnarray}$$
The initial condition is thus an axisymmetric, twisted magnetic flux tube lying in the
$R$
,
$\unicode[STIX]{x1D703}$
-plane. We can change the twist of the magnetic field lines in the initial condition by tuning the parameter
$\imath _{0}^{\ast }$
, and the toroidal magnetic field strength with the parameter
$B_{0}$
.
3 Time evolution
We simulate the time evolution of the helical magnetic fields numerically using the resistive, viscous, compressible isothermal MHD equations. The equations are solved using the PENCIL-CODE (http://pencil-code.nordita.org/). This is a highly used solver of the MHD equations on a fixed Eulerian grid often used for astrophysical applications (Brandenburg & Dobler Reference Brandenburg and Dobler2002; Haugen, Brandenburg & Dobler Reference Haugen, Brandenburg and Dobler2004; Johansen et al. Reference Johansen, Oishi, Mac Low, Klahr, Henning and Youdin2007).
The PENCIL-CODE solves the MHD equations in terms of the vector potential
$\boldsymbol{A}$
, ensuring that the magnetic field remains divergence free. The equation of motion solved is:
$$\begin{eqnarray}\frac{\text{D}\boldsymbol{v}}{\text{D}t}=-c_{s}^{2}\unicode[STIX]{x1D735}\ln \unicode[STIX]{x1D70C}+\boldsymbol{j}\times \boldsymbol{B}/\unicode[STIX]{x1D70C}+\boldsymbol{F}_{\text{visc}}/\unicode[STIX]{x1D70C},\end{eqnarray}$$
where
$\boldsymbol{B}$
is calculated through
$\boldsymbol{B}=\unicode[STIX]{x1D735}\times \boldsymbol{A}$
and the current
$\boldsymbol{j}=\unicode[STIX]{x1D735}\times \boldsymbol{B}$
. The fluid velocity is
$\boldsymbol{v}$
and the convective derivative is denoted by
$\text{D}/\text{D}t\equiv (\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t)+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$
. The simulation is isothermal, so the pressure is related to the density by
$p=\unicode[STIX]{x1D70C}c_{s}^{2}$
where
$c_{s}^{2}$
is the sound speed (set to unity). The viscous force
$\boldsymbol{F}_{\text{visc}}$
is calculated using the rate of strain tensor
$\unicode[STIX]{x1D64E}$
whose indices are given by
$\unicode[STIX]{x1D61A}_{ij}=(1/2)(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})-(1/3)\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}.$
through the equation:
$$\begin{eqnarray}\boldsymbol{F}_{\text{visc}}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }2\unicode[STIX]{x1D708}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D64E}.\end{eqnarray}$$
The continuity equation is implemented in terms of the logarithm of the density
$\unicode[STIX]{x1D70C}$
as follows
$$\begin{eqnarray}\frac{\text{D}\ln \unicode[STIX]{x1D70C}}{\text{D}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}.\end{eqnarray}$$
Since the plasma is modelled as isothermal the equation of state does not need to be solved, so the final equation is the induction equation, which in terms of the vector potential becomes
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{A}}{\unicode[STIX]{x2202}t}=\boldsymbol{v}\times \boldsymbol{B}+\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{A},\end{eqnarray}$$
where we have chosen the Weyl gauge for
$\boldsymbol{A}$
to simplify the equation and
$\unicode[STIX]{x1D702}$
is the resistivity.
These equations are solved on an Eulerian grid of
$256^{3}$
grid points in a simulation box of size 5. This puts the boundary at
$R=2.5$
. The simulation is run using perpendicular-field boundary conditions, by imposing a vanishing parallel component of the magnetic field on the boundary. This boundary condition allows magnetic field to escape from the simulation volume. The full field information is saved every 5 simulation time steps. The simulation is initialized with a constant pressure
$p=1$
throughout the volume and the velocity field is zero. The Alfvén speed
$v_{A}$
in the initial field ranges from 0.28 to 1.1 (on axis)
$B_{0}$
is varied between 0.05 and 0.2 and this makes it of the same order of magnitude as the sound speed
$c_{s}^{2}=1$
. We scale the time to a resistive time scale using
$t_{\unicode[STIX]{x1D702}}=R_{\text{char}}^{2}/\unicode[STIX]{x1D702}$
. For the characteristic length scale
$R_{\text{char}}$
we choose the distance of the magnetic axis from the origin in the initial field,
$R^{\ast }(0)=1$
. The viscosity parameter
$\unicode[STIX]{x1D708}$
is set to
$2\times 10^{-4}$
whereas the resistivity
$\unicode[STIX]{x1D702}$
is varied from
$2\times 10^{-4}$
to almost an order of magnitude lower. This makes the magnetic Prandtl number equal to unity or larger.
In our analysis of the simulation results we wish to extract the topological properties of the magnetic field structures. We do this by means of a Runge–Kutta field line integration. From the resultant field line traces we can find the magnetic axis and the rotational transform of the field line if it lies on a toroidal surface. One implementation of this field line tracing is described in the supplemental material of Smiet et al. (Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015). This method was used in §§ 4 and 5.
Furthermore we have developed an alternative implementation in the parallel computing platform CUDA to run on graphics hardware. The hardware accelerated trilinear interpolation and massive parallelization allow for a high speed up compared to CPU-based field line tracing. This code has been made publicly available (de Jong, Kok & Smiet Reference de Jong, Kok and Smiet2018). In the
$(R,0,z)$
-plane field lines are traced from a
$1024\times 1024$
grid and for every field line the rotational transform is calculated.
The evolution of the rotational transform profile can be seen in supplemental movie 1 available at https://doi.org/10.1017/S0022377818001290. This video shows the colour-coded rotational transform, similar to figure 2(a), for all times sequentially. Additionally, this video gives a good indication of how the magnetic structure evolves in time.
Resistive evolution of the magnetic structure. Parameters are:
$\imath _{0}^{\ast }=3$
,
$B_{0}=0.05$
and
$\unicode[STIX]{x1D702}=2\times 10^{-4}$
. (a) Cross-sections in the
$R$
,
$z$
-plane at different times. The colour indicates the rotational transform of a field line starting at that position. The configuration is seen to first contract onto the
$z$
-axis and then slowly expand. The horizontal lines indicates the location where the rotational transform, shown in (b), is taken. The cross-sections shown are at times which correspond with the top six coloured lines in (b) and are exponentially spaced. (b) Evolution of the rotational transform profile and the location of the magnetic axis. The black line along the centre is the location of the magnetic axis. The coloured part around the black line indicates the time at which the magnetic axis was at that location. The top six coloured lines correspond to the six cross-sections shown in (a), and their colour again indicates the time at which that rotational transform profile was calculated. The red dashed line is the analytical result of (2.11). We can see that the configuration quickly reaches a nearly flat rotational transform profile and subsequently evolves self-similarly.
The results of this analysis are shown in figure 2(a), for the field with
$\imath _{0}^{\ast }=3$
,
$B_{0}=0.05$
and
$\unicode[STIX]{x1D702}=2\times 10^{-4}$
. Every single pixel in these images is the result of calculating the rotational transform of a field line trace. In this run the magnetic field remains axisymmetric and the magnetic axis remains in the plane defined by
$z=0$
. This is the case for all runs presented in this paper, even though these symmetries are in no way enforced by the computational procedure. With higher values of
$\imath _{0}^{\ast }$
the structure can become susceptible to a non-axisymmetric kink instability, which will be the subject of a future publication.
Figure 2(b) shows the rotational transform profile at different times during the evolution of the configuration. The rotational transform profile quickly shifts from positively curved to nearly flat with a slight negative curvature. After this initial phase the rotational transform profile remains nearly constant in space but decreases in time. In the next sections we will study this decay.
In the rest of the paper we will analyse the change in rotational transform and the location of the magnetic axis in time. We will look at the effect of changing the resistivity
$\unicode[STIX]{x1D702}$
and the initial rotational transform
$\imath _{0}^{\ast }$
.
4 Resistive growth and decay of rotational transform
In figure 2 it is seen that the magnetic axis first shifts inwards and then slowly moves back out. We follow this dynamic of the magnetic structure by extracting
$R^{\ast }$
, the distance from the origin to the magnetic axis as a function of resistive time. This is done for three different values of
$\imath _{0}^{\ast }$
at constant
$\unicode[STIX]{x1D702}=4\times 10^{-4}$
. The results are shown in figure 3. The structure relaxes to a radius which depends on
$\imath _{0}^{\ast }$
, and then slowly increases in size.
As the structure grows, it has reached its equilibrium: the magnetic pressure pushes outwards, but the expansion is halted by the strong external pressure. A lowering of the pressure is observed in a toroidal region as we show here in figure 6, and is described in more detail in Smiet et al. (Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015, Reference Smiet, Candelaresi and Bouwmeester2017).
Position of the magnetic axis in time for different values of
$\imath _{0}^{\ast }$
. The magnetic axis performs a damped oscillatory motion towards an equilibrium position which depends on the initial rotational transform, and then slowly grows. The parameters
$B_{0}=0.05$
and
$\unicode[STIX]{x1D702}=2\times 10^{-4}$
are fixed with
$\imath _{0}^{\ast }$
varied as shown.
We can understand this initial relaxation qualitatively from the interplay between magnetic tension and magnetic pressure. Since the initial pressure is constant and the velocity is zero the initial motion of the fluid is purely due to the Lorentz force
$\boldsymbol{j}\times \boldsymbol{B}$
. A high value of
$\imath _{0}^{\ast }$
results in a high poloidal field, and magnetic tension along the field lines squeezes the configuration into an expanding ring. The case of a low rotational transform will result in stronger toroidal field, and the magnetic tension will cause the structure to contract. In figure 3 we see that higher
$\imath _{0}^{\ast }$
leads to an initial expansion, and an equilibrium with a value of
$R^{\ast }$
larger than
$1$
, whereas low values of
$\imath _{0}^{\ast }$
lead to an initial contraction.
$R^{\ast }$
performs a damped Alfvénic oscillation to the equilibrium position, and then slowly grows.
The later evolution of the structure proceeds on a purely resistive time scale. This is tested by simulating the evolution of the field with the parameters
$\imath _{0}^{\ast }=3$
,
$B_{0}=0.05$
and varying resistivity
$\unicode[STIX]{x1D702}=2\times 10^{-4},1\times 10^{-4}$
and
$5\times 10^{-5}$
. When rescaled to resistive time the growth of the structure and the change in rotational transform all collapse to a single curve indicating a universal growth mechanism, as shown in figure 4.
Values of
$\imath ^{\ast }(t)$
and
$R^{\ast }(t)$
as functions of resistive time for several different values of resistivity
$\unicode[STIX]{x1D702}$
. The initial rotational transform is set to
$\imath _{0}^{\ast }=3$
, and magnetic field strength
$B_{0}=0.05$
. The change of the rotational transform and the radius of the structure all behave identically on a resistive time scale.
The magnetic field strength
$B_{0}$
does not affect the equilibrium reached or the rate of growth and change in rotational transform exhibited by these configurations. This is shown in figure 5, where the fields with
$B_{0}=0.05$
and
$B_{0}=0.2$
are compared for
$\imath _{0}^{\ast }=3$
,
$\unicode[STIX]{x1D702}=2\times 10^{-4}$
. Despite a factor 4 difference in the magnetic field strength the structures behave identically except for the initial reconfiguration towards the equilibrium. As this reconfiguration is mediated by magnetic forces, it proceeds on an Alfvénic time scale linear in the magnetic field,
$\unicode[STIX]{x1D70F}_{A}=B/\sqrt{\unicode[STIX]{x1D70C}}$
. It is therefore not surprising that the oscillation to the equilibrium
$R^{\ast }$
lasts approximately 4 times longer for the field where the magnetic amplitude is a quarter of the strength.
Magnetic decay of topologically identical structures with different initial magnetic field strength
$B_{0}$
. Despite the difference in magnetic field strength the change in rotational transform proceeds at exactly the same rate with identical equilibrium. Note that the initial oscillations towards the equilibrium radius occur on the Alfvénic time scale: the oscillations to the equilibrium configuration proceed at a four times faster rate when
$B_{0}=0.2$
then when
$B_{0}=0.05$
. In these runs
$\imath _{0}^{\ast }=3$
and
$\unicode[STIX]{x1D702}=2\times 10^{-4}$
.
5 Pfirsch–Schlüter diffusion
We can understand the structure growth and change in rotational transform through the effect of finite resistivity on the plasma and the lowering of magnetic field strength through resistivity. In a perfectly conducting plasma a magnetic field is effectively ‘frozen-in’ and moves with the fluid motion (Alfvén Reference Alfvén1943; Batchelor Reference Batchelor1950; Priest & Forbes Reference Priest and Forbes2000), thus there can be no net flow of fluid perpendicular to the field lines if the magnetic configuration is static. When resistivity is included this restriction is lifted and the fluid can slip against the static magnetic field lines. Field line slip is observed in many different scenarios and is one of the driving mechanisms behind two-dimensional reconnection (Kulsrud Reference Kulsrud2011). In the toroidal geometry of an operating tokamak, field line slip gives rise to slow diffusion out of the toroidal flux surfaces. This process is called Pfirsch–Schlüter diffusion (Wesson & Campbell Reference Wesson and Campbell2011). This Pfirsch–Schlüter flow is directed outwards, in the direction of the pressure gradient.
In the self-organized structures considered here a similar magnetic slip causes a diffusion of plasma fluid into the magnetic structure. This is shown in figure 6, where the flow field is plotted along the
$x$
-axis together with the pressure profile. The magnetic axis is located at the minimum of the pressure, and it is clearly seen how there is net fluid flow directed towards the magnetic axis. We suggest that the slight discrepancy between the location of the magnetic axis and the zero of the velocity is due to the axis itself being in motion.
Fluid velocity
$v_{\bot }$
(red), pressure
$p$
(black) and total pressure
$p+p_{M}$
along the
$x$
-axis at three different times corresponding to the second, third and fifth images in figure 2(a). The resistivity of this run was
$\unicode[STIX]{x1D702}=2\times 10^{-4}$
and
$\imath _{0}^{\ast }=3$
. During the initial reconfiguration, the pressure profile is irregular, but the total pressure is smooth. In (c), the equilibrium has reached the state where the magnetic axis is located at the minimum in pressure. The flow profile shows a net flow towards the magnetic axis which is similar to Pfirsch–Schlüter diffusion.
Whilst the fluid flow slowly penetrates the magnetic structure, the magnetic energy in the structure is decreasing. The decrease of total magnetic energy for the simulations with
$\imath _{0}^{\ast }=3$
and
$B_{0}=0.05$
is shown in figure 7.
Decay of magnetic energy and decay of total helicity versus time for the runs with
$\imath _{0}^{\ast }=3$
and
$B_{0}=0.05$
for different values of resistivity
$\unicode[STIX]{x1D702}$
.
One important result to note is that the magnetic field strength decays fast compared to the resistive decay time
$t_{\unicode[STIX]{x1D702}}$
, whereas in general the magnetic field strength is expected to evolve as
$\langle |\boldsymbol{B}|\rangle \sim \langle B_{0}\rangle \text{e}^{-t/t_{\unicode[STIX]{x1D702}}}$
. Here the magnetic energy has already decreased an order of magnitude in only
$0.1t_{\unicode[STIX]{x1D702}}$
. This is because the resistive losses are not the only mechanism through which the magnetic field is lowered: during the evolution the entire configuration also expands. Even with zero resistivity such an expansion leads to a lowering of the magnetic field strength. This can be seen as follows: since it is the flux through a co-moving surface that is conserved, and if that surface expands, the magnetic field strength lowers. This effect can also be seen in the zero resistance simulations presented in (Smiet et al.
Reference Smiet, Candelaresi and Bouwmeester2017).
Figure 7 also shows the evolution of magnetic helicity, which decreases at a slower rate than the magnetic energy. The slower decay of helicity can be seen as the result of the expansion, as the structure evolves with a self-similar shape. The magnetic energy is the integral of
$(\unicode[STIX]{x1D735}\times \boldsymbol{A})\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{A})$
, whereas the integrand of the helicity integral,
$\boldsymbol{A}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{A})$
, involves one less spatial derivation. For a similar structure of a larger size, the magnetic helicity is thus larger. Note that we evolve these structures on a time scale larger than the time scales on which the helicity can be considered conserved, and that our initial condition is intentionally very regular. Therefore the localized reconnections which transform helicity between linking, writhe and twist, which reconfigure the magnetic topology whilst leaving helicity mostly unchanged in turbulent Woltjer–Taylor-type relaxation, are absent in these runs. See Smiet et al. (Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015) for simulations where this equilibrium is achieved through this more chaotic reconnection.
Time dependence of the rotational transform for different values of
$\imath _{0}^{\ast }$
. After the initial period where the magnetic structure reorganizes on an Alfvénic time scale, the rotational transform decays following a power law with a characteristic exponent between
$-2/3$
and
$-1/2$
.
Finally, we look at the change of rotational transform in time, and how it depends on the initial value for the rotational transform. The results are shown in figure 8. From the asymptotic behaviour in the log–log plot we can see that the rotational transform decays according to a power law instead of exponentially. The characteristic exponent of this decay is different for runs with different
$\imath _{0}^{\ast }$
. The rotational transform between
$t_{\unicode[STIX]{x1D702}}=0.05$
and
$t_{\unicode[STIX]{x1D702}}=0.1$
is fitted with a power law
$\imath ^{\ast }(t)=a{t_{\unicode[STIX]{x1D702}}}^{-b}$
and a characteristic exponent of
$b=0.664$
is found for the run with
$\imath _{0}^{\ast }=10$
and
$b=0.48$
for the run with
$\imath _{0}^{\ast }=3$
. Guides are drawn in figure 8 showing
$t_{\unicode[STIX]{x1D702}}^{-2/3}$
and
$t_{\unicode[STIX]{x1D702}}^{-1/2}$
decay.
The lowering rotational transform is caused by the poloidal field decreasing faster than the toroidal field. This is another indication that the lowering of field strength is primarily caused by the expansion of the structure caused by field line slip and not necessarily by resistive decay. With resistive decay we indicate the exponential decay with characteristic decay time
$t_{\unicode[STIX]{x1D702}}$
caused by resistivity on a static field configuration, solutions of (3.4) with
$\boldsymbol{v}=0$
. Decay brought on by field line slip is caused by the non-zero Pfirsch–Schlüter flow, or conversely the motion of the field lines against the plasma. This causes the term
$\boldsymbol{v}$
to be non-zero and the term
$\unicode[STIX]{x1D6FB}^{2}\boldsymbol{A}$
to decrease (through expansion the gradients get smaller) in (3.4), and thus the change in magnetic field strength can be fast compared to the resistive decay time.
As the structure expands the poloidal field strength is lowered by expansion in the horizontal plane, as it is the poloidal flux passing through the circle defined by the magnetic axis which is conserved. The increase in area through which this flux passes is proportional to
$(R^{\ast })^{2}$
. The lowering of the toroidal field meanwhile is governed by the change in area of the poloidal cross-section of the flux tube. Expansion in this plane is constrained to the positive
$R$
-direction and as such the area should go approximately linear in
$R^{\ast }$
. This constrained expansion is attested by the D-shaped magnetic surfaces seen in figure 2(a). Even though this explanation is quite qualitative and does not take into account the shape of the surfaces and the distribution of magnetic flux, it does explain why the rotational transform lowers according to a power law. Furthermore, the difference between poloidal and toroidal expansion correctly predicts a characteristic exponent of around
$-1/2$
.
It should be noted that the decay in rotational transform is fast compared to the increase of the major radius of the structure; the rotational transform changes by a factor of three in the time
$R^{\ast }$
only changes a few per cent. This is important when considering this equilibrium as a model for magnetic clouds.
6 Relation to magnetic clouds
The presented simulations of axisymmetric equilibria are more idealized than the situation encountered in the solar wind. As noted in the introduction, the plasma
$\unicode[STIX]{x1D6FD}$
in the solar wind is approximately 1, and a signature of a magnetic cloud is that the
$\unicode[STIX]{x1D6FD}$
drops below 1 (Zurbuchen & Richardson Reference Zurbuchen and Richardson2006). In these simulations, the initial
$\unicode[STIX]{x1D6FD}$
(taking the magnetic field strength and pressure on the axis) is approximately 25 when
$B_{0}=0.05$
and 1.65 when
$B_{0}=0.2$
. However, as can be seen from figure 5, when rescaled to resistive time the evolution of the fields is near identical. Simulations (not presented in this paper) have been performed at
$\unicode[STIX]{x1D6FD}$
down to 0.7 which show the same evolution, but this is not a lower limit. The simulations show the set-up of a toroidal equilibrium against a background plasma with lower field strength, similar to the situation encountered in the solar wind.
In the solar wind the Alfvén speed and the sound speed are the same order of magnitude with
$v_{s}=37~\text{km}~\text{s}^{-1}$
and
$v_{A}=47~\text{km}~\text{s}^{-1}$
(Goedbloed & Poedts Reference Goedbloed and Poedts2004). Assuming a magnetic cloud size of approximately
$10^{6}~\text{km}$
(1 day passage for a probe travelling at
$15~\text{km}~\text{s}^{-1}$
), the Alfvénic transit time is approximately 7 h, so one dozen to several dozen Alfvén transit times pass between coronal mass ejection and observation, a similar regime as is probed in the simulation. The magnetic Prandtl number (ratio of kinematic viscosity to magnetic diffusivity) for a hot thin plasma such as the interplanetary solar wind is much higher than unity;
${\sim}10^{14}$
(see Brandenburg & Subramanian Reference Brandenburg and Subramanian2005). The plasma in the solar wind is thus in a regime here the viscous forces act faster than the resistivity to allow for a similar self-organizing process as observed in the simulations. One major difference is that
$t_{\unicode[STIX]{x1D702}}$
for a magnetic cloud is approximately
$10^{9}$
years, much longer than the months for which a cloud can be observed before it leaves the solar system. When the cloud is just ejected,
$R_{\text{char}}$
is smaller and reconnection can occur (it must to trigger the ejection). Furthermore the more chaotic solar wind allows for small-scale reconnection events which increase the effective reconnection rate. Nevertheless, the extent to which the rotational transform changes must be much smaller than the full evolution presented in this paper.
There are several models for magnetic clouds discussed in the literature, see for example Burlaga (Reference Burlaga1991) for an overview. One of the more common approaches is to model a magnetic cloud as a long flux rope extending from, and still magnetically connected to, the surface of the Sun. Nevertheless there are several models that consider magnetic clouds as localized magnetic excitations within the solar wind, with their magnetic field generated by internal currents.
Kumar and Rust describe a model for a magnetic cloud as an isolated, net current carrying toroidal flux ring (Kumar & Rust Reference Kumar and Rust1996). The magnetic field inside the ring is based on the force-free Lundquist solution valid for an infinite cylinder (Lundquist Reference Lundquist1950). As they themselves note, this cannot be an exact description, as the toroidal geometry necessitates the existence of a hoop force such as described in Garren & Chen (Reference Garren and Chen1994). In this model the plasma current is zero outside of the toroid, but the net current through the toroid is non-zero such that it generates a force-free field in the surrounding plasma.
There are several differences between Kumar and Rust’s model and a magnetic cloud as a self-organized structure we describe. Firstly the magnetic field in their model is force free. Such a field is not possible, as they note themselves because a current ring will always experience a hoop force. The rotational transform profile in their model is also very different. In Lundquist solutions the axial and tangential fields (which, when the cylinder is translated to a torus, correspond to the toroidal and poloidal directions respectively) are given by Bessel functions. As such the rotational transform profile changes significantly from the magnetic axis to the edge (Bellan Reference Bellan2000). The magnetic field outside the ring is purely toroidal, which implies that the rotational transform goes to infinity. Even though their model resembles the configuration we describe superficially, the rotational transform profile is drastically different. Our simulations show that the profile quickly flattens.
Another magnetic cloud model which resembles the configurations we observe is the flare-generated spheromak model by Ivanov & Harshiladze (Reference Ivanov and Harshiladze1985) and further explored by Vandas et al. (Reference Vandas, Fischer, Pelant and Geranios1992). They describe the clouds using the spherical force-free solution of Chandrasekhar & Kendall (Reference Chandrasekhar and Kendall1957). The magnetic topology in this solution also consists of field lines lying on nested toroidal surfaces. In this model the rotational transform profile of these force-free solutions is also non-constant.
The resistive decay time for a structure with the characteristic length scale that is reasonable for magnetic clouds,
$R_{\text{char}}=10^{6}~\text{km}$
, is approximately
$1.7\times 10^{9}$
years (Goedbloed & Poedts Reference Goedbloed and Poedts2004), so change in the rotational transform due to resistive processes is expected to be small. The process leading to the formation of the self-organized localized equilibrium however takes place on a much faster time scale. Note that not just the Alfvénic oscillation towards equilibrium is fast; as seen in figure 2 much of the change in rotational transform towards equilibrium has already occurred at
$0.006t_{\unicode[STIX]{x1D702}}$
. This change, though fast, scales with resistivity. If the resistivity is many orders of magnitude lower, as in the solar wind, this change in rotational transform can be expected to be much less rapid, and the cloud will still carry much of the topology it organized into when it was ejected. The evolution of the magnetic structure, after it is generated, can therefore be considered to be approximately ideal, as is also assumed in the model by Ivanov and Harshiladze and the model by Kumar and Rust.
Because the structure we describe does not rely on the assumption of force-free fields, an assumption that is not warranted in the
$\unicode[STIX]{x1D6FD}\sim 1$
solar wind plasma, we speculate that the magnetic structures described in this paper are a more realistic model for localized magnetic clouds than the two others described above.
7 Conclusions
In this paper we have shown how a self-organizing equilibrium evolves on a resistive time scale. In agreement with our previous studies we find that the initially twisted flux tube reconfigures on an Alfvénic time scale into an axisymmetric Grad–Shafranov equilibrium characterized by a lowered pressure on the magnetic axis.
In this paper we have described how the configuration evolves subsequently; the major radius
$R^{\ast }$
grows, and the rotational transform on the magnetic axis
$\imath ^{\ast }(t)$
lowers. The rotational transform profile, which initially had a high positive curvature, quickly evolves to a almost flat and slightly negatively curved profile.
With the exception of the initial reconfiguration which proceeds on an Alfvénic time scale, the evolution is rather independent of the resistivity when scaled to a resistive time scale. The growth of the structure can be understood as a Pfirsch–Schlüter-type slip of the field lines against the plasma fluid background, or conversely the fluid slip against the field.
It is because of this growth of the structure that the magnetic field strength decays faster than the resistive time scale. It is also this growth which allows the poloidal field to decay faster than the toroidal field. We have also given a simple geometrical argument that explains why the decay of the rotational transform behaves as a power law with characteristic exponent of the order of
$1/2$
.
In this study we have limited ourselves to isothermal (constant resistivity) MHD evolution. The inclusion of temperature would result in a spatial variation of the Spitzer resistivity, which would quantitatively change the exact evolution, but the general aspects of the equilibrium and its evolution are underlined by geometrical principles and would remain unchanged. It is an interesting question whether the generation of a flat rotational transform profile would remain robust under these conditions
These results could help make predictions for the evolution of self-organized magnetic equilibria in nature. In this paper we relate this structure to magnetic clouds. Other applications for this model include ejecta from active galactic nuclei (Braithwaite Reference Braithwaite2010), and one could possibly devise a scheme for pulsed nuclear fusion power generation in which the plasma is confined in such a magnetic structure, embedded in an extremely high fluid pressure environment. When devising such a scheme it should be important to note that the decay time of the magnetic field strength proceeds on a time scale much faster than the resistive decay time, and as such that the ‘confinement’ is a very transient phenomenon.
Acknowledgements
We wish to thank the anonymous referees for their careful review of our article. This work is part of the Rubicon programme with project number 680-50-1532, which is (partly) financed by the Netherlands Organization for Scientific Research (NWO). DIFFER is part of the Netherlands Organization for Scientific Research (NWO). Notice: this manuscript is based upon work supported by the US Department of Energy, Office of Science, Office of Fusion Energy Sciences, and has been authored by Princeton University under contract no. DE-AC02-09CH11466 with the US Department of Energy. The publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/S0022377818001290.
Open access
