Inverse aspect-ratio expanded tokamak equilibria

Fitzpatrick, R.

doi:10.1063/5.0215345

Following Greene et al. [Phys. Fluids 14, 671 (1971)] and Connor et al. [Phys. Plasmas 31, 577 (1988); Plasma Phys. Control. Fusion 34, 161 (1992); and Nucl. Fusion 33, 1533 (1993)], the Grad-Shafranov equation for an axisymmetric tokamak plasma equilibrium is solved via an expansion in the, supposedly small, inverse aspect-ratio of the plasma, ϵ. The displacements of equilibrium magnetic flux-surfaces due to plasma shaping are assumed to be $O (ϵ)$ smaller than the minor radii of the surfaces, but no other restriction is placed on the nature of the shaping. The solution of the Grad-Shafranov equation is matched to a vacuum solution that extends to infinity, and consists of an expansion in toroidal functions. The external poloidal magnetic field generated by a finite set of discrete external poloidal magnetic field-coils is calculated, and incorporated into the toroidal function expansion. In this manner, the shape of a large aspect-ratio tokamak plasma is directly related to the currents flowing in the external poloidal field-coils. Finally, a pedestal in the plasma pressure, and the associated spike in the bootstrap current, are incorporated into the model.

I. INTRODUCTION

In a celebrated paper, published in 1971, Greene, Johnson, and Weimer (GJW) developed a theory of axisymmetric tokamak equilibria that is based on an expansion of the Grad–Shafranov equation (see Sec. II D) in the inverse aspect-ratio of the plasma (i.e., the ratio of the minor to the major radius of the plasma torus).¹ In addition, GJW demonstrated how to match the plasma magnetic field to a curl-free magnetic field, in the vacuum region exterior to the plasma, by distinguishing between the field generated in the vacuum region by currents flowing within the plasma and the field generated by external magnetic field-coils.

GJW assumed that the horizontal (Shafranov) shift² of the plasma magnetic flux-surfaces is $O (ϵ)$ smaller than the minor radius, where $ϵ ≪ 1$ is the inverse aspect-ratio. Furthermore, GJW assumed that the ellipticity of the flux-surfaces is $O (ϵ)$ smaller than the horizontal shift. This latter assumption was reasonable for the tokamak plasmas, characterized by almost circular boundaries in the poloidal plane, that were prevalent in the 1960s and 1970s. Nowadays, however, tokamak plasmas are highly elongated and strongly shaped in the poloidal plane.³ Starting in 1985, in a series of papers, Connor et al. (CHA) generalized the GJW analysis to allow for flux-surface elongation and triangularity that is of the same order of magnitude as the horizontal shift.^4–7 However, CHA did not perform the vacuum matching procedure. Moreover, CHA assumed that the plasma is up-down symmetric. Nowadays, however, tokamak plasmas tend to have a lower magnetic null on the plasma boundary that render them significantly up-down asymmetric.³ Tokamak plasmas also generally feature a pedestal in the plasma pressure, close to the plasma boundary, with an associated spike in the parallel current density due to the bootstrap current.³

The aim of this paper is to generalize the analysis of GJW and CHA so as to allow for strong magnetic flux-surface shaping of a general nature, lack of up-down symmetry, and the presence of a pedestal in the plasma pressure profile. Furthermore, we wish to perform the matching to the vacuum solution, with the eventual aim of directly relating the shape of the plasma to currents flowing in a set of discrete external poloidal magnetic field-coils.

Despite the widespread availability of fast numerical solvers for the Grad–Shafranov equation,^8,9 analytic model equilibria are still widely used by fusion scientists,^10–12 particularly in preliminary design studies. Most of these model equilibria are extensions of the model equilibria first discovered by Solov'ev in the 1960s.¹⁰ Solov'ev-type equilibria have the property that the plasma current profile extends all the way to the plasma boundary, necessitating a large discontinuous jump in the current density across the plasma/vacuum interface. Such a jump has a highly destabilizing influence on ideal external-kink modes,¹³ which generally means that Solov'ev-type equilibria are unsuitable for free-boundary plasma stability calculations. On the other hand, the aspect-ratio expanded model equilibria of GJW and CHA allow the plasma current density to go to zero at the plasma boundary, eliminating the need for a discontinuous jump in the plasma current density across the plasma/vacuum interface. Hence, such equilibria are quite suitable for free-boundary plasma stability calculations.

Conventional numerical solvers for the Grad–Shafranov equation, such as HELENA and CHEASE,^8,9 only calculate the solution inside the plasma boundary, and, hence, do not provide the vacuum solution outside the boundary. Reference 14 demonstrates how a vacuum solution can be added to a Solov'ev solution by means of a Green's function. Moreover, Ref. 15 shows how this vacuum solution can be related to currents flowing in a finite set of discrete poloidal magnetic field-coils surrounding the plasma. Such a vacuum solution can easily be added to a HELENA or a CHEASE solution using similar techniques. Unfortunately, the numerical implementation of a general Green's function solution is extremely time consuming. One of the main advantages of the GJW and CHA model equilibria is that most of the hard work entailed in the calculation of the vacuum solution can be performed by means of analysis, so that the numerical implementation of the solution can be carried out very rapidly.

This paper is organized as follows: In Sec. II, we discuss the equilibrium solution inside the plasma. In Sec. III, we match the plasma solution to a generic vacuum solution outside the plasma. In Sec. IV, we match the plasma solution to the vacuum solution generated by a set of discrete external poloidal magnetic field-coils. Pedestal physics is discussed in Sec. V. In Sec. VI, we describe some example calculations. (Note that, given the great efficiency of the analytic expansion method described in this paper, all of these calculations can be performed in a fraction of a second on an ordinary laptop computer.) Finally, the paper is summarized, and conclusions are drawn, in Sec. VII.

II. PLASMA SOLUTION

A. Coordinates

Consider a tokamak plasma confined on a set of nested, axisymmetric, magnetic flux-surfaces. Let R,

ϕ

⁠, Z be right-handed cylindrical coordinates whose symmetry axis corresponds to the common symmetry axis of the flux-surfaces. Let us define right-handed flux-coordinates, r, θ,

ϕ

⁠, where r(R, Z) is a magnetic flux-surface label with the dimensions of length, and

θ (R, Z)

a poloidal angle. Let^4,5,7

J \equiv {(\nabla r \times \nabla θ \cdot \nabla ϕ)}^{- 1} \equiv R (\frac{\partial R}{\partial θ} \frac{\partial Z}{\partial r} - \frac{\partial R}{\partial r} \frac{\partial Z}{\partial θ}) = \frac{r R^{2}}{R_{0}}

(1)

be the Jacobian of the flux-coordinate system. Here,

(R_{0}, 0)

are the (R, Z) coordinates of the magnetic axis of the nested magnetic flux-surfaces. The axis corresponds to r = 0 in the flux-coordinate system. Note that

| \nabla ϕ | = R^{- 1}

⁠.

B. Equilibrium magnetic field and current

The divergence-free equilibrium magnetic field is written¹

B = B_{0} R_{0} [f (r) \nabla ϕ \times \nabla r + g (r) \nabla ϕ],

(2)

where f(r) and g(r) are free (dimensionless) flux-surface functions, and B₀ is the toroidal magnetic field-strength at the magnetic axis.

It is helpful to define the toroidal magnetic flux (divided by

2 π

⁠) contained within the magnetic flux-surface whose label is r,

Ψ_{t} (r) = B_{0} \int_{0}^{r} r' g (r') d r',

(3)

and the corresponding poloidal magnetic flux (divided by

2 π

⁠),

Ψ_{p} (r) = B_{0} R_{0} \int_{0}^{r} f (r') d r' .

(4)

The safety-factor of the flux-surface (which is the inverse of the rotational transform),¹⁶

q (r) \equiv \frac{d Ψ_{t}}{d Ψ_{p}} = \frac{r g}{R_{0} f},

(5)

is directly related to these two fluxes.

It follows from Eqs. , , and that

B = B_{0} R_{0} f (r) \nabla r \times \nabla (q θ - ϕ),

(6)

which implies that

\frac{B \cdot \nabla θ}{B \cdot \nabla ϕ} = \frac{1}{q (r)} .

(7)

Consequently, if the equilibrium magnetic field-lines within a given magnetic flux-surface are plotted in the θ-

ϕ

plane then they appear as parallel straight lines whose pitch is determined by the surface's safety-factor (or rotational transform).

It is also helpful to define the toroidal plasma current flowing within the flux-surface whose label is r,

I_{t} (r) = \frac{1}{2 π} \int_{0}^{r} j \cdot \nabla ϕ J d r d θ d ϕ = \frac{2 π B_{0}}{μ_{0}} r f (r) ⟨ | \nabla r |^{2} ⟩ (r),

(8)

and the corresponding poloidal plasma current,

I_{p} (r) = \frac{1}{2 π} \int_{0}^{r} j \cdot \nabla θ J d r d θ d ϕ = \frac{2 π B_{0} R_{0}}{μ_{0}} [1 - g (r)] .

(9)

Here,

j

is the equilibrium plasma current density. Moreover,

⟨ X ⟩ (r) \equiv \oint X (r, θ) d θ / 2 π

⁠. Note that

g (0) = 1

⁠, by convention, so as to ensure that B₀ is indeed the toroidal magnetic field-strength on the magnetic axis.

C. Equilibrium magnetic flux-surfaces

The loci of the equilibrium magnetic flux-surfaces are written in the parametric form^1,4–7

\begin{matrix} R (r, ω) = R_{0} - r cos ω + \sum_{n > 0} H_{n} (r) cos [(n - 1) ω] + \sum_{n > 1} V_{n} (r) sin [(n - 1) ω] + L (r) cos ω, \end{matrix}

(10)

\begin{matrix} Z (r, ω) = r sin ω + \sum_{n > 1} H_{n} (r) sin [(n - 1) ω] - \sum_{n > 1} V_{n} (r) cos [(n - 1) ω] - L (r) sin ω . \end{matrix}

(11)

Here,

H_{1} (r)

controls the relative horizontal locations of the flux-surface centroids, respectively,

H_{2} (r)

and

V_{2} (r)

control the up-down symmetric and asymmetric flux-surface ellipticities, respectively, and

H_{3} (r)

V_{3} (r)

control the up-down symmetric and asymmetric flux-surface triangularities, respectively, etc., whereas L(r) is a relabeling parameter. Moreover,

ω (R, Z)

is a poloidal angle that is distinct from θ.

Let

J (r, ω) = \frac{\partial R}{\partial ω} \frac{\partial Z}{\partial r} - \frac{\partial R}{\partial r} \frac{\partial Z}{\partial ω}

(12)

be the Jacobian of the r, ω coordinate system. We can transform to the r, θ coordinate system by writing

θ (r, ω) = 2 π \int_{0}^{ω} \frac{J (r, ω')}{R (r, ω')} d ω' / \oint \frac{J (r, ω)}{R (r, ω)} d ω,

(13)

\frac{r}{R_{0}} = \frac{1}{2 π} \oint \frac{J (r, ω)}{R (r, ω)} d ω .

(14)

This transformation ensures that

\frac{\partial θ}{\partial ω} = \frac{J R_{0}}{r R},

(15)

and, hence, that

J \equiv R (\frac{\partial R}{\partial θ} \frac{\partial Z}{\partial r} - \frac{\partial R}{\partial r} \frac{\partial Z}{\partial θ}) = R J \frac{\partial ω}{\partial θ} = \frac{r R^{2}}{R_{0}},

(16)

in accordance with Eq. .

D. Grad–Shafranov equation

Equilibrium force balance within the plasma, and the surrounding vacuum, is governed by the Grad–Shafranov equation,^1,4–7

\frac{f}{r} \frac{\partial}{\partial r} (r f | \nabla r |^{2}) + \frac{f}{r} \frac{\partial}{\partial θ} (r f \nabla r \cdot \nabla θ) = - g \frac{d g}{d r} - \frac{μ_{0}}{B_{0}^{2}} {(\frac{R}{R_{0}})}^{2} \frac{d p}{d r},

(17)

where p(r) is the equilibrium plasma pressure. Note that we can easily extend the Grad–Shafranov equation to include the centrifugal effect of plasma flows on the equilibrium.^17,18 We choose not to do this in the present calculation because such effects are generally negligible in a reactor-sized tokamak.

E. Normalization

Let r = a correspond to the plasma/vacuum boundary, and let

ϵ = \frac{a}{R_{0}}

(18)

be the inverse-aspect ratio of the plasma. It is assumed that

ϵ ≪ 1

⁠. Furthermore, let

R = R_{0} \hat{R}, Z = R_{0} \hat{Z}, r = a \hat{r}, \nabla = a^{- 1} \hat{\nabla}

⁠,

H_{n} = ϵ a {\hat{H}}_{n}, V_{n} = ϵ a {\hat{V}}_{n}

⁠, and

L = ϵ^{2} a \hat{L}

⁠. The parametric equations of the flux-surfaces, and , become^6,7

\begin{matrix} \hat{R} (\hat{r}, ω) = 1 - ϵ \hat{r} cos ω + ϵ^{2} \sum_{n > 0} {\hat{H}}_{n} (\hat{r}) cos [(n - 1) ω] \\ + ϵ^{2} \sum_{n > 1} {\hat{V}}_{n} (\hat{r}) sin [(n - 1) ω] + ϵ^{3} \hat{L} (\hat{r}) cos ω, \end{matrix}

(19)

\begin{matrix} \hat{Z} (\hat{r}, ω) = ϵ \hat{r} sin ω + ϵ^{2} \sum_{n > 1} {\hat{H}}_{n} (\hat{r}) sin [(n - 1) ω] \\ - ϵ^{2} \sum_{n > 1} {\hat{V}}_{n} (\hat{r}) cos [(n - 1) ω] - ϵ^{3} \hat{L} (\hat{r}) sin ω . \end{matrix}

(20)

It is assumed that, other than ϵ, all quantities appearing in the previous two expressions are

O (1)

⁠.

F. Metric elements

We can determine the metric elements of the flux-coordinate system by combining the previous two equations with Eqs. and . Evaluating the elements up to

O (ϵ)

⁠, but retaining

O (ϵ^{2})

contributions to terms that are independent of ω, we obtain,^6,7

\hat{L} (\hat{r}) = \frac{{\hat{r}}^{3}}{8} - \frac{\hat{r} {\hat{H}}_{1}}{2} - \frac{1}{2} \sum_{n > 1} (n - 1) \frac{{\hat{H}}_{n}^{2}}{\hat{r}} - \frac{1}{2} \sum_{n > 1} (n - 1) \frac{{\hat{V}}_{n}^{2}}{\hat{r}},

(21)

\begin{matrix} θ = ω + ϵ \hat{r} sin ω - ϵ \sum_{n > 0} \frac{1}{n} [{\hat{H}}_{n}^{'} - (n - 1) \frac{{\hat{H}}_{n}}{\hat{r}}] sin (n ω) \\ + ϵ \sum_{n > 1} \frac{1}{n} [{\hat{V}}_{n}^{'} - (n - 1) \frac{{\hat{V}}_{n}}{\hat{r}}] cos (n ω), \end{matrix}

(22)

\begin{matrix} | \hat{\nabla} \hat{r} |^{2} = 1 + 2 ϵ \sum_{n > 0} {\hat{H}}_{n}^{'} cos (n θ) + 2 ϵ \sum_{n > 1}' {\hat{V}}_{n}^{'} sin (n θ) \\ + ϵ^{2} (\frac{3 {\hat{r}}^{2}}{4} - {\hat{H}}_{1} + \frac{1}{2} \sum_{n > 0} [{\hat{H}}_{n}^{' 2} + (n^{2} - 1) \frac{{\hat{H}}_{n}^{2}}{{\hat{r}}^{2}}] \\ + \frac{1}{2} \sum_{n > 1} [{\hat{V}}_{n}^{' 2} + (n^{2} - 1) \frac{{\hat{V}}_{n}^{2}}{{\hat{r}}^{2}}]), \end{matrix}

(23)

\begin{matrix} \hat{\nabla} \hat{r} \cdot \hat{\nabla} θ = ϵ sin θ - ϵ \sum_{n > 0} \frac{1}{n} [{\hat{H}}_{n}^{″} + \frac{{\hat{H}}_{n}^{'}}{\hat{r}} + (n^{2} - 1) \frac{{\hat{H}}_{n}}{{\hat{r}}^{2}}] sin (n θ) \\ + ϵ \sum_{n > 1} \frac{1}{n} [{\hat{V}}_{n}^{″} + \frac{{\hat{V}}_{n}^{'}}{\hat{r}} + (n^{2} - 1) \frac{{\hat{V}}_{n}}{{\hat{r}}^{2}}] cos (n θ) \\ + ϵ^{2} \frac{1}{2} [\sum_{n > 1} (n - 1) (\frac{{\hat{V}}_{n}^{'} {\hat{H}}_{n}}{{\hat{r}}^{2}} - \frac{{\hat{H}}_{n}^{'} {\hat{V}}_{n}}{{\hat{r}}^{2}}) + \sum_{n > 1} (\frac{n - 1}{n}) (\frac{{\hat{V}}_{n}^{″} {\hat{H}}_{n}}{\hat{r}} - \frac{{\hat{H}}_{n}^{″} {\hat{V}}_{n}}{\hat{r}}) + \sum_{n > 1} \frac{1}{n} ({\hat{V}}_{n}^{″} {\hat{H}}_{n}^{'} - {\hat{H}}_{n}^{″} {\hat{V}}_{n}^{'})], \end{matrix}

(24)

{\hat{R}}^{2} = 1 - 2 ϵ \hat{r} cos θ - ϵ^{2} (\frac{{\hat{r}}^{2}}{2} - \hat{r} {\hat{H}}_{1}^{'} - 2 {\hat{H}}_{1}) .

(25)

Here,

' \equiv d / d \hat{r}

⁠, and n is a non-negative integer.

G. Expansion of Grad–Shafranov equation

Let us write

\hat{r} f (\hat{r}) = ϵ f_{1} (\hat{r}) + ϵ^{3} {\hat{f}}_{3} (\hat{r}) + \dots,

(26)

g (\hat{r}) = 1 + ϵ^{2} g_{2} (\hat{r}) + ϵ^{4} g_{4} (\hat{r}) + \dots,

(27)

\frac{μ_{0} p (\hat{r})}{B_{0}^{2}} = ϵ^{2} p_{2} (\hat{r}) + ϵ^{4} p_{4} (\hat{r}) + \dots,

(28)

where f₁, f₃, g₂, g₄, p₂, and p₄ are all

O (1)

⁠. It follows from Eq. that, to lowest order in our expansion, the safety-factor profile is given by

q_{0} (\hat{r}) = \frac{{\hat{r}}^{2}}{f_{1} (\hat{r})},

(29)

whereas, to second order, the profile takes the modified form

q_{2} (\hat{r}) = \frac{{\hat{r}}^{2}}{f_{1} (\hat{r})} [1 + ϵ^{2} g_{2} (\hat{r}) - ϵ^{2} \frac{f_{3} (\hat{r})}{f_{1} (\hat{r})}] .

(30)

Expanding the Grad–Shafranov equation, , order by order in the small parameter ϵ, making use of Eqs. , we obtain^1,4–6

g_{2}^{'} = - p_{2}^{'} - \frac{f_{1} f_{1}^{'}}{{\hat{r}}^{2}},

(31)

{\hat{H}}_{1}^{″} = - (\frac{2 f_{1}^{'}}{f_{1}} - \frac{1}{\hat{r}}) {\hat{H}}_{1}^{'} - 1 + \frac{2 {\hat{r}}^{3} p_{2}^{'}}{f_{1}^{2}},

(32)

{\hat{H}}_{n}^{″} = - (\frac{2 f_{1}^{'}}{f_{1}} - \frac{1}{\hat{r}}) {\hat{H}}_{n}^{'} + (n^{2} - 1) \frac{{\hat{H}}_{n}}{{\hat{r}}^{2}} for n > 1,

(33)

{\hat{V}}_{n}^{″} = - (\frac{2 f_{1}^{'}}{f_{1}} - \frac{1}{\hat{r}}) {\hat{V}}_{n}^{'} + (n^{2} - 1) \frac{{\hat{V}}_{n}}{{\hat{r}}^{2}} for n > 1,

(34)

{\hat{r}}^{2} p_{4}^{'} + {\hat{r}}^{2} g_{4}^{'} + (f_{1} f_{3})' = F_{3} (\hat{r}),

(35)

where

\begin{matrix} F_{3} (\hat{r}) = - \frac{f_{1}^{2}}{\hat{r}} (\frac{3 {\hat{r}}^{2}}{2} - 2 \hat{r} {\hat{H}}_{1}^{'} \\ + \sum_{n > 0} [{\hat{H}}_{n}^{' 2} + 2 (n^{2} - 1) \frac{{\hat{H}}_{n}^{'} {\hat{H}}_{n}}{\hat{r}} - (n^{2} - 1) \frac{{\hat{H}}_{n}^{2}}{{\hat{r}}^{2}}] \\ + \sum_{n > 1} [{\hat{V}}_{n}^{' 2} + 2 (n^{2} - 1) \frac{{\hat{V}}_{n}^{'} {\hat{V}}_{n}}{\hat{r}} - (n^{2} - 1) \frac{{\hat{V}}_{n}^{2}}{{\hat{r}}^{2}}]) + f_{1} f_{1}^{'} (g_{2} - \frac{3 {\hat{r}}^{2}}{4} + {\hat{H}}_{1} + \frac{1}{2} \sum_{n > 0} [3 {\hat{H}}_{n}^{' 2} - (n^{2} - 1) \frac{{\hat{H}}_{n}^{2}}{{\hat{r}}^{2}}] + \frac{1}{2} \sum_{n > 1} [3 {\hat{V}}_{n}^{' 2} - (n^{2} - 1) \frac{{\hat{V}}_{n}^{2}}{{\hat{r}}^{2}}]) \\ + {\hat{r}}^{2} p_{2}^{'} (g_{2} + \frac{{\hat{r}}^{2}}{2} - 3 \hat{r} {\hat{H}}_{1}^{'} - 2 {\hat{H}}_{1}) . \end{matrix}

(36)

Note that the relative horizontal shift of magnetic flux-surfaces,

{\hat{H}}_{1}

⁠, otherwise known as the Shafranov shift,² is driven by toroidicity [the second term on the right-hand side of Eq. ], and plasma pressure gradients (the third term). All of the other shaping terms (i.e., the

{\hat{H}}_{n}

for n > 1, and the

{\hat{V}}_{n}

for n > 1) are driven by currents flowing in external magnetic field-coils.

In the present study, it is convenient to set $g_{4} = p_{4} = 0$ ⁠. However, other choices are possible.

H. Plasma interior

Let us regard

f_{1} (\hat{r})

and

p_{2} (\hat{r})

as the two free flux-surface functions that determine the equilibrium. In the plasma interior, which corresponds to

0 \leq \hat{r} \leq 1

⁠, Eq. yields

g_{2} (\hat{r}) = p_{2} (0) - p_{2} (\hat{r}) - \int_{0}^{\hat{r}} \frac{f_{1} (\hat{r}') f_{1}^{'} (\hat{r}')}{{\hat{r}}^{' 2}} d \hat{r}' .

(37)

At small

\hat{r}

⁠, assuming that

f_{1} (\hat{r}) = f_{1 c} {\hat{r}}^{2}

and

p_{2} (\hat{r}) = p_{2 c} + (1 / 2) p_{2 c}^{″} {\hat{r}}^{2}

⁠, where

f_{1 c}, p_{2 c}, p_{2 c}^{″}

are constants that can be determined from the given

f_{1} (\hat{r})

and

p_{2} (\hat{r})

functions, we deduce that

g_{2} (\hat{r}) = - (f_{1 c}^{2} + p_{2 c}^{″} / 2) {\hat{r}}^{2}

⁠. Note that we have set

g_{2} = 0

at

\hat{r} = 0

in order to ensure that B₀ is indeed the toroidal magnetic field-strength on the magnetic axis. Equation must be solved subject to the boundary conditions that

{\hat{H}}_{1} (\hat{r}) = \frac{1}{8} (\frac{2 p_{2 c}^{″}}{f_{1 c}^{2}} - 1) {\hat{r}}^{2}

(38)

at small

\hat{r}

⁠. Likewise, Eqs. and must be solved subject to the small-

\hat{r}

boundary conditions that

{\hat{H}}_{n} (\hat{r}) = {\hat{H}}_{n c} {\hat{r}}^{n - 1} for n > 1,

(39)

{\hat{V}}_{n} (\hat{r}) = {\hat{V}}_{n c} {\hat{r}}^{n - 1} for n > 1.

(40)

Here,

{\hat{H}}_{n c}

and

{\hat{V}}_{n c}

are arbitrary constants. It turns out that the only solution of Eq. for n = 1 that does not blow up at

\hat{r} = 0

is the trivial solution

{\hat{V}}_{1} (\hat{r}) = {\hat{V}}_{1 c}

⁠. However, we wish to set

{\hat{H}}_{1} = {\hat{V}}_{1} = 0

at

\hat{r} = 0

⁠, in order to ensure that

\hat{r} = 0

is indeed the magnetic axis. Hence, we are forced to conclude that

{\hat{V}}_{1} (\hat{r}) = 0

⁠, which accounts for the absence of this term in our analysis. Once we have determined

g_{2} (\hat{r})

⁠, as well as the

{\hat{H}}_{n} (\hat{r})

and the

{\hat{V}}_{n} (\hat{r})

functions, we can integrate Eq. to give

f_{3} (\hat{r}) = \frac{1}{f_{1} (\hat{r})} \int_{0}^{\hat{r}} F_{3} (\hat{r}') d \hat{r}' .

(41)

Note that

f_{3} (\hat{r}) = - f_{1 c} ({\hat{H}}_{2 c}^{2} + {\hat{V}}_{2 c}^{2}) {\hat{r}}^{2}

at small

\hat{r}

⁠. Furthermore, we have conveniently set

{\hat{f}}_{3} = 0

at

\hat{r} = 0

⁠, which ensures that

q (0) = [1 + ϵ^{2} ({\hat{H}}_{2 c}^{2} + {\hat{V}}_{2 c}^{2})] / f_{1 c}

⁠.

I. Near-vacuum region

We expect

p_{2} (\hat{r}) = 0

(42)

in the vacuum region immediately surrounding the plasma, which corresponds to

1 < \hat{r} ≪ ϵ^{- 1}

⁠. We also expect this so-called near-vacuum region to be current-free, which implies that

I_{t} (\hat{r}) = I_{t} (1)

and

I_{p} (\hat{r}) = I_{p} (1)

for

\hat{r} > 1

⁠. It follows from Eq. that

g_{2} (\hat{r}) = g_{2 a}

(43)

in the near-vacuum region, where

g_{2 a} = g_{2} (1)

is a constant determined from Eq. . Likewise, Eq. requires that

\frac{d}{d \hat{r}} (\hat{r} f ⟨ | \hat{\nabla} r |^{2} ⟩) = 0

(44)

in the region

\hat{r} > 1

⁠. Making use of the orderings introduced in Sec. II G, we deduce that

f_{1} (\hat{r}) = f_{1 a}

(45)

to lowest order, where

f_{1 a} = f_{1} (1)

is a constant that can be determined from the given

f_{1} (\hat{r})

profile, and

f_{3}^{'} (\hat{r}) = - f_{1 a} \frac{d {⟨ | \hat{\nabla} \hat{r} |^{2} ⟩}_{2}}{d \hat{r}}

(46)

to next order. Here,

{⟨ | \hat{\nabla} \hat{r} |^{2} ⟩}_{2} = \frac{3 {\hat{r}}^{2}}{4} - {\hat{H}}_{1} + \frac{1}{2} \sum_{n > 0} [{\hat{H}}_{n}^{' 2} + (n^{2} - 1) \frac{{\hat{H}}_{n}^{2}}{{\hat{r}}^{2}}] + \frac{1}{2} \sum_{n > 1} [{\hat{V}}_{n}^{' 2} + (n^{2} - 1) \frac{{\hat{V}}_{n}^{2}}{{\hat{r}}^{2}}]

(47)

is the component of

⟨ | \hat{\nabla} \hat{r} |^{2} ⟩

that is second order in ϵ. [See Eq. .]

The toroidal and poloidal current profiles in the plasma can be written

I_{t} (r) = ϵ^{2} (B_{0} R_{0} / μ_{0}) {\hat{I}}_{t} (\hat{r})

and

I_{p} (r) = ϵ^{2} (B_{0} R_{0} / μ_{0}) {\hat{I}}_{p} (\hat{r})

⁠, respectively. It follows from Eqs. and that the total normalized toroidal and poloidal plasma currents flowing in the plasma are

{\hat{I}}_{t a} = 2 π {[f_{1} (1 + ϵ^{2} {⟨ | \hat{\nabla} \hat{r} |^{2} ⟩}_{2}) + ϵ^{2} f_{3}]}_{\hat{r} = 1},

(48)

{\hat{I}}_{p a} = 2 π g_{2 a},

(49)

respectively.

In this paper, we shall assume that the equilibrium plasma current density is continuous across the plasma/vacuum interface, as is generally the case in a realistic tokamak plasma (assuming that the plasma current profile has been given sufficient time to diffuse across the plasma). Given that the current density is zero in the vacuum region, it follows that the current density must be zero just inside the plasma boundary. In other words, we require that

{\hat{I}}_{t}^{'} (1) = {\hat{I}}_{p}^{'} (1) = 0

⁠. Thus, making use of Eqs. and , we deduce that

g_{2}^{'} (1) = 0

⁠,

f_{1}^{'} (1) = 0,

(50)

and also that Eq. must be satisfied at

\hat{r} = 1

⁠. Finally, Eq. yields

p_{2}^{'} (1) = 0.

(51)

Note that Eq. is automatically satisfied in the near-vacuum region. Equations , combined with Eq. and , yield

{\hat{H}}_{1}^{″} = \frac{{\hat{H}}_{1}^{'}}{\hat{r}} - 1,

(52)

{\hat{H}}_{n}^{″} = \frac{{\hat{H}}_{n}^{'}}{\hat{r}} + (n^{2} - 1) \frac{{\hat{H}}_{n}}{{\hat{r}}^{2}} for n > 1,

(53)

{\hat{V}}_{n}^{″} = \frac{{\hat{V}}_{n}^{'}}{\hat{r}} + (n^{2} - 1) \frac{{\hat{V}}_{n}}{{\hat{r}}^{2}} for n > 1,

(54)

which can be solved to give¹

{\hat{H}}_{1} (\hat{r}) = {\hat{H}}_{1 a} - \frac{{\hat{r}}^{2}}{2} ln \hat{r} + \frac{1}{4} (2 {\hat{H}}_{1 a}^{'} + 1) ({\hat{r}}^{2} - 1),

(55)

{\hat{H}}_{n} (\hat{r}) = {\hat{H}}_{n a -} {\hat{r}}^{n + 1} - {\hat{H}}_{n a +} {\hat{r}}^{1 - n} for n > 1,

(56)

{\hat{V}}_{n} (\hat{r}) = {\hat{V}}_{n a -} {\hat{r}}^{n + 1} - {\hat{V}}_{n a +} {\hat{r}}^{1 - n} for n > 1,

(57)

in the near-vacuum region. Here,

{\hat{H}}_{n a -} = \frac{{\hat{H}}_{n a}^{'} + (n - 1) {\hat{H}}_{n a}}{2 n},

(58)

{\hat{H}}_{n a +} = \frac{{\hat{H}}_{n a}^{'} - (n + 1) {\hat{H}}_{n a}}{2 n},

(59)

{\hat{V}}_{n a -} = \frac{{\hat{V}}_{n a}^{'} + (n - 1) {\hat{V}}_{n a}}{2 n},

(60)

{\hat{V}}_{n a +} = \frac{{\hat{V}}_{n a}^{'} - (n + 1) {\hat{V}}_{n a}}{2 n} .

(61)

Moreover,

{\hat{H}}_{n a} = {\hat{H}}_{n} (1), {\hat{H}}_{n a}^{'} = {\hat{H}}_{n}^{'} (1), {\hat{V}}_{n a} = {\hat{V}}_{n} (1)

⁠, and

{\hat{V}}_{n a}^{'} = {\hat{V}}_{n}^{'} (1)

are constants determined from the solutions of Eqs. in the plasma interior, subject to the boundary conditions . It follows from Eqs. , , , and that

\begin{matrix} f_{3}^{'} (\hat{r}) = - \frac{f_{1 a}}{\hat{r}} (\frac{3 {\hat{r}}^{2}}{2} - 2 \hat{r} {\hat{H}}_{1}^{'} + \sum_{n > 0} [{\hat{H}}_{n}^{' 2} + 2 (n^{2} - 1) \frac{{\hat{H}}_{n}^{'} {\hat{H}}_{n}}{\hat{r}} - (n^{2} - 1) \frac{{\hat{H}}_{n}^{2}}{{\hat{r}}^{2}}] + \sum_{n > 1} [{\hat{V}}_{n}^{' 2} + 2 (n^{2} - 1) \frac{{\hat{V}}_{n}^{'} {\hat{V}}_{n}}{\hat{r}} - (n^{2} - 1) \frac{{\hat{V}}_{n}^{2}}{{\hat{r}}^{2}}]) . \end{matrix}

(62)

However, Eqs. and can be combined to give exactly the same relation, which confirms that zero toroidal current flows in the near-vacuum region, and also that the toroidal current density is zero at the plasma boundary. Equations and can be integrated to give

\begin{matrix} f_{3} (\hat{r}) = f_{3 a} - \frac{f_{1 a}}{2} [(1 - 2 {\hat{H}}_{1 a}^{'}) {\hat{r}}^{2} ln \hat{r} + {\hat{r}}^{2} {ln}^{2} \hat{r} + (1 - {\hat{H}}_{1 a}^{'} + {\hat{H}}_{1 a}^{' 2}) ({\hat{r}}^{2} - 1) + \sum_{n > 1} 2 n (n + 1) {\hat{H}}_{n a -}^{2} ({\hat{r}}^{2 n} - 1) + \sum_{n > 1} 2 n (n - 1) {\hat{H}}_{n a +}^{2} ({\hat{r}}^{- 2 n} - 1) \\ + \sum_{n > 1} 2 n (n + 1) {\hat{V}}_{n a -}^{2} ({\hat{r}}^{2 n} - 1) + \sum_{n > 1} 2 n (n - 1) {\hat{V}}_{n a +}^{2} ({\hat{r}}^{- 2 n} - 1)] \end{matrix}

(63)

in the near-vacuum region. Here,

f_{3 a} = f_{3} (1)

is a constant determined from Eq. .

Let

Ψ_{p} (r) = ϵ^{2} B_{0} R_{0}^{2} {\hat{Ψ}}_{p} (\hat{r})

⁠. It follows from Eqs. , , , and that

\begin{matrix} {\hat{Ψ}}_{p} (\hat{r}) = F_{1 a} + C_{1} + f_{1 a} ln \hat{r} + ϵ^{2} (F_{3 a} + C_{3}) + ϵ^{2} f_{3 a} ln \hat{r} \\ - \frac{ϵ^{2} f_{1 a}}{2} [- {1 - {\hat{H}}_{1 a}^{'} + {\hat{H}}_{1 a}^{' 2} + \sum_{n > 1} 2 n (n + 1) {\hat{H}}_{n a -}^{2} + \sum_{n > 1} 2 n (n - 1) {\hat{H}}_{n a +}^{2} + \sum_{n > 1} 2 n (n + 1) {\hat{V}}_{n a -}^{2} + \sum_{n > 1} 2 n (n - 1) {\hat{V}}_{n a +}^{2}} ln \hat{r} - {\hat{H}}_{1 a}^{'} {\hat{r}}^{2} ln \hat{r} + \frac{1}{2} {\hat{r}}^{2} {ln}^{2} \hat{r} \\ + \frac{1}{2} (1 + {\hat{H}}_{1 a}^{' 2}) ({\hat{r}}^{2} - 1) + \sum_{n > 1} (n + 1) {\hat{H}}_{n a -}^{2} ({\hat{r}}^{2 n} - 1) + \sum_{n > 1} (n - 1) {\hat{H}}_{n a +}^{2} (1 - {\hat{r}}^{- 2 n}) + \sum_{n > 1} (n + 1) {\hat{V}}_{n a -}^{2} ({\hat{r}}^{2 n} - 1) + \sum_{n > 1} (n - 1) {\hat{V}}_{n a +}^{2} (1 - {\hat{r}}^{- 2 n})] \end{matrix}

(64)

in the near-vacuum region, where

F_{1 a} = \int_{0}^{1} \frac{f_{1} (\hat{r})}{\hat{r}} d \hat{r},

(65)

F_{3 a} = \int_{0}^{1} \frac{f_{3} (\hat{r})}{\hat{r}} d \hat{r} .

(66)

Here, C₁ and C₃ are arbitrary constants.

III. MATCHING TO VACUUM SOLUTION

A. Toroidal coordinates

The r, θ, $ϕ$ coordinate system becomes singular in the vacuum region at large $\hat{r}$ [i.e., $\hat{r} \sim O (ϵ^{- 1})$ ]. Thus, in order to find a global vacuum solution, we need to match our previous near-vacuum solution to a vacuum solution calculated using a coordinate system that is nonsingular throughout the vacuum region.

An appropriate choice of a nonsingular vacuum coordinates is orthogonal toroidal coordinates,^1,7 μ, η,

ϕ

⁠, which are defined such that¹⁹

\hat{R} = \frac{sinh μ}{cosh μ - cos η},

(67)

\hat{Z} = \frac{sin η}{cosh μ - cos η} .

(68)

Here,

μ (\hat{R}, \hat{Z}) \to 0

corresponds to either

\hat{R} \to 0

or

{({\hat{R}}^{2} + {\hat{Z}}^{2})}^{1 / 2} \to \infty

(i.e., an approach to the toroidal symmetry axis or to infinity), whereas

μ (\hat{R}, \hat{Z}) \to \infty

corresponds to

(\hat{R}, \hat{Z}) \to (1, 0)

(i.e., an approach to the magnetic axis). Furthermore,

η (\hat{R}, \hat{Z})

is an angular variable in the poloidal plane.

B. Vacuum solution

According to Eqs. , , , and , the equilibrium magnetic field in the vacuum region can be written

B = \nabla ϕ \times \nabla Ψ_{p} + I_{a} \nabla ϕ,

(69)

where

I_{a} = B_{0} B_{0} (1 + ϵ^{2} g_{2 a})

is a constant. It follows that

\nabla \times B = R^{2} \nabla \cdot (\frac{\nabla Ψ_{p}}{R^{2}}) \nabla ϕ .

(70)

Hence, the condition that must be satisfied in order to ensure that no current flows in the vacuum is

\nabla \cdot (\nabla Ψ_{p} / R^{2}) = 0

⁠, which implies that

\hat{\nabla} \cdot (\frac{\hat{\nabla} {\hat{Ψ}}_{p}}{{\hat{R}}^{2}}) = 0.

(71)

When expressed in the toroidal coordinate system, the previous equation yields¹^,

\frac{\partial}{\partial μ} (\frac{cosh μ - cos η}{sinh μ} \frac{\partial {\hat{Ψ}}_{p}}{\partial μ}) + \frac{\partial}{\partial η} (\frac{cosh μ - cos η}{sinh μ} \frac{\partial {\hat{Ψ}}_{p}}{\partial η}) = 0.

(72)

If we write

{\hat{Ψ}}_{p} (μ, η) = \sum_{n \geq 0} \frac{sinh μ}{{(cosh μ - cos η)}^{1 / 2}} χ_{n} (cosh μ) e^{i n η},

(73)

where n is an integer (as must be the case to render the normalized poloidal magnetic flux,

{\hat{Ψ}}_{p}

⁠, single-valued in η), then we obtain

(1 - z^{2}) \frac{d^{2} χ_{n}}{d z^{2}} - 2 z \frac{d χ_{n}}{d z} + (n^{2} - \frac{1}{4} - \frac{1}{1 - z^{2}}) χ_{n} = 0,

(74)

where

z = cosh μ

⁠. The previous equation can be recognized as an associated Legendre function of degree

n - 1 / 2

and order 1.²⁰ The independent solutions of this equation are the so-called toroidal functions

P_{n - 1 / 2}^{1} (z)

and

Q_{n - 1 / 2}^{1} (z)

⁠.²¹ Note that

P_{n - 1 / 2}^{1} (z)

is well behaved in the limit

z \to 1

⁠, whereas

Q_{n - 1 / 2}^{1} (z)

is badly behaved. Conversely,

P_{n - 1 / 2}^{1} (z)

is badly behaved in the limit

z \to \infty

⁠, whereas

Q_{n - 1 / 2}^{1} (z)

is well behaved.

According to the previous analysis, the most general expression for the normalized poloidal magnetic flux in the vacuum region is¹

\begin{matrix} {\hat{Ψ}}_{p} (z, η) = {(\frac{z^{2} - 1}{z - cos η})}^{1 / 2} {\sum_{n \geq 0} [p_{c n} P_{n - 1 / 2}^{1} (z) + q_{c n} Q_{n - 1 / 2}^{1} (z)] cos (n η) \\ + \sum_{n > 0} [p_{s n} P_{n - 1 / 2}^{1} (z) + q_{s n} Q_{n - 1 / 2}^{1} (z)] sin (n η)}, \end{matrix}

(75)

where n is a non-negative integer. Here, the

p_{c n}, q_{c n}, p_{s n}

⁠, and

q_{s n}

are arbitrary coefficients. Moreover, it is clear from the asymptotic behaviors of the

P_{n - 1 / 2}^{1} (z)

and the

Q_{n - 1 / 2}^{1} (z)

functions that the terms in the previous series that involve the

p_{c n}

and

p_{s n}

represent the magnetic field generated in the vacuum region by currents flowing in the plasma, whereas the terms involving the

q_{c n}

and

q_{s n}

represent the magnetic field generated in the vacuum region by currents flowing in external magnetic field-coils. We now need to match expression to expression . Note that the latter expression is only valid in the near-vacuum region,

1 \leq \hat{r} ≪ ϵ_{a}^{- 1}

⁠, whereas the former expression is valid throughout the whole vacuum region.

C. Transformation of coordinates

In the immediate vicinity of the plasma boundary, Eqs. , , , , and can be combined to give

ω = π - η + ϵ F_{η} (\hat{r}, η) + O (ϵ^{2}),

(76)

where

F_{η} (\hat{r}, η) = - \frac{\hat{r}}{2} sin η + \sum_{n > 0} cos (n π) \frac{{\hat{H}}_{n}}{\hat{r}} sin (n η) + \sum_{n > 1} cos (n π) \frac{{\hat{V}}_{n}}{\hat{r}} cos (n η) .

(77)

Furthermore,

z = \frac{ξ (\hat{r}, η)}{ϵ \hat{r}},

(78)

where

\begin{matrix} ξ (\hat{r}, η) = 1 + ϵ [\frac{\hat{r}}{2} cos η + \sum_{n > 0} cos (n π) \frac{{\hat{H}}_{n}}{\hat{r}} cos (n η) - \sum_{n > 1} cos (n π) \frac{{\hat{V}}_{n}}{\hat{r}} sin (n η)] \\ + ϵ^{2} (\frac{5 {\hat{r}}^{2}}{16} - \frac{{\hat{H}}_{1}}{4} + \frac{3}{4} \sum_{n > 0} \frac{{\hat{H}}_{n}^{2}}{{\hat{r}}^{2}} + \frac{3}{4} \sum_{n > 1} \frac{{\hat{V}}_{n}^{2}}{{\hat{r}}^{2}}) . \end{matrix}

(79)

Note that

\hat{r} \sim O (1)

in the near-vacuum region, which implies that

z ≫ 1

⁠.

D. Toroidal functions

In the limit

z ≫ 1

⁠, the toroidal functions

P_{n - 1 / 2}^{1} (z)

and

Q_{n - 1 / 2}^{1} (z)

have the following asymptotic behaviors:^22,23

\frac{{(z^{2} - 1)}^{1 / 2} P_{- 1 / 2}^{1} (z)}{z^{1 / 2}} = \frac{1}{\sqrt{2} π} [2 - ln (8 z) - \frac{1 - ln (8 z)}{16 z^{2}} + O (z^{- 4})],

(80)

\frac{{(z^{2} - 1)}^{1 / 2} P_{n - 1 / 2}^{1} (z)}{z^{1 / 2}} = \frac{(n - 1)! 2^{n} z^{n}}{\sqrt{2 π} Γ (n - 1 / 2)} + O (z^{n - 2}) for n > 0,

(81)

\frac{{(z^{2} - 1)}^{1 / 2} Q_{- 1 / 2}^{1} (z)}{z^{1 / 2}} = \frac{\sqrt{π} Γ (3 / 2)}{\sqrt{2}} [1 - \frac{1}{16 z^{2}} + O (z^{- 4})],

(82)

\frac{{(z^{2} - 1)}^{1 / 2} Q_{n - 1 / 2}^{1} (z)}{z^{1 / 2}} = \frac{\sqrt{π} Γ (n + 3 / 2) z^{- n}}{2^{n + 1 / 2} n!} + O (z^{- n - 2}) for n > 0.

(83)

Note that there is a factor

\sqrt{- 1}

difference between the definition of the

P_{n - 1 / 2}^{1} (z)

used in this paper and that employed in Ref. 22. Moreover,

Γ (z)

is a Gamma function.²⁴

E. Poloidal magnetic flux in near-vacuum region

Let

p_{c 0} = \frac{π}{\sqrt{2}} {\hat{p}}_{c 0},

(84)

p_{c n} = \frac{\sqrt{2 π} Γ (n - 1 / 2) ϵ^{1 + n}}{(n - 1)! 2^{n}} {\hat{p}}_{c n} for n > 0,

(85)

p_{s n} = \frac{\sqrt{2 π} Γ (n - 1 / 2) ϵ^{1 + n}}{(n - 1)! 2^{n}} {\hat{p}}_{s n} for n > 0,

(86)

q_{c 0} = \frac{\sqrt{2}}{\sqrt{π} Γ (3 / 2)} {\hat{q}}_{0},

(87)

q_{c n} = \frac{2^{n + 1 / 2} n! ϵ^{1 - n}}{\sqrt{π} Γ (n + 3 / 2)} {\hat{q}}_{c n} for n > 0,

(88)

q_{s n} = \frac{2^{n + 1 / 2} n! ϵ^{1 - n}}{\sqrt{π} Γ (n + 3 / 2)} {\hat{q}}_{s n} for n > 0.

(89)

Here, the

{\hat{p}}_{c n}, {\hat{p}}_{s n}, {\hat{q}}_{c n}

⁠, and

{\hat{q}}_{s n}

coefficients are assumed to be

O (1)

⁠.

It follows from Eqs. , , and that

\begin{matrix} {\hat{Ψ}}_{p} (\hat{r}, η) = {(1 - \frac{ϵ \hat{r}}{ξ} cos η)}^{- 1 / 2} {{\hat{p}}_{c 0} [1 - \frac{1}{2} ln (\frac{8}{ϵ}) + \frac{1}{2} ln \hat{r} - \frac{1}{2} ln ξ] \\ - \frac{{\hat{p}}_{c 0}}{32} [1 - ln (\frac{8}{ϵ}) + ln \hat{r} - ln ξ] {(\frac{ϵ \hat{r}}{ξ})}^{2} \\ + {\hat{q}}_{c 0} [1 - \frac{1}{16} {(\frac{ϵ \hat{r}}{ξ})}^{2}] + ϵ \sum_{n > 0} {\hat{p}}_{c n} {(\frac{ξ}{\hat{r}})}^{n} cos (n η) \end{matrix}

\begin{matrix} + ϵ \sum_{n > 0} {\hat{q}}_{c n} {(\frac{\hat{r}}{ξ})}^{n} cos (n η) + ϵ \sum_{n > 0} {\hat{p}}_{s n} {(\frac{ξ}{\hat{r}})}^{n} sin (n η) + ϵ \sum_{n > 0} {\hat{q}}_{s n} {(\frac{\hat{r}}{ξ})}^{n} sin (n η)} \end{matrix}

(90)

in the near-vacuum region.

Making use of Eq. , the most general expression for the normalized poloidal magnetic flux in the near-vacuum region becomes

{\hat{Ψ}}_{p} (\hat{r}, η) = {\hat{Ψ}}_{p 0} (\hat{r}) + ϵ \sum_{n > 0} {\hat{Ψ}}_{p 1 cn} (\hat{r}) cos (n η) + ϵ \sum_{n > 0} {\hat{Ψ}}_{p 1 sn} (\hat{r}) sin (n η) + ϵ^{2} {\hat{Ψ}}_{p 2} (\hat{r}),

(91)

where

{\hat{Ψ}}_{p 0} (\hat{r}) = {\hat{p}}_{c 00} [1 - ln (\frac{8}{ϵ}) + \frac{ln \hat{r}}{2}] + {\hat{q}}_{c 00},

(92)

{\hat{Ψ}}_{pc 1} (\hat{r}) = \frac{{\hat{p}}_{c 00}}{4} [\hat{r} + \frac{2 {\hat{H}}_{1}}{\hat{r}} - ln (\frac{8}{ϵ}) \hat{r} + \hat{r} ln \hat{r}] + \frac{{\hat{p}}_{c 1}}{\hat{r}} + (\frac{{\hat{q}}_{c 00}}{2} + {\hat{q}}_{c 1}) \hat{r},

(93)

{\hat{Ψ}}_{pcn} (\hat{r}) = - {\hat{p}}_{c 00} cos (n π) \frac{{\hat{H}}_{n}}{2 \hat{r}} + {\hat{p}}_{c n} {\hat{r}}^{- n} + {\hat{q}}_{c n} {\hat{r}}^{n} for n > 1,

(94)

{\hat{Ψ}}_{psn} (\hat{r}) = {\hat{p}}_{c 00} cos (n π) \frac{{\hat{V}}_{n}}{2 \hat{r}} + {\hat{p}}_{s n} {\hat{r}}^{- n} + {\hat{q}}_{s n} {\hat{r}}^{n} for n > 0,

(95)

\begin{matrix} {\hat{Ψ}}_{ps 2} (\hat{r}) = {\hat{p}}_{c 02} [1 - ln (\frac{8}{ϵ}) + \frac{ln \hat{r}}{2}] + {\hat{q}}_{c 02} \\ + {\hat{p}}_{c 00} (- \frac{5 {\hat{r}}^{2}}{32} - [\frac{3}{8} + ln (\frac{8}{ϵ})] {\hat{H}}_{1} + \frac{{\hat{H}}_{1}}{8} ln \hat{r} - \sum_{n > 0} \frac{{\hat{H}}_{n}^{2}}{4 {\hat{r}}^{2}} - \sum_{n > 1} \frac{{\hat{V}}_{n}^{2}}{4 {\hat{r}}^{2}}) + (\frac{{\hat{q}}_{c 00}}{2} + {\hat{q}}_{c 1}) \frac{{\hat{H}}_{1}}{2} + \frac{{\hat{p}}_{c 1}}{2} + \sum_{n > 0} {\hat{p}}_{c n} cos (n π) \frac{n {\hat{H}}_{n}}{2} {\hat{r}}^{- 1 - n} - \sum_{n > 1} {\hat{q}}_{c n} cos (n π) \frac{n {\hat{H}}_{n}}{2} {\hat{r}}^{- 1 + n} - \sum_{n > 1} {\hat{p}}_{s n} cos (n π) \frac{n {\hat{V}}_{n}}{2} {\hat{r}}^{- 1 - n} + \sum_{n > 1} {\hat{q}}_{s n} cos (n π) \frac{n {\hat{V}}_{n}}{2} {\hat{r}}^{- 1 + n} . \end{matrix}

(96)

Here, we have written

{\hat{p}}_{c 0} = {\hat{p}}_{c 00} + ϵ^{2} {\hat{p}}_{c 02},

(97)

{\hat{q}}_{c 0} = {\hat{q}}_{c 00} + ϵ^{2} {\hat{q}}_{c 02},

(98)

where

{\hat{p}}_{c 00}, {\hat{p}}_{c 02}, {\hat{q}}_{c 00}

⁠, and

{\hat{q}}_{c 02}

are all

O (1)

⁠. Note that we have evaluated

{\hat{Ψ}}_{p} (\hat{r}, η)

up to

O (ϵ)

⁠, and retained

O (ϵ^{2})

contributions to those terms that are independent of η.

F. Comment

We now have two expressions for the poloidal magnetic flux in the near-vacuum region. The first, (64), is the extension of the Grad–Shafranov solution inside the plasma into the near-vacuum region. The second, (91), is the near-vacuum limit of the general vacuum solution (75). We require these two solutions to be identically equal to one another throughout the whole of the near-vacuum region. At first sight, it might seem that the problem is over-determined,²⁵ because there are far more distinct functional forms that require matching than there are free parameters in the problem. Fortunately, as described in GJW, if the matching is performed in systematic manner then exact matching is achieved, showing that the apparent over-determined nature of the problem is illusory. In fact, the achievement of exact matching is a very powerful internal self-consistency check on the whole calculation, because any algebraic error in the various terms that makeup the two solutions would lead to residual unmatched terms.

G. Matching

We must now match expression (91) to expression (64), assuming that the ${\hat{H}}_{n} (\hat{r})$ and ${\hat{V}}_{n} (\hat{r})$ functions takes the forms specified in Eqs. (55)–(57) in the near-vacuum region.

To zeroth order in ϵ, the matching yields

F_{1 a} + C_{1} + f_{1 a} ln \hat{r} = {\hat{p}}_{c 00} [1 - \frac{1}{2} ln (\frac{8}{ϵ}) + \frac{1}{2} ln \hat{r}] + {\hat{q}}_{c 00} .

(99)

Separately equating the coefficients of

ln \hat{r}

and

{\hat{r}}^{0}

⁠, we obtain¹

{\hat{p}}_{c 00} = 2 f_{1 a},

(100)

{\hat{q}}_{c 00} = F_{1 a} + C_{1} - 2 f_{1 a} [1 - \frac{1}{2} ln (\frac{8}{ϵ})] .

(101)

To first order in ϵ, matching the coefficients of

cos η

yields

[{\hat{p}}_{c 1} - f_{1 a} (\frac{1}{4} + \frac{{\hat{H}}_{1 a}^{'}}{2} - {\hat{H}}_{1 a})] {\hat{r}}^{- 1} + ({\hat{q}}_{c 1} + \frac{{\hat{q}}_{c 00}}{2} - \frac{f_{1 a}}{2} [ln (\frac{8}{ϵ}) - \frac{3}{2} - {\hat{H}}_{1 a}^{'}]) \hat{r} = 0.

(102)

Separately equating the coefficients of

{\hat{r}}^{- 1}

and

{\hat{r}}^{1}

⁠, we obtain¹

{\hat{p}}_{c 1} = f_{1 a} (\frac{1}{4} + \frac{{\hat{H}}_{1 a}^{'}}{2} - {\hat{H}}_{1 a}),

(103)

{\hat{q}}_{c 1} + \frac{{\hat{q}}_{c 00}}{2} = \frac{f_{1 a}}{2} [ln (\frac{8}{ϵ}) - \frac{3}{2} - {\hat{H}}_{1 a}^{'}] .

(104)

To first order in ϵ, matching the coefficients of

cos (n η)

⁠, where n > 1, yields

[{\hat{p}}_{c n} + f_{1 a} cos (n π) {\hat{H}}_{n a +}] {\hat{r}}^{- n} + [{\hat{q}}_{c n} - f_{1 a} cos (n π) {\hat{H}}_{n a -}] {\hat{r}}^{n} = 0.

(105)

Separately equating the coefficients of

{\hat{r}}^{- n}

and

{\hat{r}}^{n}

⁠, we obtain

{\hat{p}}_{c n} = - f_{1 a} cos (n π) {\hat{H}}_{n a +},

(106)

{\hat{q}}_{c n} = f_{1 a} cos (n π) {\hat{H}}_{n a -} .

(107)

To first order in ϵ, matching the coefficients of

sin η

yields

{\hat{p}}_{s 1} {\hat{r}}^{- 1} + {\hat{q}}_{c 1} \hat{r} = 0,

(108)

where we have made use of the fact that

{\hat{V}}_{1} = 0

⁠. Separately equating the coefficients of

{\hat{r}}^{- 1}

and

\hat{r}

⁠, we obtain

{\hat{p}}_{s 1} = 0,

(109)

{\hat{q}}_{s 1} = 0.

(110)

To first order in ϵ, matching the coefficients of

sin (n η)

⁠, where n > 1, yields

[{\hat{p}}_{s n} - f_{1 a} cos (n π) {\hat{V}}_{n a +}] {\hat{r}}^{- n} + [{\hat{q}}_{c n} + f_{1 a} cos (n π) {\hat{V}}_{n a -}] {\hat{r}}^{n} = 0.

(111)

Separately equating the coefficients of

{\hat{r}}^{- n}

and

{\hat{r}}^{n}

⁠, we obtain

{\hat{p}}_{s n} = f_{1 a} cos (n π) {\hat{V}}_{n a +},

(112)

{\hat{q}}_{s n} = - f_{1 a} cos (n π) {\hat{V}}_{n a -} .

(113)

Finally, to second order in ϵ, the matching yields

{\hat{p}}_{c 02} [1 - ln (\frac{8}{ϵ}) + \frac{ln \hat{r}}{2}] + {\hat{q}}_{c 02} = α_{c 02} ln \hat{r} + β_{c 02},

(114)

where

\begin{matrix} α_{c 02} = F_{3 a} + C_{3} + f_{1 a} (\frac{3}{16} - \frac{{\hat{H}}_{1 a}^{'}}{4} + \frac{{\hat{H}}_{1 a}^{'} {\hat{H}}_{1 a}}{2} + \frac{{\hat{H}}_{1 a}}{4} \\ - \sum_{n > 1} [(\frac{n - 1}{2}) {\hat{H}}_{n a +}^{2} + {\hat{H}}_{n a +} {\hat{H}}_{n a -} - (\frac{n + 1}{2}) \frac{{\hat{H}}_{n a -}^{2}}{2}] \\ - \sum_{n > 1} [(\frac{n - 1}{2}) {\hat{V}}_{n a +}^{2} + {\hat{V}}_{n a +} {\hat{V}}_{n a -} - (\frac{n + 1}{2}) \frac{{\hat{V}}_{n a -}^{2}}{2}]), \end{matrix}

(115)

\begin{matrix} β_{c 02} = f_{3 a} + f_{1 a} (\frac{5}{8} - \frac{{\hat{H}}_{1 a}^{'}}{4} + \frac{{\hat{H}}_{1 a}^{' 2}}{2} - \frac{{\hat{H}}_{1 a}}{2} \\ + \sum_{n > 1} [n (n - 1) {\hat{H}}_{n a +}^{2} + n (n + 1) {\hat{H}}_{n a -}^{2}] \\ + \sum_{n > 1} [n (n - 1) {\hat{V}}_{n a +}^{2} + n (n + 1) {\hat{V}}_{n a -}^{2}]) . \end{matrix}

(116)

Note that many terms proportional to

{\hat{r}}^{2} {ln}^{2} \hat{r}

and

{\hat{r}}^{2} ln \hat{r}

and

{\hat{r}}^{- 2 n}

and

{\hat{r}}^{2 n}

have canceled one another in expression on the basis of the previous matching (see Sec. III F). Separately equating the coefficients of

ln \hat{r}

and

{\hat{r}}^{0}

⁠, we obtain

{\hat{p}}_{c 02} = 2 α_{c 02},

(117)

{\hat{q}}_{c 02} = β_{c 02} - 2 α_{c 02} [1 - \frac{1}{2} ln (\frac{8}{ϵ})] .

(118)

It is convenient to set the value of the arbitrary constant C₃ so as to ensure that

{\hat{q}}_{c 02} = 0

⁠, which implies that

α_{c 02} = β_{c 02} / (2 [1 - \frac{1}{2} ln (\frac{8}{ϵ})]),

(119)

F_{3 a} + C_{3} = β_{c 02} / (2 [1 - \frac{1}{2} ln (\frac{8}{ϵ})]) - γ_{c 02},

(120)

where

\begin{matrix} γ_{c 02} = f_{1 a} (\frac{3}{16} - \frac{{\hat{H}}_{1 a}^{'}}{4} + \frac{{\hat{H}}_{1 a}^{'} {\hat{H}}_{1 a}}{2} + \frac{{\hat{H}}_{1 a}}{4} \\ - \sum_{n > 1} [(\frac{n - 1}{2}) {\hat{H}}_{n a +}^{2} + {\hat{H}}_{n a +} {\hat{H}}_{n a -} - (\frac{n + 1}{2}) \frac{{\hat{H}}_{n a -}^{2}}{2}] \\ - \sum_{n > 1} [(\frac{n - 1}{2}) {\hat{V}}_{n a +}^{2} + {\hat{V}}_{n a +} {\hat{V}}_{n a -} - (\frac{n + 1}{2}) \frac{{\hat{V}}_{n a -}^{2}}{2}]) . \end{matrix}

(121)

We have now completed the matching process.

H. Externally generated magnetic field

Let

\hat{R} = 1 + ϵ ρ cos Ω,

(122)

\hat{Z} = ϵ ρ sin Ω,

(123)

where

ρ (R, Z)

measures radial distance from the magnetic axis, and

Ω (R, Z)

is a geometric poloidal angle that is distinct from both θ and ω. Using similar analysis to that employed in Secs. III C–III E, we deduce that the normalized poloidal magnetic flux generated by the currents flowing in the external magnetic field-coils takes the form

\begin{matrix} {\hat{Ψ}}_{x} (ρ, Ω) = {\hat{q}}_{c 00} + ϵ (\frac{{\hat{q}}_{c 00}}{2} + {\hat{q}}_{c 1}) ρ cos Ω + ϵ \sum_{n > 1} {\hat{q}}_{c n} ρ^{n} cos (n Ω) \\ + ϵ \sum_{n > 1} {\hat{q}}_{s n} ρ^{n} sin (n Ω) + ϵ^{2} {\hat{q}}_{c 02} + ϵ^{2} (\frac{{\hat{q}}_{c 00}}{2} + {\hat{q}}_{c 1}) \frac{ρ^{2}}{4} \end{matrix}

(124)

in the immediate vicinity of the magnetic axis. If we write

\hat{x} = ρ cos Ω

(125)

\hat{z} = ρ sin Ω,

(126)

then we obtain

\begin{matrix} {\hat{Ψ}}_{x} (\hat{x}, \hat{z}) = {\hat{q}}_{c 00} + ϵ \frac{f_{1 a}}{2} [ln (\frac{8}{ϵ}) - \frac{3}{2} + {\hat{H}}_{1 a}^{'}] \hat{x} + ϵ f_{1 a} \sum_{n > 1} cos (n π) {\hat{H}}_{n a -} P_{c n} (\hat{x}, \hat{z}) - ϵ f_{1 a} \sum_{n > 1} cos (n π) {\hat{V}}_{n a -} P_{s n} (\hat{x}, \hat{z}) + ϵ^{2} {\hat{q}}_{c 02} + ϵ^{2} \frac{f_{1 a}}{8} [ln (\frac{8}{ϵ}) - \frac{3}{2} + {\hat{H}}_{1 a}^{'}] ({\hat{x}}^{2} + {\hat{z}}^{2}), \end{matrix}

(127)

where²⁶

P_{c n} (\hat{x}, \hat{z}) = {\hat{x}}^{n} - (\begin{matrix} n \\ 2 \end{matrix}) {\hat{x}}^{n - 2} {\hat{z}}^{2} + (\begin{matrix} n \\ 4 \end{matrix}) {\hat{x}}^{n - 4} {\hat{z}}^{4} + \dots,

(128)

P_{s n} (\hat{x}, \hat{z}) = n {\hat{x}}^{n - 1} \hat{z} - (\begin{matrix} n \\ 3 \end{matrix}) {\hat{x}}^{n - 3} {\hat{z}}^{3} + (\begin{matrix} n \\ 5 \end{matrix}) {\hat{x}}^{n - 5} {\hat{z}}^{5} + \dots,

(129)

and use has been made of the results of Sec. III G. Here,

\hat{x}

and

\hat{z}

measure horizontal and vertical distance from the magnetic axis, respectively. In particular, if

{\hat{H}}_{n a -} = {\hat{V}}_{n a -} = 0

for n > 3 then expression reduces to

\begin{matrix} {\hat{Ψ}}_{x} (\hat{x}, \hat{z}) = {\hat{q}}_{c 00} + ϵ^{2} {\hat{q}}_{c 02} + ϵ \frac{f_{1 a}}{2} [ln (\frac{8}{ϵ}) - \frac{3}{2} + {\hat{H}}_{1 a}^{'}] \hat{x} + ϵ^{2} \frac{f_{1 a}}{8} [ln (\frac{8}{ϵ}) - \frac{3}{2} + {\hat{H}}_{1 a}^{'}] ({\hat{x}}^{2} + {\hat{z}}^{2}) \\ + ϵ f_{1 a} {\hat{H}}_{2 a -} ({\hat{x}}^{2} - {\hat{z}}^{2}) - ϵ f_{1 a} {\hat{V}}_{2 a -} 2 \hat{x} \hat{z} \\ - ϵ f_{1 a} {\hat{H}}_{3 a -} ({\hat{x}}^{3} - 3 \hat{x} {\hat{z}}^{2}) + ϵ f_{1 a} {\hat{V}}_{3 a -} (3 {\hat{x}}^{2} \hat{z} - {\hat{z}}^{3}) . \end{matrix}

(130)

Let

B_{x} = B_{0} {\hat{B}}_{x}

be the magnetic field generated by the external field-coils in the immediate vicinity of the magnetic axis. It is easily seen from Eqs. , , , , , and that

{\hat{B}}_{x} = (\frac{\partial {\hat{Ψ}}_{x}}{\partial \hat{z}} e_{R} - \frac{\partial {\hat{Ψ}}_{x}}{\partial \hat{x}} e_{Z}) {(1 + ϵ \hat{x})}^{- 1},

(131)

where

e_{R} = \nabla R

and

e_{Z} = \nabla Z

⁠. Hence,

\begin{matrix} {\hat{B}}_{x} = - ϵ \frac{f_{1 a}}{2} [ln (\frac{8}{ϵ}) - \frac{3}{2} + {\hat{H}}_{1 a}^{'}] e_{Z} + ϵ^{2} \frac{f_{1 a}}{4} [ln (\frac{8}{ϵ}) - \frac{3}{2} + {\hat{H}}_{1 a}^{'}] (\hat{z} e_{R} + \hat{x} e_{Z}) - ϵ 2 f_{1 a} {\hat{H}}_{2 a -} (\hat{z} e_{R} + \hat{x} e_{Z}) - ϵ 2 f_{1 a} {\hat{V}}_{2 a -} (\hat{x} e_{R} - \hat{z} e_{Z}) + ϵ 3 f_{1 a} {\hat{H}}_{3 a -} [2 \hat{x} \hat{z} e_{R} + ({\hat{x}}^{2} - {\hat{z}}^{2}) e_{Z}] + ϵ 3 f_{1 a} {\hat{V}}_{3 a -} [({\hat{x}}^{2} - {\hat{z}}^{2}) e_{R} - 2 \hat{x} \hat{z} e_{Z}] . \end{matrix}

(132)

Here, the first term on the right-hand side of the previous equation is the vertical field that counteracts the hoop stress in the plasma torus.^1,13 The second term is the higher order quadrupole correction to the vertical field. The third and the fourth terms are quadrupole magnetic fields that control the up-down symmetric and asymmetric ellipticities of the plasma equilibrium magnetic flux-surfaces, respectively. Note that these fields are mutually orthogonal. In fact, one can be transformed into the other via a rotation in the poloidal plane through

90 °

⁠. Finally, the fifth and the sixth terms are sextupole magnetic fields that control the up-down symmetric and asymmetric triangularities of the equilibrium magnetic flux-surfaces, respectively. As before, the fields are mutually orthogonal, and one can be transformed into the other via a rotation in the poloidal plane through

90 °

⁠.

IV. EXTERNAL MAGNETIC FIELD-COILS

A. External magnetic field-coil

Suppose that the set of magnetic field-coils that generate the external poloidal magnetic field, discussed in Sec. III H, that supports the plasma equilibrium are discrete and located a finite distance from the plasma.

Let us crudely model a given magnetic field-coil as made up of a number of toroidal loops, of zero cross-section in the poloidal plane, that we shall term strands. See Fig. 1. Suppose that a given strand carries a toroidal current I_x, and is located at coordinates (R_x, Z_x) in the cylindrical system. The strand is located at coordinates

(μ_{x}, η_{x})

in the toroidal system, where²⁷

z_{x} \equiv cosh μ_{x} = \frac{d_{1}^{2} + d_{2}^{2}}{2 d_{1} d_{2}},

(133)

cos η_{x} = \frac{d_{1}^{2} + d_{2}^{2} - 4}{2 d_{1} d_{2}},

(134)

d_{1}^{2} = {({\hat{R}}_{x} + 1)}^{2} + {\hat{Z}}_{x}^{2},

(135)

d_{2}^{2} = {({\hat{R}}_{x} - 1)}^{2} + {\hat{Z}}_{x}^{2} .

(136)

Here,

{\hat{R}}_{x} = R_{x} / R_{0}

and

{\hat{Z}}_{x} = Z_{x} / R_{0}

⁠. Obviously, the strand is assumed to lie outside the plasma.

FIG. 1.

View large Download slide

The blue curves show the internal magnetic flux-surfaces of an axisymmetric tokamak plasma equilibrium characterized by $ϵ = 0.25, q_{c} = 1.0, ν = 3.0, p_{c} = 0.25, μ = 2.1$ ⁠, and $Δ_{ped} = 0$ ⁠. The red curve denotes the plasma boundary, and the red dot corresponds to the magnetic axis. Here, OH is the Ohmic heating coil, and P₁ through P₆ are the six poloidal field-coils. Each green dot corresponds to a constituent strand (of zero poloidal cross-sectional area) of a given field-coil that carries the same toroidal current. The relative weights of the net toroidal currents flowing in OH and P₁ to P₆ are 0.80, 0.02, 0.01, −0.01, −0.02, 0.65, and 0.75, respectively.

B. Ampère–Maxwell equation

Combining Eq. with the Ampère–Maxwell equation, we obtain

R \nabla \cdot (\frac{\nabla Ψ_{x}}{R^{2}}) = \frac{μ_{0} I_{x} δ (μ - μ_{x}) δ (η - η_{x})}{h^{2}},

(137)

where

Ψ_{x} (R, Z)

is the poloidal magnetic flux generated by the strand, and

h = \frac{R_{0}}{cosh μ - cos η}

(138)

is the poloidal scale-factor for the toroidal coordinate system.²⁷ Now,²⁸

\nabla \cdot (\frac{\nabla Ψ_{x}}{R^{2}}) = \frac{1}{R h^{2}} [\frac{\partial}{\partial μ} (\frac{1}{R} \frac{\partial Ψ_{x}}{\partial μ}) + \frac{\partial}{\partial η} (\frac{1}{R} \frac{\partial Ψ_{x}}{\partial η})] .

(139)

Hence, if we write

Ψ_{x} = ϵ^{2} B_{0} R_{0}^{2} {\hat{Ψ}}_{x}

and

I_{x} = ϵ^{2} (B_{0} R_{0} / μ_{0}) {\hat{I}}_{x}

then we obtain

\frac{\partial}{\partial z} [(z - cos η) \frac{\partial {\hat{Ψ}}_{x}}{\partial z}] + \frac{1}{z^{2} - 1} \frac{\partial}{\partial η} [(z - cos η) \frac{\partial {\hat{Ψ}}_{x}}{\partial η}] = {\hat{I}}_{x} δ (z - z_{x}) δ (η - η_{x}),

(140)

where

z = cosh μ

⁠.

C. Determination of matching coefficients

Let

{\hat{Ψ}}_{x} (z, η) = {(\frac{z^{2} - 1}{z - cos η})}^{1 / 2} χ_{x} (z, η) .

(141)

Equation transforms to give

\begin{matrix} (1 - z^{2}) \frac{\partial^{2} χ_{x}}{\partial z^{2}} - 2 z \frac{\partial χ_{x}}{\partial z} - (\frac{1}{4} + \frac{1}{1 - z^{2}}) χ_{x} - \frac{\partial^{2} χ_{x}}{\partial η^{2}} \\ = - {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} δ (z - z_{x}) δ (η - η_{x}) . \end{matrix}

(142)

Now,²⁹^,

δ (η - η_{x}) = \frac{1}{2 π} + \frac{1}{π} \sum_{n = 1, \infty} cos [n (η - η_{x})] .

(143)

Let

χ_{x} (z, η) = \sum_{n = 0, \infty} χ_{x n} (z) cos [n (η - η_{x})] .

(144)

The previous three equations yield

\begin{matrix} (1 - z^{2}) \frac{d^{2} χ_{x n}}{d z^{2}} - 2 z \frac{d χ_{x n}}{d z} + (n^{2} - \frac{1}{4} - \frac{1}{1 - z^{2}}) χ_{x n} \\ = - \frac{ζ_{n}}{π} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} δ (z - z_{x}), \end{matrix}

(145)

where

ζ_{n} = {\begin{array}{l} 1 / 2, & n = 0, \\ 1, & n > 0. \end{array}

(146)

The solution of Eq. that is well-behaved as

z \to 1

and

z \to \infty

⁠, and is continuous at z = z_x, is^20,21

χ_{x n} (z) = A_{n} {\begin{matrix} Q_{n - 1 / 2}^{1} (z_{x}) P_{n - 1 / 2}^{1} (z), & z < z_{x}, \\ Q_{n - 1 / 2}^{1} (z) P_{n - 1 / 2}^{1} (z_{x}), & z > z_{x} . \end{matrix}

(147)

Integration of Eq. across z = z_x reveals that

A_{n} (1 - z_{x}^{2}) {[P_{n - 1 / 2} \frac{d Q_{n - 1 / 2}^{1}}{d z} - Q_{n - 1 / 2}^{1} \frac{d P_{n - 1 / 2}^{1}}{d z}]}_{z = z_{x}} = - \frac{ζ_{n}}{π} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} .

(148)

However,³⁰

P_{n - 1 / 2} \frac{d Q_{n - 1 / 2}^{1}}{d z} - Q_{n - 1 / 2}^{1} \frac{d P_{n - 1 / 2}^{1}}{d z} = \frac{n^{2} - 1 / 4}{1 - z^{2}},

(149)

which yields

A_{n} = - \frac{ζ_{n}}{π (n^{2} - 1 / 4)} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} .

(150)

The previous analysis implies that, in the region interior to the strand,

\begin{matrix} {\hat{Ψ}}_{x} (z, η) = - {\hat{I}}_{x} {(\frac{z^{2} - 1}{z - cos η})}^{1 / 2} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{x}) \\ \times \sum_{n = 0, \infty} \frac{ζ_{n}}{π (n^{2} - 1 / 4)} Q_{n - 1 / 2}^{1} (z) [cos (n η_{x}) cos (n η) + sin (n η_{x}) sin (n η)] . \end{matrix}

(151)

A comparison with expression reveals that the contributions of the strand to the matching coefficients that characterize the external magnetic field are

q_{c n} = - \frac{ζ_{n}}{π (n^{2} - 1 / 4)} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{x}) cos (n η_{x}),

(152)

q_{s n} = - \frac{ζ_{n}}{π (n^{2} - 1 / 4)} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{x}) sin (n η_{x}) .

(153)

Finally, making use of the normalizations , we obtain

{\hat{q}}_{c 0} = \frac{2 Γ (3 / 2)}{\sqrt{2 π}} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} P_{- 1 / 2}^{1} (z_{x}),

(154)

{\hat{q}}_{c n} = - \frac{Γ (n + 3 / 2) ϵ^{n - 1}}{\sqrt{2 π} 2^{n} n! (n^{2} - 1 / 4)} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{x}) cos (n η_{x}) for n > 0,

(155)

{\hat{q}}_{s n} = - \frac{Γ (n + 3 / 2) ϵ^{n - 1}}{\sqrt{2 π} 2^{n} n! (n^{2} - 1 / 4)} {\hat{I}}_{x} {(\frac{z_{x}^{2} - 1}{z_{x} - cos η_{x}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{x}) sin (n η_{x}) for n > 0.

(156)

Suppose, finally, that the external magnetic field-coils are made up of K strands. Let the kth strand carry the toroidal current

I_{k} = ϵ^{2} (B_{0} R_{0} / μ_{0}) {\hat{I}}_{k}

and be located at z = z_k and

η = η_{k}

⁠. It follows that

{\hat{q}}_{c 0} = \frac{2 Γ (3 / 2)}{\sqrt{2 π}} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{- 1 / 2}^{1} (z_{k}),

(157)

{\hat{q}}_{c n} = - \frac{Γ (n + 3 / 2) ϵ^{n - 1}}{\sqrt{2 π} 2^{n} n! (n^{2} - 1 / 4)} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{k}) cos (n η_{k}) for n > 0,

(158)

{\hat{q}}_{s n} = - \frac{Γ (n + 3 / 2) ϵ^{n - 1}}{\sqrt{2 π} 2^{n} n! (n^{2} - 1 / 4)} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{k}) sin (n η_{k}) for n > 0.

(159)

D. Matching

We must now match our coil-generated poloidal magnetic flux to our previous plasma solution.

Given that we adjusted the value of the arbitrary constant C₃ to ensure that

{\hat{q}}_{c 02} = 0

⁠, in Sec. III G, it follows from Eq. that

{\hat{q}}_{c 00} = {\hat{q}}_{c 0}

⁠. Equations and yield

F_{1 a} + C_{1} = \frac{1}{\sqrt{2}} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{- 1 / 2}^{1} (z_{k}) + 2 f_{1 a} [1 - \frac{1}{2} ln (\frac{8}{ϵ})] .

(160)

The previous equation implies that, adopting the convention that the poloidal magnetic flux is zero at infinity, the normalized poloidal magnetic flux at the plasma boundary (including both plasma and vacuum contributions) is

\begin{matrix} {\hat{Ψ}}_{p a} = \frac{1}{\sqrt{2}} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{- 1 / 2}^{1} (z_{k}) + 2 f_{1 a} [1 - \frac{1}{2} ln (\frac{8}{ϵ})] \\ + ϵ^{2} [β_{c 02} / (2 [1 - \frac{1}{2} ln (\frac{8}{ϵ})]) - γ_{c 02}] . \end{matrix}

(161)

[See Eqs. , , and .]

Equations and give

{\hat{q}}_{c 1} + \frac{{\hat{q}}_{c 00}}{2} = \frac{1}{\sqrt{8}} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} [P_{- 1 / 2}^{1} (z_{k}) - P_{1 / 2}^{1} (z_{k}) cos η_{k}] .

(162)

It follows from Eq. that

\sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} [P_{- 1 / 2}^{1} (z_{k}) - P_{1 / 2}^{1} (z_{k}) cos η_{k}] = \sqrt{2} f_{1 a} [ln (\frac{8}{ϵ}) - \frac{3}{2} - {\hat{H}}_{1 a}^{'}] .

(163)

The previous equation is the condition that must be satisfied in order for the plasma to be in a state of horizontal force balance.

Equations and yield

\sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{1 / 2}^{1} (z_{k}) sin η_{k} = 0.

(164)

The previous equation is the condition that must be satisfied in order for the plasma to be in a state of vertical force balance.

Equations and give

{\hat{H}}_{n a -} = - \frac{Γ (n + 3 / 2) ϵ^{n - 1}}{\sqrt{2 π} 2^{n} n! (n^{2} - 1 / 4)} \frac{cos (n π)}{f_{1 a}} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{k}) cos (n η_{k})

(165)

for n > 1. The previous equation determines the up-down symmetric shaping of the plasma equilibrium.

Finally, Eqs. and yield

{\hat{V}}_{n a -} = \frac{Γ (n + 3 / 2) ϵ^{n - 1}}{\sqrt{2 π} 2^{n} n! (n^{2} - 1 / 4)} \frac{cos (n π)}{f_{1 a}} \sum_{k = 1, K} {\hat{I}}_{k} {(\frac{z_{k}^{2} - 1}{z_{k} - cos η_{k}})}^{1 / 2} P_{n - 1 / 2}^{1} (z_{k}) sin (n η_{k})

(166)

for n > 1. The previous equation determines the up-down asymmetric shaping of the plasma equilibrium. We have now completed the matching process.

V. PEDESTAL PHYSICS

A. Introduction

Tokamak plasmas possessing magnetic nulls on the plasma/vacuum interface generally operate in the so-called H-mode regime.^3,31 In this regime, a transport barrier forms just inside the plasma boundary, giving rise to a narrow edge region in which the (negative) equilibrium pressure gradient is much larger than that in the plasma interior. This high-gradient region is known as the pedestal. The non-inductive bootstrap current³² is a parallel plasma current that is driven by pressure gradients.³³ Hence, as a consequence of the bootstrap current, as well as the strong pressure gradients present in the pedestal, H-mode tokamak plasmas generally contain a localized spike in the equilibrium parallel plasma current in the pedestal. We wish to incorporate these features of realistic tokamak plasma equilibria into our model.

B. Pedestal model

In order to incorporate the pedestal in the plasma pressure into our model, we can modify our normalized pressure profile as follows:

p_{2} (\hat{r}) = p_{2 original} (\hat{r}) + Δ p_{2} (\hat{r}),

(167)

where

Δ p_{2} (\hat{r}) = \frac{Δ_{ped}}{2 {\hat{δ}}_{ped}} \int_{\hat{r}}^{1} (\frac{1 - {\hat{r}}^{' 2}}{1 - {\hat{r}}_{ped}^{2}}) \frac{d \hat{r}'}{{cosh}^{2} [(\hat{r}' - {\hat{r}}_{ped}) / {\hat{δ}}_{ped}]} .

(168)

Here,

Δ_{ped}, {\hat{δ}}_{ped}

⁠, and

{\hat{r}}_{ped}

are the normalized height, width, and radius of the pedestal, respectively. Note that it is still the case that

p_{2}^{'} (1) = 0

⁠, assuming that

p_{2 original}^{'} (1) = 0

⁠, which ensures that the poloidal current density at the plasma boundary remains zero. [See Eq. .] In the limit that

{\hat{δ}}_{ped} ≪ 1

⁠, the previous two expressions reduce to

p_{2} (\hat{r}) ≃ p_{2 original} (\hat{r}) + \frac{Δ_{ped}}{2} [1 - tanh (\frac{\hat{r} - {\hat{r}}_{ped}}{{\hat{δ}}_{ped}})] .

(169)

Clearly, we have effectively introduced into the plasma pressure profile a negative step of height

Δ_{ped}

at the pedestal radius.

C. Bootstrap current model

The flux-surface averaged parallel plasma current density can be written,⁷

⟨ {(\frac{R}{R_{0}})}^{2} j \cdot B ⟩ = \frac{B_{0}^{2}}{μ_{0}} \frac{g^{2}}{r} \frac{d}{d r} (\frac{r f}{g} ⟨ | \nabla r |^{2} ⟩) .

(170)

To lowest order in our expansion, the previous equation yields

j_{∥} ≃ \frac{B_{0} ϵ}{μ_{0} a} \frac{1}{\hat{r}} \frac{d f_{1}}{d \hat{r}} .

(171)

Now, a fairly primitive expression for the bootstrap current density is³⁴

j_{bs} = - 1.46 \frac{{(ϵ \hat{r})}^{1 / 2}}{B_{θ}} \frac{d p}{d r} .

(172)

To lowest order in our expansion, the previous equation gives

j_{bs} ≃ - 1.46 \frac{B_{0} ϵ^{3 / 2}}{μ_{0} a} \frac{{\hat{r}}^{3 / 2}}{f_{1}} \frac{d p_{2}}{d \hat{r}} .

(173)

A comparison between Eqs. and reveals that we can incorporate the change in the bootstrap current profile associated with the pedestal in the plasma pressure into our model by modifying our

f_{1} (\hat{r})

profile as follows:

\frac{1}{2} f_{1}^{2} (\hat{r}) = \frac{1}{2} f_{1 original}^{2} (\hat{r}) - 1.46 ϵ^{1 / 2} \int_{0}^{\hat{r}} {\hat{r}}^{' 5 / 2} \frac{d Δ p_{2}}{d \hat{r}'} d \hat{r}' .

(174)

Note that it is still the case that

f_{1}^{'} (1) = 0

⁠, assuming that

f_{1 original}^{'} (1) = 0

⁠, and given that

Δ p_{2}^{'} (1) = 0

⁠, which ensures that the toroidal current density at the plasma boundary remains zero [See Eq. ].

VI. EXAMPLE CALCULATIONS

A. Introduction

The example calculations described in this sections are made for illustrative purpose only, and are not intended to be particularly realistic.

B. Model profiles

Our model lowest-order toroidal current and pressure profiles are³⁵

f_{1} (\hat{r}) = \frac{1}{ν q_{c}} [1 - {(1 - {\hat{r}}^{2})}^{ν}],

(175)

p_{2} (\hat{r}) = p_{c} {(1 - {\hat{r}}^{2})}^{μ},

(176)

respectively. Here,

ϵ^{2} (B_{0}^{2} / μ_{0}) p_{c}

is the central plasma pressure, and q_c the safety-factor on the magnetic axis. Note that the model profiles satisfy the constraints and provided that

ν > 2

and

μ > 2

⁠.

C. Coil set

Our model coil set consists of a 29-strand Ohmic heating coil located within the plasma torus, plus six 5-strand poloidal field-coils arrayed around the outer side of the plasma torus. See Fig. 1. Each constituent strand of a given field-coil carries the same toroidal current.

D. Example 1

Figure 1 shows an example axisymmetric tokamak equilibrium characterized by $ϵ = 0.25, q_{c} = 1.0, ν = 3.0, p_{c} = 0.25, μ = 2.1$ ⁠, and $Δ_{ped} = 0$ ⁠.

The equilibrium pictured in Fig. 1 is calculated as follows. First, the total toroidal currents flowing in the Ohmic heating coil, OH, and the six poloidal field-coils, P₁ to P₆, are given the relative weights 0.80, 0.02, 0.01, −0.01, −0.02, 0.65, and 0.75, respectively. Next, Eq. (163) is used to determine the total toroidal current that must flow in the whole coil-set in order for the plasma to be in horizontal force balance. Next, the coil-set is shifted vertically (which is equivalent to shifting the plasma vertically) until the vertical force-balance constraint (164) is satisfied [adjusting the total toroidal current flowing in the coil-set such that Eq. (163) is always satisfied]. Finally, the shape of the plasma boundary is deduced from Eqs. (165) and (166).

In this particular example, the normalized vertical field is ${\hat{B}}_{v} = - 0.168$ ⁠, the normalized total toroidal plasma current is ${\hat{I}}_{t a} = 2.11$ ⁠, the normalized total poloidal plasma current is ${\hat{I}}_{p a} = 6.98 \times 10^{- 2}$ ⁠, and the value of the safety-factor at the plasma boundary is $q_{2} (1) = 4.91$ ⁠.

Figure 2 shows the values of the shaping functions at the plasma boundary, ${\hat{H}}_{n a}$ and ${\hat{V}}_{n a}$ ⁠, for the example equilibrium pictured in Fig. 1. It can be seen that ${\hat{H}}_{1 a} < 0$ ⁠, which corresponds to the usual outward horizontal shift of the plasma axis with respect to the plasma boundary. Moreover, ${\hat{H}}_{2 a} > 0$ ⁠, which corresponds to a vertical elongation of the plasma boundary. It can be seen that the ${\hat{V}}_{n a}$ are finite, but generally much smaller in magnitude than the ${\hat{H}}_{n a}$ ⁠, indicating a modest up-down asymmetry of the plasma (due to the slight up-down asymmetry of the currents flowing in the poloidal field-coils). Finally, it is clear that the ${\hat{H}}_{n a}$ and the ${\hat{V}}_{n a}$ are negligible for n > 10. This implies that, despite the fairly strong flux-surface shaping evident in Fig. 1, the plasma equilibrium can be accurately described in terms of a few shaping harmonics.

FIG. 2.

View large Download slide

Values of the shaping functions at the plasma boundary for the example tokamak plasma equilibrium shown in Fig. 1. Here, the $H_{n a}$ control the up-down symmetric shaping of magnetic flux-surfaces, whereas the $V_{n a}$ control the up-down asymmetric shaping. Moreover, n is a poloidal mode number.

Finally, Fig. 3 shows various characteristic profiles for the example equilibrium pictured in Fig. 1. It can be seen that the corrected safety-factor profile, $q_{2} (\hat{r})$ ⁠, climbs more steeply in the edge regions of the plasma than the lowest-order profile, $q_{0} (\hat{r})$ ⁠. Moreover, the normalized toroidal and poloidal plasma current profiles, ${\hat{I}}_{t} (\hat{r})$ and ${\hat{I}}_{p} (\hat{r})$ ⁠, respectively, have zero gradients at the plasma boundary, implying that the plasma current density is continuous across the plasma/vacuum interface.

FIG. 3.

View large Download slide

Various characteristic profiles for the example tokamak plasma equilibrium shown in Fig. 1. Here, p₂ is the normalized plasma pressure profile, g₂ controls the variation of the toroidal magnetic field-strength, f₁ and f₃ control the variation of the poloidal magnetic field-strength, q₀ and q₂ are the lowest-order and corrected safety-factor profiles, respectively, and ${\hat{I}}_{t}$ and ${\hat{I}}_{p}$ are the normalized toroidal and poloidal currents, respectively, flowing within flux-surfaces. Finally, $\hat{r}$ is a flux-surface label that is zero at the magnetic axis, and unity at the plasma boundary.

E. Example 2

In our second example calculation, we add a pedestal characterized by $Δ_{ped} = 0.05, {\hat{δ}}_{ped} = 0.025$ ⁠, and ${\hat{r}}_{ped} = 0.95$ to the plasma equilibrium shown in Fig. 1. The resulting equilibrium is illustrated in Figs. 4–6. The pedestal in the plasma pressure is clearly evident in the top left-hand panel of Fig. 6, whereas the spike in the bootstrap current can be seen in the bottom left-hand panel. Note that, as a consequence of the bootstrap current spike, the safety-factor profile becomes slightly non-monotonic in the pedestal. In this particular example, the normalized vertical field is ${\hat{B}}_{v} = - 0.199$ ⁠, the normalized total toroidal plasma current is ${\hat{I}}_{a} = 2.80$ ⁠, the normalized total poloidal plasma current is ${\hat{I}}_{p a} = - 2.19 \times 10^{- 2}$ ⁠, and the value of the safety-factor at the plasma boundary is $q_{2} (1) = 3.67$ ⁠.

FIG. 4.

Internal magnetic flux-surfaces of a axisymmetric tokamak plasma equilibrium characterized by ϵ=0.25, qc=1.0, ν=3.0, pc=0.25, μ=2.1, Δped=0.05, δ̂ped=0.025, and r̂ped=0.95. See Fig. 1 caption.

View large Download slide

Internal magnetic flux-surfaces of a axisymmetric tokamak plasma equilibrium characterized by $ϵ = 0.25, q_{c} = 1.0, ν = 3.0, p_{c} = 0.25, μ = 2.1$ ⁠, $Δ_{ped} = 0.05, {\hat{δ}}_{ped} = 0.025$ ⁠, and ${\hat{r}}_{ped} = 0.95$ ⁠. See Fig. 1 caption.

FIG. 5.

View large Download slide

Values of the shaping functions at the plasma boundary for the example tokamak equilibrium shown in Fig. 4. See Fig. 2 caption.

FIG. 6.

View large Download slide

Various characteristic profiles for the example tokamak equilibrium shown in Fig. 4. See Fig. 3 caption.

VII. SUMMARY AND DISCUSSION

Following Greene et al. (GJW)¹ and Connor et al.^4–7 we have solved the Grad–Shafranov equation for an axisymmetric tokamak plasma equilibrium via an expansion in the, supposedly small, inverse aspect-ratio of the plasma, ϵ. We have assumed that the displacements of equilibrium magnetic flux-surfaces due to plasma shaping are $O (ϵ)$ smaller than the minor radii of the surfaces, but have, otherwise, placed no restriction on the nature of the shaping. In particular, we allow for an infinite number of shaping harmonics, and also for a lack of up-down symmetry of the plasma. Following GJW, we have matched our solution of the Grad–Shafranov equation to a vacuum solution that extends to infinity, and consists of an expansion in toroidal functions. We have also calculated the external poloidal magnetic field generated by a finite set of discrete external poloidal magnetic field-coils, and incorporated that calculation into our toroidal function expansion. In this manner, we are able to directly relate the shape of a large aspect-ratio tokamak plasma to the currents flowing in the external poloidal field-coils. Finally, we have incorporated a pedestal in the plasma pressure, located in the outer regions of the plasma, and the associated spike in the bootstrap current, into our model.

The main value of our calculation lies in the fact that it can directly determine the metric elements of the plasma equilibrium via Eqs. (21)–(25). This implies that the calculation can be used as the basis for an aspect-ratio-expanded determination of the tearing mode stability of the plasma equilibrium, as described in Ref. 7, or an aspect-ratio-expanded calculation of the response of the equilibrium to an externally applied resonant magnetic perturbation (RMP), along the lines of Ref. 36. In future work, we intend to perform these calculations.

ACKNOWLEDGMENTS

This research was directly funded by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, under Contract No. DE-SC0021156. The calculations described this paper have been verified using the REDUCE computer algebra system.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Richard Fitzpatrick: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The digital data used in the figures in this paper can be obtained from the author upon reasonable request.

REFERENCES

1.

J. M.

Greene

,

J. L.

Johnson

, and

K. E.

Weimer

,

Phys. Fluids

14

,

671

(

1971

).

https://doi.org/10.1063/1.1693488

Google Scholar

Crossref

2.

V. D.

Shafranov

,

Sov. At. Energ.

13

,

1149

(

1963

).

https://doi.org/10.1007/BF01312317

Google Scholar

Crossref

3.

J. A.

Wesson

and

D. J.

Campbell

,

Tokamaks

, 4th ed. (

Oxford University Press

,

Oxford, UK

,

2011

).

Google Scholar

4.

J. W.

Connor

and

R. J.

Hastie

, “

The effect of shaped plasma cross sections on the ideal kink mode in a tokamak

,” Report No. CLM-M106 (

Culham Laboratory

,

Abingdon, UK

,

1985

).

Google Scholar

5.

J. W.

Connor

,

S. C.

Cowley

,

R. J.

Hastie

,

T. C.

Hender

,

A.

Hood

, and

T. J.

Martin

,

Phys. Fluids

31

,

577

(

1988

).

https://doi.org/10.1063/1.866840

Google Scholar

Crossref

6.

R.

Fitzpatrick

,

C. G.

Gimblett

, and

R. J.

Hastie

,

Plasma Phys. Controlled Fusion

34

,

161

(

1992

).

https://doi.org/10.1088/0741-3335/34/2/003

Google Scholar

Crossref

7.

R.

Fitzpatrick

,

R. J.

Hastie

,

T. J.

Martin

, and

C. M.

Roach

,

Nucl. Fusion

33

,

1553

(

1993

).

https://doi.org/10.1088/0029-5515/33/10/l11

Google Scholar

8.

G. T. A.

Huysmans

,

J. P.

Goedbloed

, and

W.

Kerner

,

Int. J. Mod. Phys. C

2

,

371

(

1991

).

https://doi.org/10.1142/S0129183191000512

Google Scholar

Crossref

9.

H.

Lütjens

,

A.

Bondeson

, and

O.

Sauter

,

Comput. Phys. Commun.

97

,

219

(

1996

).

https://doi.org/10.1016/0010-4655(96)00046-X

Google Scholar

Crossref

10.

L. S.

Solov'ev

,

Sov. Phys. JETP

26

,

400

(

1968

).

Google Scholar

11.

A. J.

Cerfon

and

J. P.

Freidberg

,

Phys. Plasmas

17

,

032502

(

2010

).

https://doi.org/10.1063/1.3328818

Google Scholar

Crossref

12.

L.

Guazzotto

and

J. P.

Freidberg

,

J. Plasma Phys.

87

,

905870305

(

2021

).

https://doi.org/10.1017/S0022377821000118

Google Scholar

Crossref

13.

J. F.

Freidberg

,

Ideal Magnetohydrodynamics

(

Plenum Press

,

New York

,

1987

).

Google Scholar

Crossref

14.

T.

Xu

and

R.

Fitzpatrick

,

Nucl. Fusion

59

,

064002

(

2019

).

https://doi.org/10.1088/1741-4326/ab1ce3

Google Scholar

Crossref

15.

J.

Liu

,

P.

Zhu

, and

H.

Li

,

Phys. Plasmas

29

,

084502

(

2022

).

https://doi.org/10.1063/5.0096256

Google Scholar

Crossref

16.

A. H.

Boozer

,

Rev. Mod. Phys.

76

,

1071

(

2005

).

https://doi.org/10.1103/RevModPhys.76.1071

Google Scholar

Crossref

17.

R.

Iacono

,

A.

Bondeson

,

F.

Troyon

, and

R.

Gruber

,

Phys. Fluids B

2

,

1794

(

1990

).

https://doi.org/10.1063/1.859451

Google Scholar

Crossref

18.

L.

Guazzotto

,

R.

Betti

,

J.

Manickam

, and

S.

Kaye

,

Phys. Plasmas

11

,

604

(

2004

).

https://doi.org/10.1063/1.1637918

Google Scholar

Crossref

19.

P. M.

Morse

and

H.

Feshbach

,

Methods of Theoretical Physics

(

McGraw-Hill

,

New York

,

1953

), p.

1301

.

Google Scholar

20.

M.

Abramowitz

and

I. A.

Stegun

,

Handbook of Mathematical Functions

(

Dover

,

New York

,

1964

), p.

332

.

Google Scholar

21.

M.

Abramowitz

and

I. A.

Stegun

,

Handbook of Mathematical Functions

(

Dover

,

New York

,

1964

), Sec. 8.11.

Google Scholar

22.

P. M.

Morse

and

H.

Feshbach

,

Methods of Theoretical Physics

(

McGraw-Hill

,

New York

,

1953

), pp.

1302

–

1329

.

Google Scholar

23.

A.

Erdélyi

,

W.

Magnus

,

F.

Oberhettinger

, and

F. C.

Tricomi

,

Higher Transcendental Functions

(

McGraw-Hill

,

New York

,

1953

), Vol.

1

, Sec. 1.7.3.

Google Scholar

24.

M.

Abramowitz

and

I. A.

Stegun

,

Handbook of Mathematical Functions

(

Dover

,

New York

,

1964

), Chap. 6.

Google Scholar

25.

J. P.

Goedbloed

,

R.

Keppens

, and

S.

Poedts

,

Magnetohydrodynamics of Laboratory and Astrophysical Plasmas

(

Cambridge University Press

.

Cambridge, UK

,

2019

), Chap. 16.

Google Scholar

Crossref

26.

I. S.

Gradshteyn

and

I. M.

Ryzhik

,

Table of Integrals, Series, and Products

(

Academic Press

,

San Diego

,

1980

), p.

27

.

Google Scholar

27.

G.

Arfken

,

Mathematical Methods for Physicists

(

Academic Press

,

New York

,

1968

), p.

112

.

Google Scholar

28.

P. M.

Morse

and

H.

Feshbach

,

Methods of Theoretical Physics

(

McGraw-Hill

,

New York

,

1953

), p.

36

.

Google Scholar

29.

J. D.

Jackson

,

Classical Electrodynamics

, 2nd ed. (

John Wiley & Sons

,

New York

,

1975

), p.

117

.

Google Scholar

30.

P. M.

Morse

and

H.

Feshbach

,

Methods of Theoretical Physics

(

McGraw-Hill

,

New York

,

1953

), p.

1330

.

Google Scholar

31.

F.

Wagner

,

G.

Becker

,

K.

Behringer

,

D.

Campbell

,

A.

Eberhagen

,

W.

Engelhardt

,

G.

Fussmann

,

O.

Gehre

,

J.

Gernhardt

,

G. V.

Gierke

et al,

Phys. Rev. Lett.

49

,

1408

(

1982

).

https://doi.org/10.1103/PhysRevLett.49.1408

Google Scholar

Crossref

32.

R. J.

Bickerton

,

J. W.

Connor

, and

J. B.

Taylor

,

Nat. Phys. Sci.

229

,

110

(

1971

).

https://doi.org/10.1038/physci229110a0

Google Scholar

Crossref

33.

R.

Fitzpatrick

,

Tearing Mode Dynamics in Tokamak Plasmas

(

IOP

,

Bristol, UK

,

2023

).

Google Scholar

Crossref

34.

M. N.

Rosenbluth

,

R. D.

Hazeltine

, and

F. L.

Hinton

,

Phys. Fluids

15

,

116

(

1972

).

https://doi.org/10.1063/1.1693728

Google Scholar

Crossref

35.

J. A.

Wesson

,

Nucl. Fusion

18

,

87

(

1978

).

https://doi.org/10.1088/0029-5515/18/1/010

Google Scholar

Crossref

36.

R.

Fitzpatrick

,

Phys. Plasmas

24

,

072506

(

2017

).

https://doi.org/10.1063/1.4990701

Google Scholar

Crossref

Inverse aspect-ratio expanded tokamak equilibria

I. INTRODUCTION

II. PLASMA SOLUTION

A. Coordinates

B. Equilibrium magnetic field and current

C. Equilibrium magnetic flux-surfaces

D. Grad–Shafranov equation

E. Normalization

F. Metric elements

G. Expansion of Grad–Shafranov equation

H. Plasma interior

I. Near-vacuum region

III. MATCHING TO VACUUM SOLUTION

A. Toroidal coordinates

B. Vacuum solution

C. Transformation of coordinates

D. Toroidal functions

E. Poloidal magnetic flux in near-vacuum region

F. Comment

G. Matching

H. Externally generated magnetic field

IV. EXTERNAL MAGNETIC FIELD-COILS

A. External magnetic field-coil

B. Ampère–Maxwell equation

C. Determination of matching coefficients

D. Matching

V. PEDESTAL PHYSICS

A. Introduction

B. Pedestal model

C. Bootstrap current model

VI. EXAMPLE CALCULATIONS

A. Introduction

B. Model profiles

C. Coil set

D. Example 1

E. Example 2

VII. SUMMARY AND DISCUSSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only