**3.1 Basic assumptions of the MHD model**

The MHD model regards the plasma as a quasi-neutral electrically conducting fluid. The first and most fundamental assumption of this description is to regard the ensemble of ions and electrons conforming the plasma as a single continuum medium. This is valid when the length scales associated with the magnetic field gradients is much larger than the internal length scales of the plasma (such as the ionic and electronic gyroradii). This condition holds in virtually every laboratory plasma dedicated to fusion research.

The second important assumption is to consider that the plasma is in thermodynamic equilibrium so the particles have a Maxwellian distribution of velocities. This is a good approximation as long as the shortest time scale of the process under consideration is much longer than the collision time and the shortest length scale of the system is larger than the mean free path of the particles. In other words, the plasma should be in a collisional regime (this condition is required to derive the fluid equations from the kinetic equations (Braginskii, 1965)). The collisionality hypothesis is usually not satisfied at the highest temperatures obtained in modern tokamak experiments. However, spheromak plasmas are much colder (*<sup>T</sup>* <sup>∼</sup> 102 eV) so that this assumption is still reasonable. Moreover, there are several arguments supporting the validity of the MHD model even in collisionless systems (Friedberg, 1987); (Priest & Forbes, 2000).

Finally, in the context of the MHD model the plasma is assumed to be electrically neutral (or quasi-neutral since the charges are present but exactly balanced). This is approximately true when the length scales under consideration are larger than the Debye shielding of electrons.

### **3.2 MHD equations**

Now we seek for the equations that describe the evolution of the two main quantities that govern the dynamics of such an MHD system: the velocity field and the magnetic field. The equation for the evolution of the plasma velocity **u**, expresses the balance of linear momentum

$$
\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \mathbf{J} \times \mathbf{B} + \mu \nabla \cdot \Pi \tag{7}
$$

where *ρ* is the mass density and *p* is the thermodynamic pressure. The second term on the right hand side is the Lorentz force, where **J** is the current density and **B** is the magnetic field. We note that due to quasi-neutrality the current density is produced by the relative motion

**3.3 Diffusion of magnetic field lines and frozen-in-flux condition**

*∂B*

may be found elsewhere (Biskamp, 2000);(Priest & Forbes, 2000).

∇*p* from Eq. (7) and assume that *ρ* = *ρ*<sup>0</sup> is a constant.

in the resistive time scale *τ<sup>r</sup>* = *L*2/*η*.

the induction equation reduces to,

**3.4 Closing the system of equations**

*β* = 2*p*0/*B*<sup>2</sup>

low *β* plasmas.

The quotient between the magnitudes of these effects can be estimated as

<sup>|</sup>*η*∇2**B**<sup>|</sup> <sup>∼</sup> *<sup>u</sup>*0/*<sup>L</sup>*


The two terms in the right hand side of Eq. (12) account for two very different physical effects.

Dynamics of Magnetic Relaxation in Spheromaks 93

where *u*<sup>0</sup> and *L* are typical velocity and length scales and *Rm* is the magnetic Reynolds number. Thus, when *Rm* ∼ 0 the magnetic field simply diffuses and the configuration decays

The opposite limit (*Rm* � 1) is more representative of the actual situation in most laboratory (in the context of magnetic confinement) and space plasmas. In this limit (called ideal limit)

This equation implies the conservation of the magnetic flux through any closed surface that moves with the local velocity of the fluid. If we regard the magnetic field lines as very thin flux tubes and we imagine closed curves surrounding them that move with the fluid, we realize that the plasma drags the field lines as it moves. It is said that the field lines are frozen in the plasma (*frozen-in-flux condition*). Since each field line is simply convected by the flow (assumed to be smooth and continuous) its connectivity is preserved. This means that in the ideal MHD approximation the changes in the topology of the magnetic field are not possible. This idea, which is intimately related to the Kelvin's circulation theorem for inviscid flows, was first introduced by Hannes Alfvén in 1943. More details on the frozen in flux condition

The system formed by Eqs. (7) and (12) and the constraint ∇ · **B** = 0, has too many unknowns and can not be solved. Even if the current density can easily be expressed in terms of the magnetic field using Ampère's law, we still need to introduce some information concerning the density and the pressure. We describe four common approaches to accomplish this.

Firstly, the zero-*β* approximation gives the simplest option. The nondimensional parameter

pressure. A very low *β* value (which is the case in most confinement experiments) means that the dynamics of the plasma is mainly dictated by the magnetic field (via the Lorentz force) while the pressure gradient has little influence. Thus, we may simply remove the term

The second option is to consider an incompressible flow. In this case *ρ* is still a constant but the pressure gradient is no dropped. In this case the required information is completed by the incompressible condition ∇ · **u** = 0. The pressure can not be directly computed with this equation but it can be indirectly inferred. It plays the role of a Lagrange multiplier. Although this is a less crude and more consistent option, it is not generally a good approximation for

<sup>0</sup> measures the ratio between the thermodynamic pressure and the magnetic

*<sup>η</sup>*/*L*<sup>2</sup> <sup>=</sup> *<sup>u</sup>*0*<sup>L</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> ∇ × (**<sup>u</sup>** <sup>×</sup> **<sup>B</sup>**) (*Rm* � <sup>1</sup>). (14)

*<sup>η</sup>* <sup>≡</sup> *Rm* (13)

between ions and electrons. The last term in Eq. (7) is the viscous force where *μ* = *ρν* is the dynamic viscosity, *ν* is the cinematic viscosity and the tensor Π is given by

$$
\Pi = \nabla \mathbf{u} + \nabla \mathbf{u}^T - \frac{2}{3} (\nabla \cdot \mathbf{u}) \mathbf{I}. \tag{8}
$$

If the flow is incompressible (∇ · **u** = 0) Π reduces to ∇**u**.

Let us mention some basic aspects of the Lorentz force term. Using the low-frequency Ampère law **J** = ∇ × **B** (we rescale **B** and **J** in such a way that *μ*<sup>0</sup> = 1) and the vector identity (∇ × **<sup>B</sup>**) <sup>×</sup> **<sup>B</sup>** = (**<sup>B</sup>** · ∇)**<sup>B</sup>** − ∇(*B*2/2), we can decompose this term into two contributions

$$\mathbf{J} \times \mathbf{B} = (\mathbf{B} \cdot \nabla)\mathbf{B} - \nabla\left(\frac{B^2}{2}\right). \tag{9}$$

The first term on the right represents a magnetic tension force in the direction of **B** which has a restoring effect when the magnetic field lines are bent. The second term is regarded as a magnetic pressure that acts in all directions. Clearly, both forces must cancel out along the magnetic field lines since the term **J** × **B** can not accelerate the fluid in the direction of **B**.

The equation for the magnetic field evolution comes from the Maxwell equations and a constitutive law that relates the electric field to the magnetic field and the current density (the Ohm's law). We begin with the Faraday's law in the low-frequency limit (i.e. neglecting the displacement current)

$$
\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}.\tag{10}
$$

On the other hand, the Ohm's law relates the current density to the electric field in the frame of reference of the conducting medium **E**� = *η***J**� , where *η* is the electric resistivity (the reciprocal of the conductivity) and the prime denotes that the quantities have to be measured in the plasma's reference frame. When this equation is expressed in the lab's frame (from which the plasma moves at velocity **u**) it adopts the form

$$\mathbf{E} = -\mathbf{u} \times \mathbf{B} + \eta \mathbf{J} \tag{11}$$

where relativistic effects have been neglected (*u* � *c*, where *c* is the speed of light).

Combining Eqs. (10) and (11) together with the identity ∇×∇× **<sup>B</sup>** <sup>=</sup> <sup>∇</sup>(∇ · **<sup>B</sup>**) − ∇2**<sup>B</sup>** and the constraint ∇ · **B** = 0, we obtain the MHD induction equation

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \eta \nabla^2 \mathbf{B} \tag{12}$$

where spatial uniformity of *η* was assumed. Although not considered in this work, we point out that, whenever present, resistivity gradients may give rise to the so-called current interchange effect which constitutes an effective mechanism of current exchange between flux tubes (Zheng & Furukawa, 2010). Note that the terms **J** × **B** and **u** × **B** introduce a strong non-linear coupling between Eqs. (7) and (12).

8 Will-be-set-by-IN-TECH

between ions and electrons. The last term in Eq. (7) is the viscous force where *μ* = *ρν* is the

Let us mention some basic aspects of the Lorentz force term. Using the low-frequency Ampère law **J** = ∇ × **B** (we rescale **B** and **J** in such a way that *μ*<sup>0</sup> = 1) and the vector identity (∇ ×

The first term on the right represents a magnetic tension force in the direction of **B** which has a restoring effect when the magnetic field lines are bent. The second term is regarded as a magnetic pressure that acts in all directions. Clearly, both forces must cancel out along the magnetic field lines since the term **J** × **B** can not accelerate the fluid in the direction of **B**.

The equation for the magnetic field evolution comes from the Maxwell equations and a constitutive law that relates the electric field to the magnetic field and the current density (the Ohm's law). We begin with the Faraday's law in the low-frequency limit (i.e. neglecting

∇ × **<sup>E</sup>** <sup>=</sup> <sup>−</sup>*∂***<sup>B</sup>**

On the other hand, the Ohm's law relates the current density to the electric field in the frame of

of the conductivity) and the prime denotes that the quantities have to be measured in the plasma's reference frame. When this equation is expressed in the lab's frame (from which the

Combining Eqs. (10) and (11) together with the identity ∇×∇× **<sup>B</sup>** <sup>=</sup> <sup>∇</sup>(∇ · **<sup>B</sup>**) − ∇2**<sup>B</sup>** and

where spatial uniformity of *η* was assumed. Although not considered in this work, we point out that, whenever present, resistivity gradients may give rise to the so-called current interchange effect which constitutes an effective mechanism of current exchange between flux tubes (Zheng & Furukawa, 2010). Note that the terms **J** × **B** and **u** × **B** introduce a strong

where relativistic effects have been neglected (*u* � *c*, where *c* is the speed of light).

the constraint ∇ · **B** = 0, we obtain the MHD induction equation

*∂***B**

*∂t*

3

 *B*<sup>2</sup> 2 

(∇ · **u**)**I**. (8)

. (9)

. (10)

, where *η* is the electric resistivity (the reciprocal

**E** = −**u** × **B** + *η***J** (11)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> ∇ × (**<sup>u</sup>** <sup>×</sup> **<sup>B</sup>**) + *<sup>η</sup>*∇2**<sup>B</sup>** (12)

<sup>Π</sup> <sup>=</sup> <sup>∇</sup>**<sup>u</sup>** <sup>+</sup> <sup>∇</sup>**u***<sup>T</sup>* <sup>−</sup> <sup>2</sup>

**<sup>B</sup>**) <sup>×</sup> **<sup>B</sup>** = (**<sup>B</sup>** · ∇)**<sup>B</sup>** − ∇(*B*2/2), we can decompose this term into two contributions

**J** × **B** = (**B** · ∇)**B** − ∇

dynamic viscosity, *ν* is the cinematic viscosity and the tensor Π is given by

If the flow is incompressible (∇ · **u** = 0) Π reduces to ∇**u**.

the displacement current)

reference of the conducting medium **E**� = *η***J**�

plasma moves at velocity **u**) it adopts the form

non-linear coupling between Eqs. (7) and (12).

#### **3.3 Diffusion of magnetic field lines and frozen-in-flux condition**

The two terms in the right hand side of Eq. (12) account for two very different physical effects. The quotient between the magnitudes of these effects can be estimated as

$$\frac{|\nabla \times (\mathbf{u} \times \mathbf{B})|}{|\eta \nabla^2 \mathbf{B}|} \sim \frac{u\_0/L}{\eta/L^2} = \frac{u\_0L}{\eta} \equiv R\_m \tag{13}$$

where *u*<sup>0</sup> and *L* are typical velocity and length scales and *Rm* is the magnetic Reynolds number. Thus, when *Rm* ∼ 0 the magnetic field simply diffuses and the configuration decays in the resistive time scale *τ<sup>r</sup>* = *L*2/*η*.

The opposite limit (*Rm* � 1) is more representative of the actual situation in most laboratory (in the context of magnetic confinement) and space plasmas. In this limit (called ideal limit) the induction equation reduces to,

$$\frac{\partial B}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) \qquad (\mathcal{R}\_m \gg 1). \tag{14}$$

This equation implies the conservation of the magnetic flux through any closed surface that moves with the local velocity of the fluid. If we regard the magnetic field lines as very thin flux tubes and we imagine closed curves surrounding them that move with the fluid, we realize that the plasma drags the field lines as it moves. It is said that the field lines are frozen in the plasma (*frozen-in-flux condition*). Since each field line is simply convected by the flow (assumed to be smooth and continuous) its connectivity is preserved. This means that in the ideal MHD approximation the changes in the topology of the magnetic field are not possible. This idea, which is intimately related to the Kelvin's circulation theorem for inviscid flows, was first introduced by Hannes Alfvén in 1943. More details on the frozen in flux condition may be found elsewhere (Biskamp, 2000);(Priest & Forbes, 2000).

#### **3.4 Closing the system of equations**

The system formed by Eqs. (7) and (12) and the constraint ∇ · **B** = 0, has too many unknowns and can not be solved. Even if the current density can easily be expressed in terms of the magnetic field using Ampère's law, we still need to introduce some information concerning the density and the pressure. We describe four common approaches to accomplish this.

Firstly, the zero-*β* approximation gives the simplest option. The nondimensional parameter *β* = 2*p*0/*B*<sup>2</sup> <sup>0</sup> measures the ratio between the thermodynamic pressure and the magnetic pressure. A very low *β* value (which is the case in most confinement experiments) means that the dynamics of the plasma is mainly dictated by the magnetic field (via the Lorentz force) while the pressure gradient has little influence. Thus, we may simply remove the term ∇*p* from Eq. (7) and assume that *ρ* = *ρ*<sup>0</sup> is a constant.

The second option is to consider an incompressible flow. In this case *ρ* is still a constant but the pressure gradient is no dropped. In this case the required information is completed by the incompressible condition ∇ · **u** = 0. The pressure can not be directly computed with this equation but it can be indirectly inferred. It plays the role of a Lagrange multiplier. Although this is a less crude and more consistent option, it is not generally a good approximation for low *β* plasmas.

where *τ<sup>r</sup>* is the time scale of resistive dissipation. In the simulations presented in Sec. 5 we

Dynamics of Magnetic Relaxation in Spheromaks 95

It is common to observe magnetized fluids and plasmas as well as other continuum media to exist naturally in states with some kind of large scale order. These states are to some extend independent of the initial conditions, that is to say, they are preferred configurations toward which the system evolves if the correct boundary conditions are imposed. Moreover, if the system is perturbed it tends to return to the same preferred state recovering the large scale order. The large scale order of some quantity always comes together with the disorder at small scales of another quantity. These preferred configurations are called *self-organized* sates and the dynamical process of achieving these states is called *self-organization* (Hasegawa, 1985).

Self-organized (or relaxed) states can not be deduced from force balance or stability considerations alone. The theory of magnetic relaxation always relies on some variational principle, that is to say, the minimization of some quantity subjected to one or more constraints. Possibly the simplest and surely the most widespread option adopted to describe plasma relaxation was introduced by Taylor (1974). While a rather obvious choice was made for the quantity to minimize (the magnetic energy) a very clever option was made for the constraint. Among all the ideal MHD invariants Taylor chose the total magnetic helicity. The total magnetic helicity quantifies several topological properties of the system and even when magnetic reconnection can change the topology of the magnetic field lines, the total helicity of

the system is still a robust invariant. These ideas are further developed below.

The total magnetic helicity *H* of the magnetic field **B** within the volume *V* is

*H* = *V*

where **A** is the potential vector (**B** = ∇ × **A**). A relevant question may be posed at this time regarding how this quantity is modified by a change in the gauge of **A**. It is clear that in order to have a meaningful definition, Eq. (20) should be gauge invariant. The helicity change Δ*H*

> *V*

∇*χ* · **B** *dV* =

Δ*H* = *S*

where we have used the fact that ∇ · **B** = 0. Applying the divergence theorem in a simply

where *S* is the surface that encloses *V* and **ds** is the outward-pointing normal surface element. Therefore, the definition (20) is gauge invariant only if the normal component of the magnetic field vanishes at the boundary of *V*, which was assumed to be simply connected. We will

**A** · **B** *dV* (20)

∇ · (*χ***B**) *dV* (21)

*χ***B** · **ds** (22)

have used *Pm* <sup>=</sup> 1 and *<sup>S</sup>* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>4.

Plasma relaxation is an example of self-organization.

introduced when **A** is replaced by **A** + ∇*χ* is

connected volume *V* this becomes

Δ*H* = *V*

**4. Plasma relaxation**

**4.1 Magnetic helicity**

Thirdly, we could allow compressible flows by using the continuity equation

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \tag{15}$$

and assuming some simple relationship between the pressure and the density. For instance, *p* = *c*<sup>2</sup> *<sup>s</sup> ρ* where *cs* is the speed of sound (isothermal approximation) or *p* = *cadρ<sup>γ</sup>* where *cad* is a constant and *γ* is the polytropic index (polytropic approximation).

Finally, a more elaborated option can be obtained if we incorporate, besides Eq. (15), an equation for the energy balance

$$\frac{\partial w}{\partial t} = -\nabla \cdot \mathbf{q} - Q\_{\mathcal{E}\mathcal{I}} \tag{16}$$

where *w* is the total energy density defined as

$$w = \rho \left(\frac{\mu^2}{2} + e\right) + \frac{B^2}{2} \tag{17}$$

*e* is the internal energy per unit of mass, the term *Qc*,*<sup>r</sup>* accounts for conductive and radiative losses and **q** is the energy flux given by

$$\mathbf{q} = \left[\rho\left(\boldsymbol{\varrho} + \frac{\boldsymbol{\mu}^2}{2}\right) + p\right] \mathbf{u} + \mathbf{E} \times \mathbf{B} - \mu\left(\mathbf{u} \cdot \boldsymbol{\Pi}\right). \tag{18}$$

The three terms on the right denote respectively the energy flux due to convection, the electromagnetic energy flux (Poynting's vector) and the viscous dissipation of energy. In this context we also need an equation of state of the form *p* = *p*(*ρ*,*e*). The most common choice is the ideal gas law *p* = (*γ* − 1)*ρe* where *γ* is the ratio of specific heats.

#### **3.5 Scales and dimensionless numbers**

The results presented in Sec. 5 are obtained by numerically solving the MHD equations. It is a common practice to use a nondimensionalized version of the considered equations. The removal of the units is achieved by the choice of suitable scales that can be condensed in few nondimensional numbers. We will see which are the chosen scales and the relevant nondimensional quantities in this study.

Since spheromaks are very low-*β* plasmas, the zero-*β* approximation is used to close the system. The resulting equations can be nondimensionalized with a length scale *a* (the radius of the cylinder inside which we will solve the equations) and a velocity scale *cA* <sup>=</sup> *<sup>B</sup>*/√*ρ*, which is known as the Alfvén velocity. Perpendicular perturbations travel along the magnetic field lines at this velocity. The time will be normalized by the Alfvén time *τ<sup>A</sup>* = *a*/*cA*, which represents the typical time scale of the MHD fluctuations.

Using these scales we obtain two nondimensional numbers: the Lundquist number *S* = *acA*/*η* and *ν*/*acA* which is usually expressed in terms of the magnetic Prandtl *Pm* = *ν*/*η*. The Lundquist number can be rewritten as

$$S = \frac{a^2/\eta}{a/c\_A} = \frac{\tau\_r}{\tau\_A} \tag{19}$$

where *τ<sup>r</sup>* is the time scale of resistive dissipation. In the simulations presented in Sec. 5 we have used *Pm* <sup>=</sup> 1 and *<sup>S</sup>* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>4.

## **4. Plasma relaxation**

10 Will-be-set-by-IN-TECH

and assuming some simple relationship between the pressure and the density. For instance,

Finally, a more elaborated option can be obtained if we incorporate, besides Eq. (15), an

*<sup>s</sup> ρ* where *cs* is the speed of sound (isothermal approximation) or *p* = *cadρ<sup>γ</sup>* where *cad* is

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ · (*ρ***u**) = <sup>0</sup> (15)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> −∇ · **<sup>q</sup>** <sup>−</sup> *Qc*,*<sup>r</sup>* (16)

<sup>2</sup> (17)

. (18)

(19)

Thirdly, we could allow compressible flows by using the continuity equation

*∂ρ*

*∂w*

*w* = *ρ*

 *u*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>e</sup>* + *B*2

*e* is the internal energy per unit of mass, the term *Qc*,*<sup>r</sup>* accounts for conductive and radiative

The three terms on the right denote respectively the energy flux due to convection, the electromagnetic energy flux (Poynting's vector) and the viscous dissipation of energy. In this context we also need an equation of state of the form *p* = *p*(*ρ*,*e*). The most common choice is

The results presented in Sec. 5 are obtained by numerically solving the MHD equations. It is a common practice to use a nondimensionalized version of the considered equations. The removal of the units is achieved by the choice of suitable scales that can be condensed in few nondimensional numbers. We will see which are the chosen scales and the relevant

Since spheromaks are very low-*β* plasmas, the zero-*β* approximation is used to close the system. The resulting equations can be nondimensionalized with a length scale *a* (the radius of the cylinder inside which we will solve the equations) and a velocity scale *cA* <sup>=</sup> *<sup>B</sup>*/√*ρ*, which is known as the Alfvén velocity. Perpendicular perturbations travel along the magnetic field lines at this velocity. The time will be normalized by the Alfvén time *τ<sup>A</sup>* = *a*/*cA*, which

Using these scales we obtain two nondimensional numbers: the Lundquist number *S* = *acA*/*η* and *ν*/*acA* which is usually expressed in terms of the magnetic Prandtl *Pm* = *ν*/*η*.

> <sup>=</sup> *<sup>τ</sup><sup>r</sup> τA*

*<sup>S</sup>* <sup>=</sup> *<sup>a</sup>*2/*<sup>η</sup> a*/*cA*

**u** + **E** × **B** − *μ*

 **<sup>u</sup>** · <sup>Π</sup>

a constant and *γ* is the polytropic index (polytropic approximation).

*p* = *c*<sup>2</sup>

equation for the energy balance

where *w* is the total energy density defined as

losses and **q** is the energy flux given by

**3.5 Scales and dimensionless numbers**

nondimensional quantities in this study.

The Lundquist number can be rewritten as

represents the typical time scale of the MHD fluctuations.

**q** = *ρ e* + *u*2 2 + *p* 

the ideal gas law *p* = (*γ* − 1)*ρe* where *γ* is the ratio of specific heats.

It is common to observe magnetized fluids and plasmas as well as other continuum media to exist naturally in states with some kind of large scale order. These states are to some extend independent of the initial conditions, that is to say, they are preferred configurations toward which the system evolves if the correct boundary conditions are imposed. Moreover, if the system is perturbed it tends to return to the same preferred state recovering the large scale order. The large scale order of some quantity always comes together with the disorder at small scales of another quantity. These preferred configurations are called *self-organized* sates and the dynamical process of achieving these states is called *self-organization* (Hasegawa, 1985). Plasma relaxation is an example of self-organization.

Self-organized (or relaxed) states can not be deduced from force balance or stability considerations alone. The theory of magnetic relaxation always relies on some variational principle, that is to say, the minimization of some quantity subjected to one or more constraints. Possibly the simplest and surely the most widespread option adopted to describe plasma relaxation was introduced by Taylor (1974). While a rather obvious choice was made for the quantity to minimize (the magnetic energy) a very clever option was made for the constraint. Among all the ideal MHD invariants Taylor chose the total magnetic helicity. The total magnetic helicity quantifies several topological properties of the system and even when magnetic reconnection can change the topology of the magnetic field lines, the total helicity of the system is still a robust invariant. These ideas are further developed below.

#### **4.1 Magnetic helicity**

The total magnetic helicity *H* of the magnetic field **B** within the volume *V* is

$$H = \int\_{V} \mathbf{A} \cdot \mathbf{B} \,dV\tag{20}$$

where **A** is the potential vector (**B** = ∇ × **A**). A relevant question may be posed at this time regarding how this quantity is modified by a change in the gauge of **A**. It is clear that in order to have a meaningful definition, Eq. (20) should be gauge invariant. The helicity change Δ*H* introduced when **A** is replaced by **A** + ∇*χ* is

$$
\Delta H = \int\_{V} \nabla \chi \cdot \mathbf{B} \,dV = \int\_{V} \nabla \cdot (\chi \mathbf{B}) \,dV \tag{21}
$$

where we have used the fact that ∇ · **B** = 0. Applying the divergence theorem in a simply connected volume *V* this becomes

$$
\Delta H = \int\_{S} \chi \mathbf{B} \cdot \mathbf{ds} \tag{22}
$$

where *S* is the surface that encloses *V* and **ds** is the outward-pointing normal surface element. Therefore, the definition (20) is gauge invariant only if the normal component of the magnetic field vanishes at the boundary of *V*, which was assumed to be simply connected. We will

for the helicity due to the linking of the tubes. In general, if each tube winds around the other *N* times (in Fig. 2 (b) we have *N* = 6) one obtains *H*<sup>1</sup> = *H*<sup>2</sup> = *N*Φ1Φ<sup>2</sup> (Moffat, 1978). *N* is the

Dynamics of Magnetic Relaxation in Spheromaks 97

The helicity can also measure the self-twisting of a single tube. So far we did not pay attention to the tube's cross section shape. In order to compute the helicity of a twisted flux tube it is convenient to consider structures having elongated cross sections, like the ribbon shown in Fig. 3 (a). This ribbon is untwisted and has no helicity. If we cut this ribbon, we rotate one end

Fig. 3. (a) An untwisted ribbon-like flux tube, (b) a twisted flux tube and (c) the same twisted

by 2*π* and we join both ends together again we obtain the twisted ribbon shown in Fig. 3 (b). The helicity of this structure may be computed using Eq. (27) obtained for two linked tubes applying the following reasoning. Regard the twisted tube as two adjacent tubes carrying one half of the total magnetic flux (see Fig. 3 (c)). The helicity of this system has a contribution

twisting of the smaller tubes. If Φ is the total magnetic flux of the original twisted tube, the

Note that this mental process to convert helicity due to twisting into helicity due to linking

 Φ 2*n* 2

for each contribution due to linking. Finally, the helicity of the twisted tube is obtained by

link <sup>=</sup> <sup>Φ</sup><sup>2</sup> <sup>∞</sup>

<sup>=</sup> <sup>Φ</sup><sup>2</sup>

∑ *n*=1 1

 Φ 2 2

*H*1 link = 2

*H<sup>n</sup>* link <sup>=</sup> <sup>2</sup>*<sup>n</sup>*

> *N* ∑ *n*=1 *H<sup>n</sup>*

*H* = lim *N*→∞ link and also a contribution coming from the self

. (28)

<sup>2</sup>*<sup>n</sup>* (29)

<sup>2</sup>*<sup>n</sup>* <sup>=</sup> <sup>Φ</sup>2. (30)

*linking number* of the tubes.


coming from the linking of the two tubes *H*<sup>1</sup>

helicity due to the linking of each half is

can be recursively applied to obtain

adding all these contributions

tube marked with different colors.

respect these two conditions throughout this work. When the normal component of the magnetic field does not vanish at the boundary or the volume *V* is not simply connected a generalized definition, the so-called *relative helicity*, must be employed (Finn & Antonsen, 1985).

To see how *H* can measure topological properties of the system we will consider the concept of flux tube. The magnetic flux through an open and orientable surface *S* is

$$
\Phi = \int\_{S} \mathbf{B} \cdot \mathbf{ds} = \oint\_{\mathbb{C}} \mathbf{A} \cdot \mathbf{dl} \tag{23}
$$

where *C* is the path along the perimeter of *S* in the counterclockwise direction. Note that the Eq. (23) holds even if the gauge of **A** is changed.

We present an example given by Moffat (1978). Consider two linked flux tubes like those shown in Fig. 2 (a). We assume that there are no other contributions to the magnetic field. In

Fig. 2. Linked flux tubes.

this simplified case the total helicity can be computed as *<sup>H</sup>* = *<sup>H</sup>*<sup>1</sup> + *<sup>H</sup>*2, with *Hi* = *Vi* **A** · **B** *dV*, for *i* = 1, 2. To compute *Hi* we note that *dV* = **dl** · **ds**, where **dl** is the element of length along the tube and **ds** its cross section. By construction, **dl** and **ds** are parallel to **B**, so we can rearrange the integrand as **A** · **B** *dV* = **A** · **B dl** · **ds** = (**A** · **dl**)(**B** · **ds**), and thus

$$H\_{\bar{l}} = \oint\_{\mathcal{C}\_{l}} \int\_{S\_{l}} (\mathbf{A} \cdot \mathbf{dl})(\mathbf{B} \cdot \mathbf{ds}).\tag{24}$$

Since the magnetic flux is constant along each curve *Ci* the last equation can be written as

$$H\_i = \Phi\_i \oint\_{\mathbf{C}\_i} \mathbf{A} \cdot \mathbf{dl}.\tag{25}$$

On the other hand, the contour *C*<sup>1</sup> encloses the magnetic flux Φ<sup>2</sup> and vice versa, so from Eq. (23) it is clear that

$$
\oint\_{\mathbb{C}\_1} \mathbf{A} \cdot \mathbf{dl} = \Phi\_2 \qquad \text{and} \qquad \oint\_{\mathbb{C}\_2} \mathbf{A} \cdot \mathbf{dl} = \Phi\_1 \tag{26}
$$

and thus *H*<sup>1</sup> = *H*<sup>2</sup> = Φ1Φ<sup>2</sup> which finally gives

$$H\_{\rm link} = 2\Phi\_1 \Phi\_2 \tag{27}$$

12 Will-be-set-by-IN-TECH

respect these two conditions throughout this work. When the normal component of the magnetic field does not vanish at the boundary or the volume *V* is not simply connected a generalized definition, the so-called *relative helicity*, must be employed (Finn & Antonsen,

To see how *H* can measure topological properties of the system we will consider the concept

**B** · **ds** =

where *C* is the path along the perimeter of *S* in the counterclockwise direction. Note that the

We present an example given by Moffat (1978). Consider two linked flux tubes like those shown in Fig. 2 (a). We assume that there are no other contributions to the magnetic field. In

this simplified case the total helicity can be computed as *<sup>H</sup>* = *<sup>H</sup>*<sup>1</sup> + *<sup>H</sup>*2, with *Hi* =

rearrange the integrand as **A** · **B** *dV* = **A** · **B dl** · **ds** = (**A** · **dl**)(**B** · **ds**), and thus

*Hi* = *Ci Si*

 *C*1

and thus *H*<sup>1</sup> = *H*<sup>2</sup> = Φ1Φ<sup>2</sup> which finally gives

(a) (b)

for *i* = 1, 2. To compute *Hi* we note that *dV* = **dl** · **ds**, where **dl** is the element of length along the tube and **ds** its cross section. By construction, **dl** and **ds** are parallel to **B**, so we can

Since the magnetic flux is constant along each curve *Ci* the last equation can be written as

 *Ci*

On the other hand, the contour *C*<sup>1</sup> encloses the magnetic flux Φ<sup>2</sup> and vice versa, so from Eq.

 *C*2

*Hi* = Φ*<sup>i</sup>*

**A** · **dl** = Φ<sup>2</sup> and

 *C*

**A** · **dl** (23)

*Vi*

(**A** · **dl**)(**B** · **ds**). (24)

**A** · **dl**. (25)

*H*link = 2Φ1Φ<sup>2</sup> (27)

**A** · **dl** = Φ<sup>1</sup> (26)

**A** · **B** *dV*,

of flux tube. The magnetic flux through an open and orientable surface *S* is

Φ = *S*

Eq. (23) holds even if the gauge of **A** is changed.

**-**

Fig. 2. Linked flux tubes.

(23) it is clear that

1985).

for the helicity due to the linking of the tubes. In general, if each tube winds around the other *N* times (in Fig. 2 (b) we have *N* = 6) one obtains *H*<sup>1</sup> = *H*<sup>2</sup> = *N*Φ1Φ<sup>2</sup> (Moffat, 1978). *N* is the *linking number* of the tubes.

The helicity can also measure the self-twisting of a single tube. So far we did not pay attention to the tube's cross section shape. In order to compute the helicity of a twisted flux tube it is convenient to consider structures having elongated cross sections, like the ribbon shown in Fig. 3 (a). This ribbon is untwisted and has no helicity. If we cut this ribbon, we rotate one end

Fig. 3. (a) An untwisted ribbon-like flux tube, (b) a twisted flux tube and (c) the same twisted tube marked with different colors.

by 2*π* and we join both ends together again we obtain the twisted ribbon shown in Fig. 3 (b). The helicity of this structure may be computed using Eq. (27) obtained for two linked tubes applying the following reasoning. Regard the twisted tube as two adjacent tubes carrying one half of the total magnetic flux (see Fig. 3 (c)). The helicity of this system has a contribution coming from the linking of the two tubes *H*<sup>1</sup> link and also a contribution coming from the self twisting of the smaller tubes. If Φ is the total magnetic flux of the original twisted tube, the helicity due to the linking of each half is

$$H\_{\rm link}^1 = 2\left(\frac{\Phi}{2}\right)^2. \tag{28}$$

Note that this mental process to convert helicity due to twisting into helicity due to linking can be recursively applied to obtain

$$H\_{\rm link}^{\rm n} = 2^{\rm n} \left(\frac{\Phi}{2^n}\right)^2 = \frac{\Phi^2}{2^n} \tag{29}$$

for each contribution due to linking. Finally, the helicity of the twisted tube is obtained by adding all these contributions

$$H = \lim\_{N \to \infty} \sum\_{n=1}^{N} H\_{\text{link}}^{\eta} = \Phi^2 \sum\_{n=1}^{\infty} \frac{1}{2^n} = \Phi^2. \tag{30}$$

*W* =

*<sup>V</sup> <sup>B</sup>*2/2 *dV*, decays at the rate

*∂***A**

*dH dt* <sup>=</sup> *V*

Ohm's law **E** = *η***J**, we finally obtain

fluctuations, such as the tokamak.

injection.

*dW dt* <sup>=</sup> <sup>−</sup>*<sup>η</sup>*

*∂***B** *∂t dV* = *V*

 *V ∂***A**

*dH*

*dt* <sup>=</sup> <sup>−</sup>2*<sup>η</sup>*

current densities proportional to *Bη*<sup>−</sup>1/2. In such case, the energy decay rate becomes

*dt* <sup>∝</sup> <sup>2</sup>*η*1/2

so that as *η* tends to zero the helicity dissipation becomes negligible. Therefore, it is important to keep in mind that a plasma will relax (in the sense described here) only if there is a certain level of small scale fluctuations that gives rise to many localized current sheets. For this reason, relaxation theory does not apply to devices with a very low level of magnetic

Localized magnetic reconnection events may redistribute currents in the plasma by helicity transfer between flux tubes. Even when this idea must be applied with care because the helicity is by definition a global quantity, it is clear that the helicity of a single flux tube may change after reconnection with another flux tube. This helicity transfer process is certainly at work during toroidal current drive in spheromaks and other devices sustained by helicity

*dW dt* <sup>∝</sup> *V*

which is independent of *η*. On the other hand, the total helicity decays as

*dH*

the magnetic helicity we take the time derivative of (20) and obtain

*<sup>∂</sup><sup>t</sup>* · **<sup>B</sup>** <sup>+</sup> **<sup>A</sup>** ·

*dH dt* <sup>=</sup> <sup>2</sup>  *V J*

Dynamics of Magnetic Relaxation in Spheromaks 99

This expression can be obtained by scalar multiplying Eq. (10) by **B** and integrating over a fixed volume at whose boundary the normal component of the Poynting vector vanishes. For

Using vector identities and the divergence theorem, the last expression can be rewritten as

*<sup>∂</sup><sup>t</sup>* · **<sup>B</sup>** *dV* <sup>−</sup>

where *S* is the surface that encloses *V*. If field lines do not penetrate the volume *V* the surface integral in Eq. (33) vanishes. In the absence of charge separation **E** = −*∂***A**/*∂t*, and using

> *V*

It is clear that in the absence of resistivity (the ideal limit) *W* and *H* are conserved and, for this reason, it is said that they are ideal invariants. By contrast, in a real plasma both energy and helicity will decay at a rate proportional to *η*. However, when turbulence is present, magnetic fluctuations produce many thin current sheets with thicknesses of order *η*1/2 and

*B*<sup>2</sup> *dV*

*V B*2*dV*

*∂***A**

 *S* **A** × *∂***A**

*<sup>∂</sup><sup>t</sup>* · **<sup>B</sup>** <sup>+</sup> **<sup>A</sup>** ·∇× *<sup>∂</sup>***<sup>A</sup>**

<sup>2</sup> *dV*. (31)

*∂t* 

**J** · **B** *dV*. (34)

*<sup>∂</sup><sup>t</sup>* · **ds** (33)

*dV*. (32)

#### **4.2 Localized magnetic reconnection**

Magnetic reconnection is ubiquitous in almost all space and laboratory plasmas. It consists in a rearrangement of the topology of the magnetic field due to a change in the connectivity of the magnetic field lines. This process plays an important role in several confinement devices (such us the tokamak, the RFP and the spheromak) as well as in several astrophysical phenomena (Earth magnetosphere, solar corona, solar wind, accretion disks, etc.). Since the majority of those plasmas have a very high magnetic Reynolds number (and a high Lundquist number as well) the ideal MHD model should provide an adequate level of description for the physical system. However, as already mentioned in Sec. 3.3, topological changes in the magnetic field are not allowed in the ideal limit. What actually happens is that the coupled non linear evolution of the magnetic field and the flow inevitably develops *current sheets*, i.e. localized regions where the magnetic field gradients become very large. Within these highly localized regions the ideal approximation breaks down and the last term of Eq. (12) becomes relevant causing the magnetic field to diffuse and change the connectivity of the field lines.

The fundamental *ansatz* of the plasma relaxation theory is that these localized reconnection events do not change the total helicity of the system. Even when dissipation is involved in this process, it is assumed that only magnetic energy is affected. How can magnetic helicity be conserved during a localized reconnection event is schematically shown in Fig. 4. Two

Fig. 4. (a) Two linked ribbon-like flux tubes undergo a localized reconnection process that give rise to two separate but twisted tubes (b). The global helicity of the system is conserved.

untwisted flux tubes that are initially linked can be locally reconnected giving rise to a pair of separate but twisted tubes, in such a way that the total helicity of the system is conserved.

Moreover, plasma relaxation is based on the fact that localized reconnection events allow topological changes and dissipate magnetic energy much faster than the total helicity. There exists a number of arguments to justify this behavior (see Sec. 9.1.1 of Priest & Forbes (2000) or Montgomery et al. (1978)). Let's consider, for instance, how these two quantities (magnetic energy and helicity) decay in the presence of a small uniform resistivity *η*. Magnetic energy, *W* = *<sup>V</sup> <sup>B</sup>*2/2 *dV*, decays at the rate

14 Will-be-set-by-IN-TECH

Magnetic reconnection is ubiquitous in almost all space and laboratory plasmas. It consists in a rearrangement of the topology of the magnetic field due to a change in the connectivity of the magnetic field lines. This process plays an important role in several confinement devices (such us the tokamak, the RFP and the spheromak) as well as in several astrophysical phenomena (Earth magnetosphere, solar corona, solar wind, accretion disks, etc.). Since the majority of those plasmas have a very high magnetic Reynolds number (and a high Lundquist number as well) the ideal MHD model should provide an adequate level of description for the physical system. However, as already mentioned in Sec. 3.3, topological changes in the magnetic field are not allowed in the ideal limit. What actually happens is that the coupled non linear evolution of the magnetic field and the flow inevitably develops *current sheets*, i.e. localized regions where the magnetic field gradients become very large. Within these highly localized regions the ideal approximation breaks down and the last term of Eq. (12) becomes relevant

causing the magnetic field to diffuse and change the connectivity of the field lines.

The fundamental *ansatz* of the plasma relaxation theory is that these localized reconnection events do not change the total helicity of the system. Even when dissipation is involved in this process, it is assumed that only magnetic energy is affected. How can magnetic helicity be conserved during a localized reconnection event is schematically shown in Fig. 4. Two

(a) H = 2 Φ<sup>2</sup> (b) H = H1+H2 = 2Φ<sup>2</sup>

Fig. 4. (a) Two linked ribbon-like flux tubes undergo a localized reconnection process that give rise to two separate but twisted tubes (b). The global helicity of the system is conserved. untwisted flux tubes that are initially linked can be locally reconnected giving rise to a pair of separate but twisted tubes, in such a way that the total helicity of the system is conserved.

Moreover, plasma relaxation is based on the fact that localized reconnection events allow topological changes and dissipate magnetic energy much faster than the total helicity. There exists a number of arguments to justify this behavior (see Sec. 9.1.1 of Priest & Forbes (2000) or Montgomery et al. (1978)). Let's consider, for instance, how these two quantities (magnetic energy and helicity) decay in the presence of a small uniform resistivity *η*. Magnetic energy,

H2 = Φ<sup>2</sup>

H1 = Φ<sup>2</sup>

After reconnection

**4.2 Localized magnetic reconnection**

Reconnection

$$\frac{dW}{dt} = -\eta \int\_{V} J^2 \,dV.\tag{31}$$

This expression can be obtained by scalar multiplying Eq. (10) by **B** and integrating over a fixed volume at whose boundary the normal component of the Poynting vector vanishes. For the magnetic helicity we take the time derivative of (20) and obtain

$$\frac{dH}{dt} = \int\_{V} \left( \frac{\partial \mathbf{A}}{\partial t} \cdot \mathbf{B} + \mathbf{A} \cdot \frac{\partial \mathbf{B}}{\partial t} \right) dV = \int\_{V} \left( \frac{\partial \mathbf{A}}{\partial t} \cdot \mathbf{B} + \mathbf{A} \cdot \nabla \times \frac{\partial \mathbf{A}}{\partial t} \right) dV. \tag{32}$$

Using vector identities and the divergence theorem, the last expression can be rewritten as

$$\frac{dH}{dt} = 2 \int\_{V} \frac{\partial \mathbf{A}}{\partial t} \cdot \mathbf{B} \, dV - \int\_{S} \mathbf{A} \times \frac{\partial \mathbf{A}}{\partial t} \cdot \mathbf{ds} \tag{33}$$

where *S* is the surface that encloses *V*. If field lines do not penetrate the volume *V* the surface integral in Eq. (33) vanishes. In the absence of charge separation **E** = −*∂***A**/*∂t*, and using Ohm's law **E** = *η***J**, we finally obtain

$$\frac{dH}{dt} = -2\eta \int\_{V} \mathbf{J} \cdot \mathbf{B} \,dV. \tag{34}$$

It is clear that in the absence of resistivity (the ideal limit) *W* and *H* are conserved and, for this reason, it is said that they are ideal invariants. By contrast, in a real plasma both energy and helicity will decay at a rate proportional to *η*. However, when turbulence is present, magnetic fluctuations produce many thin current sheets with thicknesses of order *η*1/2 and current densities proportional to *Bη*<sup>−</sup>1/2. In such case, the energy decay rate becomes

$$\frac{dW}{dt} \propto \int\_{V} \, ^B \boldsymbol{B}^2 \, dV$$

which is independent of *η*. On the other hand, the total helicity decays as

$$\frac{dH}{dt} \propto 2\eta^{1/2} \int\_{V} B^2 dV$$

so that as *η* tends to zero the helicity dissipation becomes negligible. Therefore, it is important to keep in mind that a plasma will relax (in the sense described here) only if there is a certain level of small scale fluctuations that gives rise to many localized current sheets. For this reason, relaxation theory does not apply to devices with a very low level of magnetic fluctuations, such as the tokamak.

Localized magnetic reconnection events may redistribute currents in the plasma by helicity transfer between flux tubes. Even when this idea must be applied with care because the helicity is by definition a global quantity, it is clear that the helicity of a single flux tube may change after reconnection with another flux tube. This helicity transfer process is certainly at work during toroidal current drive in spheromaks and other devices sustained by helicity injection.

since

*<sup>V</sup>* **<sup>B</sup>** · ∇ *f dV* <sup>=</sup>

allowed value of *λ<sup>n</sup>* (i.e. *λ*1).

*<sup>V</sup>* ∇ · (*f***B**)*dV* <sup>=</sup>

*Br* = *B*<sup>0</sup>

*B<sup>θ</sup>* = −*B*<sup>0</sup>

which depends on the geometry of the flux conserver is

Eqs. (41) - (43) from four different positions.

meaning for the eigenvalue: *λ<sup>n</sup>* is proportional to the quotient *W*/*H*. For this reason it is clear that for a given amount of helicity, the minimum energy state will be given by the lowest

Dynamics of Magnetic Relaxation in Spheromaks 101

The most frequent model employed to describe a spheromak configuration is the relaxed state inside a cylindrical flux conserver. Using cylindrical coordinates, the condition **B** · **ds** = 0 means *Bz* = 0 at *z* = 0 and *z* = *h* and *Br* = 0 at *r* = *a*, where *h* and *a* are the height and the radius of the cylinder. In this case the solution to Eq. (39) can be found analytically (Bellan,

where *γ*<sup>1</sup> = *x*11/*a*, *k*<sup>1</sup> = *π*/*h* and *x*<sup>11</sup> is the first zero of *J*1. Note that, since this is an eigenfunction (of the curl operator) it is defined up to a constant *B*0. Note also that this solution has no toroidal dependence (i.e. it is axisymmetric). The corresponding eigenvalue

2000). In terms of Bessel functions and trigonometric functions the solution is

*π*

*λ*1 *γ*1

*λ*<sup>1</sup> =

 *x*2 11 *<sup>a</sup>*<sup>2</sup> <sup>+</sup>

Fig. 5. Four magnetic field lines showing four nested magnetic flux surfaces. This fully

relaxed state has the same value of *λ* (equal to *λ*1) on each surface.

In Fig. 5 we show the magnetic field lines obtained after following the trajectories given by

*π*2

*<sup>∂</sup>V*(*f***B**) · **ds** = 0. Eq. (40) gives us an important

*<sup>γ</sup>*1*<sup>h</sup> <sup>J</sup>*1(*γ*1*r*) cos(*k*1(*<sup>z</sup>* <sup>−</sup> *<sup>h</sup>*)) (41)

*Bz* = −*B*<sup>0</sup> *J*0(*γ*1*r*) sin(*k*1(*z* − *h*)) (43)

*J*1(*γ*1*r*) sin(*k*1(*z* − *h*)) (42)

*<sup>h</sup>*<sup>2</sup> . (44)

#### **4.3 Relaxed states**

The magnetic relaxation theory is developed for systems in which the magnetic forces are dominant, i.e. whenever the parameter *β* is low. In such cases, the MHD equilibrium Eq. (1) reduces to the force-free condition

$$
\nabla \times \mathbf{B} = \lambda(\mathbf{r}) \cdot \mathbf{B} \tag{35}
$$

where *λ*(**r**) is some scalar function. As discussed in the previous Section, magnetic fluctuations induce localized reconnection events that relax the plasma toward the state of minimum magnetic energy maintaining the total helicity of the system (Taylor, 1974). Woltjer (1958) has shown that force-free fields with *λ* equal to a constant represent the state of lowest magnetic energy under the constraint of magnetic helicity conservation in a closed system (i.e. with no field lines intercepting the boundary). The proof uses the method of Lagrange multipliers. At a constrained minimum, the variation of magnetic energy is equal to a constant (the Lagrange multiplier) times the variation of helicity

$$
\delta W = \frac{\lambda}{2} \delta H \tag{36}
$$

where *λ*/2 is the Lagrange multiplier. Substituting *W* = *<sup>V</sup> <sup>B</sup>*2/2 *dV* for the magnetic energy and Eq. (20) for *H* yields

$$\int\_{V} \left[ 2 \, \mathbf{B} \cdot \delta \mathbf{B} - \lambda (\delta \mathbf{A} \cdot \mathbf{B} + \mathbf{A} \cdot \delta \mathbf{B}) \right] dV = 0. \tag{37}$$

Using the identities

$$\mathbf{B} \cdot \delta \mathbf{B} = \delta \mathbf{A} \cdot \nabla \times \mathbf{B} - \nabla \cdot (\mathbf{B} \times \delta \mathbf{A})$$

and

$$\mathbf{A} \cdot \delta \mathbf{B} = \mathbf{B} \cdot \delta \mathbf{A} + \nabla \cdot (\delta \mathbf{A} \times \mathbf{A})$$

and the divergence theorem in Eq. (37) one obtains

$$\int\_{V} \mathbf{2} \left( \nabla \times \mathbf{B} - \lambda \mathbf{B} \right) \cdot \delta \mathbf{A} \, d^3 r = 0 \tag{38}$$

where we omitted the surface integrals because they vanish in the absence of field lines penetrating the volume under consideration. Since *δ***A** is arbitrary, the parenthesis of the integrand of Eq. (38) must be identically zero, which finally gives us the linear force-free condition

$$
\nabla \times \mathbf{B} = \lambda \mathbf{B} \tag{39}
$$

where *λ* is a constant. When we impose **B** · **ds** = 0 at the boundary, we obtain an eigenvalue problem that has non trivial solution only for certain discrete values *λ<sup>n</sup>* (which are real and positive).

Since ∇ × **B** = *λn***B**, we can write the magnetic field as **B** = *λn***A** + ∇ *f* , where *f* is an arbitrary potential. Thus, we can compute the magnetic energy as

$$W = \frac{1}{2} \int\_{V} \mathbf{B} \cdot (\lambda\_{\text{fl}} \mathbf{A} + \nabla f)dV = \frac{\lambda\_{\text{fl}}}{2} \int\_{V} \mathbf{B} \cdot \mathbf{A} \,dV = \frac{\lambda\_{\text{fl}}}{2}H \tag{40}$$

16 Will-be-set-by-IN-TECH

The magnetic relaxation theory is developed for systems in which the magnetic forces are dominant, i.e. whenever the parameter *β* is low. In such cases, the MHD equilibrium Eq. (1)

where *λ*(**r**) is some scalar function. As discussed in the previous Section, magnetic fluctuations induce localized reconnection events that relax the plasma toward the state of minimum magnetic energy maintaining the total helicity of the system (Taylor, 1974). Woltjer (1958) has shown that force-free fields with *λ* equal to a constant represent the state of lowest magnetic energy under the constraint of magnetic helicity conservation in a closed system (i.e. with no field lines intercepting the boundary). The proof uses the method of Lagrange multipliers. At a constrained minimum, the variation of magnetic energy is equal to a constant

*<sup>δ</sup><sup>W</sup>* <sup>=</sup> *<sup>λ</sup>*

2 **B** · *δ***B** − *λ*(*δ***A** · **B** + **A** · *δ***B**)

**B** · *δ***B** = *δ***A** ·∇× **B** −∇· (**B** × *δ***A**)

**A** · *δ***B** = **B** · *δ***A** + ∇ · (*δ***A** × **A**)

where we omitted the surface integrals because they vanish in the absence of field lines penetrating the volume under consideration. Since *δ***A** is arbitrary, the parenthesis of the integrand of Eq. (38) must be identically zero, which finally gives us the linear force-free

where *λ* is a constant. When we impose **B** · **ds** = 0 at the boundary, we obtain an eigenvalue problem that has non trivial solution only for certain discrete values *λ<sup>n</sup>* (which are real and

Since ∇ × **B** = *λn***B**, we can write the magnetic field as **B** = *λn***A** + ∇ *f* , where *f* is an arbitrary

2 *V*

**<sup>B</sup>** · (*λn***<sup>A</sup>** <sup>+</sup> <sup>∇</sup> *<sup>f</sup>*)*dV* <sup>=</sup> *<sup>λ</sup><sup>n</sup>*

∇ × **B** = *λ*(**r**) **B** (35)

<sup>2</sup> *<sup>δ</sup><sup>H</sup>* (36)

<sup>2</sup> (∇ × **<sup>B</sup>** <sup>−</sup> *<sup>λ</sup>***B**) · *<sup>δ</sup>***<sup>A</sup>** *<sup>d</sup>*3*<sup>r</sup>* <sup>=</sup> <sup>0</sup> (38)

∇ × **B** = *λ***B** (39)

**<sup>B</sup>** · **<sup>A</sup>** *dV* <sup>=</sup> *<sup>λ</sup><sup>n</sup>*

<sup>2</sup> *<sup>H</sup>* (40)

*<sup>V</sup> <sup>B</sup>*2/2 *dV* for the magnetic energy

*dV* = 0. (37)

**4.3 Relaxed states**

reduces to the force-free condition

and Eq. (20) for *H* yields

Using the identities

and

condition

positive).

(the Lagrange multiplier) times the variation of helicity

where *λ*/2 is the Lagrange multiplier. Substituting *W* =

 *V* 

and the divergence theorem in Eq. (37) one obtains

potential. Thus, we can compute the magnetic energy as

*<sup>W</sup>* <sup>=</sup> <sup>1</sup> 2 *V*  *V* since *<sup>V</sup>* **<sup>B</sup>** · ∇ *f dV* <sup>=</sup> *<sup>V</sup>* ∇ · (*f***B**)*dV* <sup>=</sup> *<sup>∂</sup>V*(*f***B**) · **ds** = 0. Eq. (40) gives us an important meaning for the eigenvalue: *λ<sup>n</sup>* is proportional to the quotient *W*/*H*. For this reason it is clear that for a given amount of helicity, the minimum energy state will be given by the lowest allowed value of *λ<sup>n</sup>* (i.e. *λ*1).

The most frequent model employed to describe a spheromak configuration is the relaxed state inside a cylindrical flux conserver. Using cylindrical coordinates, the condition **B** · **ds** = 0 means *Bz* = 0 at *z* = 0 and *z* = *h* and *Br* = 0 at *r* = *a*, where *h* and *a* are the height and the radius of the cylinder. In this case the solution to Eq. (39) can be found analytically (Bellan, 2000). In terms of Bessel functions and trigonometric functions the solution is

$$B\_r = B\_0 \frac{\pi}{\gamma\_1 h} J\_1(\gamma\_1 r) \cos(k\_1 (z - h))\tag{41}$$

$$B\_{\theta} = -B\_0 \frac{\lambda\_1}{\gamma\_1} f\_1(\gamma\_1 r) \sin(k\_1 (z - h))\tag{42}$$

$$B\_2 = -B\_0 l\_0(\gamma\_1 r) \sin(k\_1 (z - h))\tag{43}$$

where *γ*<sup>1</sup> = *x*11/*a*, *k*<sup>1</sup> = *π*/*h* and *x*<sup>11</sup> is the first zero of *J*1. Note that, since this is an eigenfunction (of the curl operator) it is defined up to a constant *B*0. Note also that this solution has no toroidal dependence (i.e. it is axisymmetric). The corresponding eigenvalue which depends on the geometry of the flux conserver is

$$
\lambda\_1 = \sqrt{\frac{\mathbf{x}\_{11}^2}{a^2} + \frac{\pi^2}{h^2}}.\tag{44}
$$

In Fig. 5 we show the magnetic field lines obtained after following the trajectories given by Eqs. (41) - (43) from four different positions.

Fig. 5. Four magnetic field lines showing four nested magnetic flux surfaces. This fully relaxed state has the same value of *λ* (equal to *λ*1) on each surface.

*Jz* = (∇ × **<sup>B</sup>**)*<sup>z</sup>* <sup>=</sup> <sup>1</sup>

 *r* 0

A more interesting case is obtained by setting *λ* = *λ<sup>n</sup>* (constant) which gives

*λ*(*ψ*) = *λ*¯

listed in Table 1. Note that each *α* value uniquely defines a configuration.

 1 + *α* <sup>2</sup> *<sup>ψ</sup> ψma* − 1

which is a linear *λ*(*ψ*) profile with slope *α* and mean value *λ*¯ . When this linear profile is injected in Eq. (48) a generalized non-linear eigenvalue problem is obtained. Some mathematical considerations as well as a basic numerical scheme to solve this problem were given by Kitson & Browning (1990). Note that even if one is able to solve the non-linear Grad-Shafranov equation, the profile given by Eq. (50) includes *ψma* which is not know *a priori*. The procedure adopted here is to set *ψma* = 1, fix the desired value of *α* and iterate over *λ*¯ until *ψ* is equal to one at the magnetic axis. With this procedure we obtain the values of *λ*¯

In Fig. 6 (a) we show two linear *λ*(*ψ*) profiles and (b) *ψ* contours and the *λ* colormap for the *α* = −0.4 case. The reason why we have chosen negative values for *α* is the following. It is evident that for negative values of the slope the configuration will have larger *λ* values in the outer flux surfaces (at lower *ψ* values) and *vice versa*. Since *λ* is proportional to the current

Using this, along with Eq. (3) we obtain

conditions (*ψ*|*∂*<sup>Ω</sup> = 0) are applied.

Grad-Shafranov operator.

Eq. (45) as

*<sup>B</sup><sup>θ</sup>* <sup>=</sup> <sup>1</sup> 2*πr*

*∂*2*ψ <sup>∂</sup>r*<sup>2</sup> <sup>−</sup> <sup>1</sup> *r ∂ψ ∂r* + *∂*2*ψ*

In this study we will consider initial equilibria having

*r ∂ ∂r*

Dynamics of Magnetic Relaxation in Spheromaks 103

*<sup>λ</sup>Bz* <sup>2</sup>*πrd*˜ *<sup>r</sup>*˜ <sup>=</sup> <sup>1</sup>

Thus, expressing the magnetic field in terms of *ψ* we can rewrite the toroidal component of

*<sup>∂</sup>z*<sup>2</sup> <sup>+</sup> *<sup>λ</sup>*(*ψ*)

which is the force-free version of the Grad-Shafranov equation. We are interested in solving Eq. (48) in the rectangle Ω : (*r*, *z*)=[0, *a*] × [0, *h*], i.e. a cylinder of radius *a* and height *h*.

The most simple option for *λ*(*ψ*) would be *λ* = 0, which corresponds to the vacuum solution (currentless magnetic field). The solution vanishes in this case if homogeneous boundary

<sup>−</sup> <sup>Δ</sup>∗*<sup>ψ</sup>* <sup>=</sup> *<sup>λ</sup>*<sup>2</sup>

where we have introduced the Grad-Shafranov operator defined as <sup>Δ</sup><sup>∗</sup> <sup>=</sup> *<sup>∂</sup>*2/*∂r*<sup>2</sup> <sup>−</sup> (1/*r*)*∂*/*∂r* + *∂*2/*∂z*2. If we impose homogeneous boundary conditions we obtain an eigenvalue problem which has non trivial solutions only for a discrete set of real and positive values of *λn*. The lowest value (*λ*1) is given by Eq. (44) and its associated eigenfunction is the minimum energy state described in detail in Sec. 4.3. Thus, if the appropriate boundary conditions are imposed, we can also regard the spheromak as the lowest eigenfunction of the

2*πr*

 *ψ* 0

 *ψ* 0

(*rB<sup>θ</sup>* ) = *λBz*. (46)

*λ*(*ψ*˜)*dψ*˜. (47)

*λ*(*ψ*˜)*dψ*˜ = 0 (48)

*<sup>n</sup>ψ* (49)

(50)
