Magnetic axis safety factor of finite \beta spheromaks and transition from spheromaks  to toroidal magnetic bubbles

# Magnetic axis safety factor of finite $β$ spheromaks and transition from spheromaks to toroidal magnetic bubbles

## Abstract

The value of the safety factor on the magnetic axis of a finite-beta spheromak is shown to be a function of beta in contrast to what was used in P. M. Bellan, Phys. Plasmas 9, 3050 (2002); this dependence on beta substantially reduces the gradient of the safety factor compared to the previous calculation. The method for generating finite-beta spheromak equilibria is extended to generate equilibria describing toroidal magnetic “bubbles” where the hydrodynamic pressure on the magnetic axis is less than on the toroid surface. This ”anti-confinement” configuration can be considered an  equilibrium with an inverted beta profile and is relevant to interplanetary magnetic clouds as these clouds have lower hydrodynamic pressure in their interior than on their surface.

one two three

## I Introduction

In Ref. Bellan (2002), one of the authors (PMB) examined analytic forms of finite spheromak equilibria and used a well-known expression for the value of the safety factor on the magnetic axis, denoted as ,to argue that finite causes the beneficial effect of a much larger gradient than when . However, co-author (RP) numerically calculated for these finite analytic equilibria and found numerical results substantially different from the given in Ref.Bellan (2002). The reasons for this difference are identified as resulting from a subtle misuse of an expression for Resolution of this issue revealed that the analytic equilibria presented in Ref.Bellan (2002) could be extended to give an interesting toroidal equilibria where the pressure on the magnetic axis of a toroid is lower than the pressure at the surface (edge) of the toroid rather than higher as in a tokamak; i.e., the beta profile is inverted and the configuration is bubble-like. Increase of a parameter (to be defined below) results in solutions to a Grad-Shafranov equation evolving from characterizing finite spheromak equilibria, to a conventional zero spheromak, to magnetic “bubbles” which are tokamak-like configurations having inverted profiles, and then to a tokamak with conventional profile. This evolution is characterized by the ratio of two Bessel functions changing sign as their argument is progressively increased. Interplanetary magnetic clouds are an example of the magnetic bubble situation because on the magnetic axis these clouds have lower hydrodynamic pressure than at their edge. Magnetic clouds have been previously modeled using numerical solutions to Grad-Shafranov equations Sonnerup et al. (2006),Hu et al. (2013) in a slab approximation (i.e., equations are solved in Cartesian geometry in the plane with the direction ignorable); the model presented here differs by being analytic and axisymmetric (i.e., equations are solved in cylindrical geometry in the plane with the direction ignorable) so that, in contrast to a slab approximation, toroidal geometry effects are inherently included. The analytic model has only a few parameters and so has less freedom than a numerical model but nevertheless has the useful feature of revealing parametric dependence and scaling. The analytic model also offers the possibility of providing a useful framework for other calculations, for example, calculating particle orbits in an axisymmetric cloud; the virtues of developing a repertoire of analytic solutions to the Grad-Shafranov equation has been discussed in Ref.Cerfon and Freidberg (2010).

## Ii Basic relations

We use  a cylindrical coordinate system and consider the general axisymmetric magnetic field

 B=12π∇ψ×∇ϕ+Bϕr∇ϕ (1)

where is the poloidal flux function and is the toroidal field. The direction is called the toroidal direction and any direction lying in the poloidal plane ( plane) is called a poloidal direction. From Ampere’s law the associated current density is

 J=12πμ0∇(Bϕr)×∇ϕ−r2∇ϕ2πμ0∇⋅(1r2∇ψ). (2)

We are interested in configurations where the poloidal flux function has a local extremum in the plane; both spheromaks and tokamaks are this type of configuration. The location of this extremum is called the magnetic axis and its vertical location defines the origin while its radial location is defined to be is thus at a maximum or a minimum at If is at a maximum on the magnetic axis then is positive at the axis whereas if is at a minimum on the magnetic axis then is negative at the axis.

Spheromaks and tokamaks are characterized by the safety factor which is the number of times a field line goes around toroidally for each time it goes poloidally around the magnetic axis. Tokamaks typically have near-unity on the magnetic axis with increasing with increasing distance from the magnetic axis whereas spheromaks have near-unity on the magnetic axis and decreasing on moving away from the magnetic axis. The gradient of denoted as provides stability properties and detailed calculations show that a zero spheromak has small

The safety factor at the magnetic axis is given by Bellan (2000)

 qaxis=  e1/2+e−1/2raxisBϕ,axisμ0Jϕ,axis (3)

where

 e=(ψrrψzz)axis (4)

is a measure of the ellipticity of in the vicinity of the magnetic axis such that corresponds to vertically elongated equilibria (prolate) while corresponds to vertically shortened equilibria (oblate). The force-free relation was invoked in Ref. Bellan (2000) to give but this result is valid only if the plasma is indeed force-free (i.e., has zero and equilibrium given by ). If is finite, then  and it is necessary to calculate the actual value of by consideration of the details of the finite equilibrium.

To do this, we start by defining

 βrel=μ0Paxis−PlcB2axis (5)

where and are respectively the hydrodynamic pressures on the magnetic axis and on the last closed flux surface. Positive thus corresponds to a conventional profile whereas negative corresponds to an inverted profile. This definition differs from that used in Ref. Bellan (2002) because (i) here is used and (ii) a relative rather than absolute pressure is used. The definition in Ref. Bellan (2002) used, in contrast, the average poloidal field linking the circular surface lying in the plane between the geometric axis and the magnetic axis. Because the definition of uses the relative hydrodynamic pressure it is seen that can be positive or negative. In particular, if is smaller then then will be negative The definition of is useful because it provides a simple mathematical way to distinguish toroidal equilibria with inverted profiles from those with normal profiles. The former are toroidal magnetic bubbles while the latter are toroidal confinement configurations such as spheromaks and tokamaks.

On expressing the magnetic field as

 B=12π(∇ψ(r,z)×∇ϕ+μ0I(r,z)∇ϕ) (6)

where is the poloidal current, MHD equilibrium can be expressed as the Grad-Shafranov equation Grad and Rubin (1958); Shafranov (1966)

 r∂∂r(1r∂ψ∂r)+∂2ψ∂z2+4π2μ0r2dPdψ+μ20IdIdψ=0. (7)

We assume that is a linear function of the poloidal flux and so can be expressed as

 P=Paxis−Plc(ψaxis−ψlc)ψ−Paxisψlc−Plcψaxis(ψaxis−ψlc) (8)

where is the last closed flux surface of the configuration.

The poloidal current is similarly assumed to be a linear function of the poloidal flux and can be expressed as

 μ0I=λψ. (9)

We note that the assumed linear dependence in Eq.9 differs from the assumption used in Solov’ev-type solutions such as in Ref.Cerfon and Freidberg (2010)  where it is assumed that For the linear dependence assumed here, is linear in whereas for the Solov’ev-type assumption, is a constant.

Using Eq.9, the toroidal component of Eq.6 gives

 Bϕ,axis=λψaxis2πraxis. (10)

The gradient of with respect to can then be expressed in terms of as

 dPdψ=βrelB2axisμ0ψaxis(1−ψlc/ψaxis). (11)

Since vanishes at the magnetic axis, and so

 dPdψ=βrelλ2ψaxis4π2r2axisμ0(1−ψlc/ψaxis). (12)

## Iii Cylindrical Solutions to Finite βrel Grad-Shafranov Equation

We now introduce dimensionless quantities

 ¯ψ=ψψaxis, ¯r=rraxis, ¯z=zraxis, ¯λ=λraxis (13)

so Eq.7 can be expressed as

 ¯r∂∂¯r(1¯r∂¯ψ∂¯r)+ ∂2¯ψ∂¯z2+ ¯λ2(βrel¯r2 1−¯ψlc+ ¯ψ)=0. (14)

We define

 ¯χ=βrel¯r2(1−¯ψlc)+¯ψ (15)

so Eq.14 becomes

 ¯r∂∂¯r(1¯r∂¯χ∂¯r)+ ∂2¯χ∂¯z2+ ¯λ2¯χ=0. (16)

We assume a solution of the form

 ¯χ=¯rg(¯r)cos(¯k¯z) (17)

so Eq.16 becomes

 ∂2g∂¯r2+1¯r∂g∂¯r+(¯γ2−1¯r2)g=0 (18)

where

 ¯γ2=¯λ2−¯k2. (19)

Equation 18 is Bessel’s equation with general solution for real

 g(¯r)=σJJ1(¯γ¯r)+σYY1(¯γ¯r) (20)

where and are constant coefficients to be determined by boundary conditions.

From Eqs.15 and 17 the solution to the normalized Grad-Shafranov equation is

 ¯ψ=¯r(σJJ1(¯γ¯r)+σYY1(¯γ¯r))cos(¯k¯z)−¯r2Q (21)

where

 Q=βrel1−¯ψlc. (22)

However, is required at (i.e., at the magnetic axis) so

 σJJ1(¯γ)+σYY1(¯γ)=1+Q . (23)

The following three Bessel identities where or will now be used repeatedly in the rest of the discussion:

 dC0(s)ds =−C1(s) (24a) dds(sC1(s)) =sC0(s) (24b) sC2(s) =2C1(s)−sC0(s) . (24c) The magnetic axis is also where ∂¯ψ/∂¯r vanishes and so taking the derivative of Eq.21 with respect to ¯r, using Eq.24b, and then setting ¯r=1 and ¯z=0 gives
 σJJ0(¯γ)+σYY0(¯γ)=2¯γQ. (25)

Equations 23 and 25 constitute two linear inhomogeneous algebraic equations for the coefficients and Solving these equations for and and using the Wronskian

 J1(s)Y0(s)−J0(s) Y1(s) =2πs (26)

and Eq.24c gives

 σJ =π¯γ2(Y0(¯γ)−QY2(¯γ)) (27a) σY =π¯γ2(−J0(¯γ)+QJ2(¯γ)). (27b)
 Q=J0(¯γ)J2(¯γ). (28)

Substituting for in Eq.27a and using Eqs.24c and 26 gives

 σJ=2¯γJ2(¯γ). (29)

Using Eq.24c to substitute for in Eq.28 shows that Eq.28 can alternately be written as

 Q=2J1(¯γ)¯γJ2(¯γ)−1 (30)

so one can also write as

 σJ= 1+QJ1(¯γ). (31)

Because for a spheromak Eqs.22 and 28 show that a spheromak has

 βrel=J0(¯γ)J2(¯γ)(1−¯ψlc) (32)

and

 σJ=1J1(¯γ)(1+βrel1−¯ψlc). (33)

On substituting for and in Eq.21 the solution to the normalized Grad-Shafranov equation becomes

 ¯ψ=¯rJ1(¯γ¯r)J1(¯γ)(1+βrel1−¯ψlc)cos(¯k¯z)−¯r2βrel1−¯ψlc (34)

If and are additionally assumed, the standard result for a zero-beta spheromak in a cylindrical flux conserver of radius is retrieved, namely where is the first root of Since and the last closed flux surface is at the cylinder radius, then the assumption and in Eq. 34 implies in which case   where is the first root of Thus, for a spheromak, as is well known. Equation 32 shows that spheromaks with finite positive are restricted to the range but, as will be discussed in Sec.VI, physically relevant non-spheromak configurations with negative exist when

Substitution of Eq.32 into Eq.34 gives

 ¯ψ=1J2(¯γ)(2¯r¯γJ1(¯γ¯r)cos(¯k¯z)−¯r2J0(¯γ)) (35)

which reverts to the solution when as can be seen using Eq.24c to give if .

## V Safety Factor of Spheromaks with finite βrel

The last closed flux surface of a spheromak has and in which case Eqs. 5 and 10 give

 βrel=4π2r2axisμ0Paxisλ2ψ2axis (36)

and Eq.34 becomes

 ¯ψ=¯rJ1(¯γ¯r)J1(¯γ)(1+βrel)cos(¯k¯z)−βrel¯r2 (37)

which is the same as Eq.(2) of Ref. Bellan (2002) except for the different definition of

In order to determine Eq.3  shows that it is necessary to calculate Equation 2 shows that

 μ0Jϕ=−ψaxis2π¯rr3axis[¯r∂∂¯r(1¯r∂¯ψ∂¯r)+ ∂2¯ψ∂¯z2] (38)

so, using Eq.10 and Eq.14 it is seen that

 μ0Jϕ,axisBϕ,axis=(1+βrel)λ. (39)

Thus only if Inserting Eq.39 in Eq.3 gives

 qaxis=  e1/2+e−1/2¯λ(1+βrel) (40)

which differs from Eq.(30) of Ref.Bellan (2002) by having an extra and important factor of in the denominator.

From Eq.35 and use of the Bessel identities it is seen that

 (¯ψzz)axis =−2¯k2¯γ J1(¯γ)J2(¯γ) (41a) (¯ψrr)axis =−2¯γJ1(¯γ)J2(¯γ) (41b) so the ellipticity is
 e=¯γ2¯k2. (42)

This indicates that the poloidal flux surfaces will be circular near the magnetic axis (i.e., have ) if in which case Combination of Eqs.32, 40, and 42 gives

 qaxis= ¯λ¯γ¯k11+βrel= ¯λ2¯kJ2(¯γ)J1(¯γ). (43)

Equation 43 has been validated by direct numerical integration of field lines in the vicinity of the magnetic axis of a magnetic configuration characterized by Eq.6 with given by Eq.35. In the limit, and which is Eq.(33) of Ref.Bellan (2002), but for finite positive Eq. 43 shows that is reduced from its value.

The safety factor at the wall is Bellan (2002)

 qwall=¯λ2π¯kcos−1(J0(¯γ)) (44)

and so the ratio of safety factor at the wall to that at the axis is

 qwallqaxis=J1(¯γ)J2(¯γ)cos−1(J0(¯γ))π (45)

which is plotted in Fig.1. Contrary to Ref.Bellan (2002) it is seen that the shear (difference between and decreases with increasing (i.e., with decreasing below 2.405). Using to write

 cos−1(J0(¯γ))=−iln(J0(¯γ)+i√1−(J0(¯γ))2) (46)

and then using for , it is seen that for

 cos−1(J0(¯γ))≃−iln(1−¯γ24+i¯γ√2)≃¯γ√2. (47)

Since and for Eq.45 has the limiting behavior

 qwallqaxis→4√2π=0.900 for ¯γ≪1 (48)

which is seen in Fig. 1. Furthermore, Eq.32 has the limiting behavior

 βrel→8¯γ2 for ¯γ≪1; (49)

i.e., diverges at small  which is also seen in Fig. 1.

We note that numerical calculations reported in Ref.Gautier et al. (1981) assumed and in a spherical geometry and found that the gradient of the shear had a strong dependence on The analytic solution given here would correspond approximately to the numerical solution reported in Ref.Gautier et al. (1981); the correspondence is not exact because of the different assumptions for the dependence of on , the shape of the boundary (cylinder v. sphere), and the assumption of a central hole in Ref. Gautier et al. (1981).

## Vi Toroidal magnetic bubble: Negative βrel

We now consider the situation where and We consider the case first as was assumed for spheromaks and then later consider the more general case where both and are finite.

### vi.1 σY=0 case

In the case is mathematically identical to the spheromak solution considered in Sec.IV, i.e., Eq.35 provides the relevant flux function. The difference here is that is no longer assumed to be zero. Plots of using show that has periodic maxima and minima because of its dependence. Equation 13 defined to be unity on the magnetic axis, i.e., at and the magnetic axis was defined to be where was a maximum or minimum. Because of the oscillatory behavior of Bessel functions, maxima or minima of occur not only at but also for However, the maxima and minima occurring where do not have and so do not satisfy the condition given in Eq.13. Thus, only the maximum of at will be considered since maxima or minima at larger do not satisfy the requirement stipulated in Eq.13.

Examination of Eq.35 shows that is independent of if at this radius where is the first root of We now show that this radius is infinitesimally larger than the radius of the last closed flux surface. Since is independent of when the flux surface passing through must be a straight vertical line, i.e., for all Because a straight vertical line goes to , the flux surface passing through is open. Immediately to the left of this line the flux surfaces are closed and so the last closed flux surface is at the radius where

 ¯rlc=limδ→0(x11¯γ−δ)=x11¯γ. (50)

This can also be seen graphically from the flux surface contours shown in Fig.2 (to be discussed in more detail later) where it is seen that a straight vertical line separatrix lies between the blue-purple closed flux surfaces having magnetic axis at and the green-orange flux surfaces to the right. Equation 50 gives the radial location of this vertical line.

A toroidal inverse aspect ratio (ratio of torus minor to major radius) can be defined as

 ε=rlc−raxisraxis=¯rlc−1=x11−¯γ¯γ. (51)

Using at the last closed flux surface, Eq. 35 may be evaluated at to give

 ¯ψlc=−x211¯γ2J0(¯γ)J2(¯γ). (52)

Inserting in Eq.32 gives

 βrel=(1+x211¯γ2J0(¯γ)J2(¯γ))J0(¯γ)J2(¯γ). (53)

In order to have Eq.51 shows that it is necessary to have  A plot of Eq.53 shows that is negative if ; changes sign at because the quantity in parenthesis in Eq.53 changes sign at Thus if , is negative and also

Because the minimum of occurs when identifying   it is seen that is at a minimum when in which case

 min[βrel]=−¯γ24x211. (54)

Using the Bessel identities, the magnetic field components are

 Br =−ψaxis2πr2axis¯r∂¯ψ∂¯z=ψaxis2πr2axis 2¯k ¯γJ1(¯γ¯r)J2(¯γ)sin(¯k¯z) (55a) Bϕ =λψ2πr= ψaxis2πr2axis  ¯λJ2(¯γ)(2¯γJ1(¯γ¯r)cos(¯k¯z)−¯rJ0(¯γ)) (55b) Bz =ψaxis2πr2axis¯r∂¯ψ∂¯r=ψaxis2πr2axis 2J2(¯γ)( J0(¯γ¯r)cos(¯k¯z)−J0(¯γ)). (55c) Using Eq.24c and Eq.55b it is seen that
 Bϕ,axis=ψaxis2πr2axis  ¯λ. (56)

A normalized magnetic field can be defined as with components

 ¯Br(¯r,¯z) =2¯k ¯γ¯λJ1(¯γ¯r)J2(¯γ)sin(¯k¯z) (57a) ¯Bϕ(¯r,¯z) =2¯γJ1(¯γ¯r)cos(¯k¯z)−¯rJ0(¯γ)J2(¯γ) (57b) ¯Bz(¯r,¯z) =2¯λ (J0(¯γ¯r)cos(¯k¯z)−J0(¯γ))J2(¯γ). (57c) As required, both ¯Br and ¯Bz vanish on the magnetic axis (i.e., at ¯r=1,¯z=0) and ¯Bϕ=1 on the magnetic axis. Equation 35 with 2.405<¯γ<3.736 thus gives the flux surface for a magnetic bubble, i.e., a toroidal configuration with closed field lines where the pressure on the magnetic axis is lower than the pressure at the surface of the toroid. The direction of the J×B force is thus outwards rather than inwards in contrast to a tokamak. This configuration is relevant to axisymmetric interplanetary magnetic clouds ejected from the sun by coronal mass ejections. Spacecraft measurements indicate that P is smaller in the interior of these clouds than outside  so these clouds have negative βrel. Another possible situation would be in the solar interior where a toroidal bubble configuration as described here would be a toroidal region of stronger magnetic field but reduced hydrodynamic pressure compared to the surroundings. As a concrete example of such a configuration, consider the situation where ¯γ=¯k=2.5 and ¯λ=√2¯γ. In this case e=1 so the poloidal flux surfaces are circular near the magnetic axis, the last closed flux surface is at ¯ψlc=0.25 and from Eq.53 βrel=−0.081. From Eq.51, it is seen that the inverse aspect ratio is ε=\allowbreak0.53. Figure 2 plots contours of ψ(¯r,¯z) and it is seen that the last closed flux surface intersects ¯z=0 to the right of the magnetic axis at indeed ¯rlc=x11/¯γ=\allowbreak1.\allowbreak53. Figures 3, 4, 5, and 6 plot ¯ψ(¯r,0), ¯Bϕ(¯r,0), ¯Bz(¯r,0), and ¯B2(¯r,0) respectively. From Eq.5 it is seen that
 μ0B2axisPaxis=μ0B2axisPlc+βrel (58)

so the hydrodynamic pressure on the magnetic axis is lower than on the last closed flux surface. If is set to zero, then the external pressure would be

 μ0B2axisPlc=−βrel (59)

in which case the configuration would be a vacuum at the magnetic axis (zero plasma pressure) with increasing pressure going away from the magnetic axis toward last closed flux surface.

If is further increased, the sign of can become positive again in which case the equilibrium will become tokamak-like (higher pressure on magnetic axis). Additional increase of will cause to oscillate in sign giving a sequence of bubble-like and tokamak-like configurations. Also, for a given configuration one could elect to truncate the flux at some value larger than and so obtain a smaller aspect ratio equilibrium. In accordance with the Shafranov virial theorem, any one of these configurations will involve a jump in the magnetic field at the surface of the toroid if it is assumed that at the surface the external magnetic field differs from the internal field. This jump corresponds to the existence of surface currents. In a tokamak these surface currents are provided by a set of coils immediately external to the toroidal volume and these coils are called the vertical field coils. The field produced by these coils is mainly in the direction and will be referred to here as . This field constitutes a portion of the total field inside the toroidal volume and provides equilibrium in the major radius direction. This takes place via a radial force directed towards that balances the radially outward hoop force as well as some hydrodynamic pressure forces. The hoop force is a property of any toroidal current system and occurs because a toroidal current produces a stronger poloidal field near , (inside) than at (outside) This stronger poloidal field on the inside compared to the outside corresponds to greater magnetic pressure on the inside than on the outside; for low the force resulting from magnetic pressure imbalance dominates any hydrodynamic pressure imbalance. Without the offsetting force provided by the hoop force would act to expand the torus major radius.

At first sight it might appear that the flux contours in Fig.2 are such that the magnetic pressure is higher on the outside than on the inside because the midplane poloidal flux surfaces in Fig. 2 are more tightly packed outside the magnetic axis (e.g., at than inside the magnetic axis (e.g., at ). However, the density of field lines and hence the poloidal field is nevertheless stronger inside the magnetic axis than outside because of the inverse dependence in . The twice as tight midplane flux surface packing in Fig. 2 at compared to at gives a twice as large on the outside compared to the inside. However, this twice as tight radial packing is overcome by the factor, a toroidal geometry effect that produces an approximately six-fold inside-to-outside enhancement with the net result that is about three times larger at than at . This three-fold inside-to-outside ratio of is evident in Fig. 5.

In order to have the required  for equilibrium, it would be necessary to have surface currents flowing on the surface of the toroid. Since there are no powered coils to sustain surface currents exterior to a magnetic cloud, it is unlikely that such surface currents would exist in the magnetic cloud context. Without the provided by surface currents  (and intrinsic to the equilibrium given here), the hoop force resulting from the imbalance between on the inside and on the outside will cause the major radius of magnetic clouds to increase with time. The difference between poloidal flux surfaces with and without incorporation of is of the order of the inverse aspect ratio because is a toroidal effect and so scales as

### vi.2 Finite σJ and σY case

The spheromak solution required to be zero to avoid singularity at The magnetic bubble solution discussed above used the same functional form as the spheromak solution (i.e., had and used Eq.34) and found that a tokamak-like solution with (i.e., inverted beta profile) occurred if If is excluded from the domain so the configuration is doubly-connected, the singular nature of at