1 Introduction

A routine in tokamak operation is to reconstruct magnetic field under the constraints of MHD force balance and various measurements. This kind of task can be done by various codes, e.g., EFIT. I developed a similar code called heq (github.com/youjunhu/heq).

I also discuss the fixed-boundary equilibrium problem, where the boundary magnetic surface is given and one is asked to solve the force-balance equation within that boundray.

Another topic is how to construct spatial coordinates that conform to the shape of the magnetic surfaces, based on discrete numerical equilibrium data output by the equilibrium reconstructing codes.

Let us begin with the basics.

2 Axisymmetric magnetic field

Due to the divergence-free nature of magnetic field, i.e., ∇⋅ B = 0, magnetic field can be expressed as the curl of a vector field,

(1)

where A is called the vector potential of B. In cylindrical coordinates (R,ϕ,Z), the above expression is written

(2)

We consider axisymmetric magnetic field. The axial symmetry means that, when expressed in the cylindrical coordinate system (R,ϕ,Z), the components of B, namely B_R, B_Z, and B_ϕ, are all independent of ϕ. For this case, it can be proved that an axisymmetric vector potential A suffices for expressing the magnetic field, i.e., all the components of the vector potential A can also be taken independent of ϕ. Using this, Eq. (2) is written

(3)

In tokamak literature, direction is called the toroidal direction, and (R,Z) planes (i.e., ϕ = const planes) are called poloidal planes.

2.1 Poloidal magnetic field

Equation (3) indicates that the two poloidal components of B, namely B_R and B_Z, are determined by a single component of A, namely A_ϕ. This motivates us to define a function Ψ(R,Z):

(4)

Then Eq. (3) implies the poloidal components, B_R and B_Z, can be written as

(5)

(6)

(Note that it is the property of being axisymmetric and divergence-free that enables us to express the two components of B, namely B_R and B_Z, in terms of a single function Ψ(R,Z).) Furthermore, it is ready to prove that Ψ is constant along a magnetic field line, i.e. B ⋅∇Ψ = 0. [Proof:

The function Ψ is usually called the “poloidal flux function” in tokamak literature because Ψ can be related to the poloidal magnetic flux, as we will discuss in Sec. 2.7.

Using Eqs. (5) and (6), the poloidal magnetic field B_p is written as

2.2 Toroidal magnetic field

Next, let’s examine the toroidal component B_ϕ. Equation (3) indicates that B_ϕ involves both A_R and A_Z, which means that using the potential form does not enable useful simplification for B_ϕ. Therfore we will directly use B_ϕ. Define g ≡ RB_ϕ(R,Z) (the reason that we define this quantity will become clear when we discuss the forece balance), then the toroidal magnetic field is written

(9)

2.3 General form of axisymmetric magnetic field

Combining Eqs. (8) and (9), we obtain

2.4 Gauge transform of Ψ

Let us discuss the gauge transformation of the vector potential A in the axisymmetric case. As is well known, magnetic field remains the same under the following gauge transformation:

(11)

where f is an arbitrary scalar field. Here we require that ∇f be axisymmetric because, as mentioned above, an axisymmetric vector potential suffices for describing an axisymmetric magnetic field. In cylindrical coordinates, ∇f is given by

(12)

Since ∇f is axisymmetric, it follows that all the three components of ∇f are independent of ϕ, i.e., ∂²f∕∂R∂ϕ = 0, ∂²f∕∂Z∂ϕ = 0, and ∂²f∕∂ϕ² = 0, which implies that ∂f∕∂ϕ is independent of R, Z, and ϕ, i.e., ∂f∕∂ϕ  is actually a spatial constant. Using this, the ϕ component of the gauge transformation (11) is written

(14)

which indicates Ψ has the same gauge transformation as the electric potential, i.e., adding a constant. Note that the definition Ψ(R,Z) ≡ RA_ϕ does not imply Ψ(R = 0,Z) = 0 because A_ϕ can adopt 1∕R dependence under the gauge transformation (13). However, we are usually interested in the situation where A_ϕ is finite at R = 0 (and axisymmetry indicates that A_ϕ(R = 0,Z) is a spatial constant). In this case, Ψ is zero at R = 0. This is the case we encounter in equilibrium reconstruction, where magnetic flux measurements are used to constrain the numerically reconstructed Ψ.

2.5 Contours of Ψ in the poloidal plane

Because Ψ is constant along a magnetic field line and Ψ is independent of ϕ, it follows that the projection of a magnetic field line onto (R,Z) plane is a contour of Ψ. On the other hand,  are contours of Ψ on (R,Z) plane also the projections of magnetic field lines onto the plane? Yes, they are. [Proof. A contour of Ψ on (R,Z) plane satisfies

(15)

i.e.,

(16)

(17)

Using Eqs. (5) and (6), the above equation is written

(18)

i.e.,

(19)

which is the equation of the projection of a magnetic field line in (R,Z) plane. This proves that contours of Ψ are also the projections of magnetic field lines in (R,Z) plane.]

2.6 Magnetic surfaces

For axial symmetry system, magnetic surfaces can be defined in a trivial way. The axial symmetry of tokamak magnetic field allows us to introduce a surface of revolution that is generated by rotating the projection of a magnetic field line on (R,Z) plane around the axis of symmetry, Z axis. The unique property of this revolution surface is that no field line point-intersects it and a field line with one point on it will have the whole field line on it. This revolution surface can be called a magnetic surface. For instance, consider an arbitrary magnetic field line, whose projection on the poloidal plane is shown in Fig. 1. A magnetic surface is generated by rotating the projection line around the Z axis (some people disagree with this defintion since they only consider toroidal surfaces that are closed within a tokamak vessell as magnetic surfaces.)

The value of Ψ is constant on a magnetic surface (since Ψ is constant along a magnetic field line and Ψ is independent of ϕ). On the other hand, the value of Ψ is usually different on different magnetic surfaces. The above two facts enable Ψ to be used as labels of magnetic surfaces.

2.7 Relation of Ψ with the poloidal magnetic flux

Note that Ψ is defined by Ψ = RA_ϕ, which is just a component of the vector potential A, thereby having no obvious physical meaning. Next, we try to find the physical meaning of Ψ, i.e., some simple algebraic relation of Ψ with some quantity that can be measured in experiments.

In Fig. 2, there are two magnetic surfaces labeled, respectively, by Ψ = Ψ₁ and Ψ = Ψ₂. Using Gauss’s theorem in the toroidal volume between the two magnetic surface, we know that the poloidal magnetic flux through any toroidal ribbons between the two magnetic surfaces is equal to each other. Next, we calculate this poloidal magnetic flux. To make the calculation easy, we select a plane perpendicular to the Z axis, as is shown by the dash line in Fig. 1. In this case, only B_Z contribute to the poloidal magnetic flux, which is written (the positive direction of the plane is chosen in the direction of )

2.8 Measurement of poloidal magnetic flux

By measuring the voltage around a toroidal wire loop, we can obtain the time derivative of the poloidal magnetic flux though the loop and, after integrating over time, the poloidal magnetic flux itself.

Suppose that the loop is located at (R,Z) and denote the magnetic flux through the loop by Ψ_pol(R,Z,t) (only the poloidal magnetic field contribute to this flux since the loop is in the toroidal direction). Then Faraday’s law gives

(21)

where 𝜀 is the emf. If the loop is a coil with N turns, the induced voltage V in the coil is N times the emf 𝜀, i.e., V = N𝜀. Using this, Eq. (21) is written as

(22)

Integrating the above equation over time, we obtain

(23)

The starting time t = 0 can be chosen at a time when no plasma is present and thus the initial value, Ψ_pol(R,Z,0), is easy to obtain.

As discussed in Sec. 2.7 ,the relation between Ψ and Ψ_pol is given by

(24)

where Ψ(R = 0,Z,t), as discussed above, is zero. Then Eq. (24) is written as

(25)

which tells us how to obtain Ψ form the measured Ψ_pol(R,Z,t), which in trun is inferred from the measured voltage V , as is given by Eq. (23).

Note that the positive normal direction of the surface (where the magnetic flux Ψ_pol is defined) is chosen in the Z direction.

The above toroidal loop of wire used to measure poloidal magnetic flux is often called “flux loop” by tokamak operators. There are usually many flux loops (e.g. 35 on EAST[31], 41 on DIII-D) at different locations in the poloidal plane (see Fig. 3). The measured value of poloidal flux can be used as constraints in reconstructing magnetic configuration by, e.g. EFIT. With the feedback control (by adjusting the current in the poloidal field coils), these measurement, along with the poloidal field measurement by magnetic probes, can be used to control the shape of the last-closed-flux-surface (LCFS), the distance of the LCFS from the first wall, and the X-point location. This is often called iso-flux control, gap control, or X-point control.

2.9 Closed magnetic surfaces

In most part of a tokamak plasma, contours of Ψ on (R,Z) plane are closed curves. As discussed above, the contours of Ψ are the projection of magnetic lines on the poloidal plane. Closed contours of Ψ implies closed magnetic surfaces, as shown in Fig 4.

The innermost magnetic surface reduces to a curve, which is called magnetic axis (in Fig. 4, Ψ₀ labels the magnetic axis). Because the magnetic axis is the point of maximum/minimum of Ψ(R,Z), the value of ∇Ψ is zero at the magnetic axis. As a result, the poloidal component of the equilibrium magnetic field is zero on magnetic axis (refer to Eq. (8)), i.e., the magnetic field has only toroidal component there.

As is discussed in Sec. 2.7, the poloidal magnetic flux enclosed by a magnetic surface Ψ (the poloidal magnetic flux through the toroidal surfaces S₂) is given by

(26)

where Ψ₀ is the value of Ψ at the magnetic axis. Here the positive direction of the surface S₂ is defined to be in the clockwise direction when an observer looks along the direction of . In practice, we need to pay attention to the positive direction chosen (there can be a sign difference when choosing different positive directions).

Also note, the poloidal magnetic flux enclosed by a closed magnetic surface can have two different definitions, one of which is the the poloidal magnetic flux through the surface S₂ in Fig. 4, the other one is the poloidal magnetic flux through the central hole of the magnetic surface, i.e., the poloidal flux through S₁ in Fig. 4. In the latter case, as is discussed in Sec. 2.8, the poloidal magnetic flux is related to Ψ by

(27)

where the positive direction of the surface S₁ is in the clockwise direction.

Also note that, since the poloidal magnetic field can be written as B_p = ∇Ψ ×∇ϕ, the condition Ψ_LCFS − Ψ_axis > 0 means B_p points in the anticlockwise direction (viewed along ϕ direction), and Ψ_LCFS − Ψ_axis < 0 means B_p points in the clockwise direction.

Shape parameters of a closed magnetic surface

Let us introduce parameters characterizing the shape of a magnetic surface in the poloidal plane. The “midplane” is defined as the plane that passes through the magnetic axis and is perpendicular to the symmetric axis (Z axis). For a up-down symmetric (about the midplane) magnetic surface, its shape can be roughly characterized by four parameters, namely, the R coordinate of the innermost and outermost points in the midplane, R_in and R_out; the (R,Z) coordinators of the highest point of the magnetic surface, (R_top,Z_top). These four parameters are indicated in Fig. 5.

In terms of these four parameters, we can define the major radius of a magnetic surface

(28)

(which is the R coordinate of the geometric center of the magnetic surface), the minor radius of a magnetic surface

(29)

the triangularity of a magnetic surface

(30)

and, the ellipticity (elongation) of a magnetic surface

(31)

Usually, we specify the value of R₀, a, δ_t, and κ, instead of (R_in,R_out,R_top,Z_top), to characterize the shape of a magnetic surface. The value of the triangularity δ_t is usually positive in traditional tokamak operations, but negative triangularity is achievable and potentially useful, which is under active investigation.

Besides, using a and R₀, we can define another useful parameter 𝜀 ≡ a∕R₀, which is called the inverse aspect ratio.

The four shape parameters for the typical Last-Closed-Flux-Surface (LCFS) of EAST tokamak are: major radius R₀ = 1.85m (can reach 1.9m), minor radius a = 0.45m, ellipticity κ = 1.8 (can be in the range from 1.5 to 2.0), triangularity δ_t = 0.6 (can be in the range from 0.3 to 0.6). Note that the major radius R₀ of the LCFS is usually different from R_axis (the R coordinate of the magnetic axis). Usually we have R_axis > R₀ due to the so-called Shafranov shift.

2.10 Safety factor

Magnetic field lines on closed magnetic surfaces travel closed curves in a poloidal plane. For these field lines, we can define the safety factor q, which is the number of toroidal loops a magnetic field line travels when it makes one poloidal loop, i.e.

(32)

where △ϕ as the change of the toroidal angle when a magnetic field line travels a full poloidal loop. The safety factor can also be understood as the average pitch angle of a magnetic field line in (𝜃,ϕ) plane of a closed magnetic surface.

For open field line region (where a field line touches the wall before its poloidal projection can close itself), the “connection length” is often used to characterize the magnetic field.

Expression of safety factor in terms of magnetic field

The equation of magnetic field lines is given by

(33)

where dℓ_p is the line element along the direction of B_p on the poloidal plane. Equation (33) can be arranged in the form

(34)

which can be integrated over dℓ_p to give

(35)

where the line integration is along the poloidal magnetic field (the contour of Ψ on the poloidal plane). Using this, Eq. (32) is written

(36)

Expression of safety factor in terms of magnetic flux

The safety factor given by Eq. (36) is expressed in terms of the components of the magnetic field. The safety factor can also be expressed in terms of the magnetic flux. Define ΔΨ_p as the poloidal magnetic flux enclosed by two neighboring magnetic surface, then ΔΨ_p is given by

(37)

where Δx is the length of a line segment in the poloidal plane between the two magnetic surfaces, which is perpendicular to the first magnetic surface (so perpendicular to the B_p). Note that Δx, as well as R and B_p, generally depends on the poloidal location whereas ΔΨ_p is independent of the poloidal location.

Using Eq. (37), the poloidal magnetic field is written as

(38)

Substituting Eq. (38) into Eq. (36), we obtain

(39)

We know ΔΨ_p is a constant independent of the poloidal location, so ΔΨ_p can be taken outside the integration to give

(40)

It is ready to realise that the integral appearing in Eq. (40) is the toroidal magnetic flux enclosed by the two magnetic surfaces, ΔΨ_t. Using this, Eq. (40) is written as

(41)

Equation (41) indicates that the safety factor of a magnetic surface is equal to the differential of the toroidal magnetic flux with respect to the poloidal magnetic flux enclosed by the magnetic surface.

Rational surfaces vs. irrational surfaces

If the safety factor of a magnetic surface is a rational number, i.e., q = m∕n, where m and n are integers, then this magnetic surface is called a rational surface, otherwise an irrational surface. A field line on a rational surface with q = m∕n closes itself after it travels n poloidal loops. An example of a field line on a rational surface is shown in Fig. 6.

3 Non-axisymmetric magnetic perturbations

In the above, the magnetic field is assumed to be axisymmetric. With this assumption, the poloidal magnetic field (having two components) can be expressed in terms of a single component of the vector potential A, A_ϕ (specifically via Ψ ≡ A_ϕR). This kind of simplification can not be achieved if the axisymmetricity assumption is dropped, because other components of the vector potential (namely A_R and A_Z) will appear in the expression of the poloidal magnetic field. Let us re-examine Eq. (2) for a general magnetic perturbation:

(43)

(44)

(45)

where δΨ = RδA_ϕ. Therefore this kind of magnetic perturbation can still be written in the same form of the equilibrium poloidal magnetic field:

(46)

The above approximation is widely used in practice, e.g., tearing mode theory; turbulence simulation, where δA_ϕ is replaced by δA_∥. (Do we miss some magnetic perturbations that is important for plasma transport when using the above specific form?)

The total magnetic field is then written as

(47)

**check**Can the projection of the total magnetic field line in the poloidal plane can be traced by tracing the contour of Ψ + δΨ? No. The contours of Ψ + δΨ will not show island structures in the poloidal plane. To show the expected island structures, we need to subtract non-reconnecting poloidal magnetic field from the total poloidal field. I will discuss this in Sec. (), where I will show that contours of the so-called helical flux will give the expected island structures near the resonant surfaces.**can be worng, I have not checked this.

Next, we go back to discuss the 2D case (i.e., assuming axisymmetry).

4 Plasma current density in terms of Ψ and g

When the displacement current term is neglectable (the case we consider here), the conductive current is just another representation of the magnetic field. Specifically, the current density can be written in terms of the magnetic field as (Ampère’s law):

(48)

4.1 Poloidal current density

Use Eq. (48) and the definition g ≡ RB_ϕ, the poloidal components of the current density, J_Z and J_R, can be written as

(49)

and

(50)

respectively.

4.2 Toroidal current density

Ampere’s law (48) indicates the toroidal current density J_ϕ is given by

(52)

which is the Laplace operator in cylindrical coordinates for the axisymmertic case, then Eq. (51) is written as

(53)

(Many authors call Eq. (53) the Grad-Shafranov equation. This is incorrect. This is just Ampere’s law, which has nothing to do with the force-balance.)

5 Constraint of force-balance on magnetic field

Let us consider what constraint the force balance imposes on the axisymmetric magnetic field discussed above. The MHD momentum equation is given by

(54)

where ρ, ρ_q, ℙ, J, E, and B are mass density, charge density, thermal pressure tensor, current density, electric field, and magnetic field, respectively. The electric field force ρ_qE is usually ignored due to either ρ_q = 0 or E = 0. Further assume that there is no plasma flow (u = 0, the flow effect is discussed in A.2) and the plasma pressure is isotropic, then the steady state momentum equation (force balance equation) is written as

(55)

where P is the scalar plasma pressure.

Is the force balance (55) always satisfied in a real toakamak discharge? To answer this question, we need to go back to the original momentum equation (54). The imbalance between J × B and ∇P will give rise to the compressional Alfven waves, the time-scale of which, τ_A, is much shorter than the time-scale τ we are interested in. Therefore, on the time scale τ (and for slow flow with u < C_s, where C_s is the the sound speed), the leading order of the momentum equation is the force balance (55). (to be sure, I need to think about this again). This reasoning is from Youwen Sun[25].

5.1 Parallel force balance

Consider the force balance in the direction of B. Dotting the equilibrium equation (55) by B, we obtain

(56)

which implies that P is constant along a magnetic field line. Since Ψ is also constant along a magnetic field line, P can be expressed in terms of only Ψ on a single magnetic line. Note that this does not necessarily mean P is a single-valued function of Ψ, (i.e. P = P(Ψ)). This is because P still has the freedom of taking different values on different magnetic field lines with the same value of Ψ while still satisfying B ⋅∇P = 0. This situation can appear when there are saddle points (X points) in Ψ contours (refer to Sec. A.10) and P takes different functions of Ψ in islands of Ψ sepearated by a X point. For pressure within a single island of Ψ, P = P(Ψ) can be safely assumed.

On the other hand, if P = P(Ψ), then we obtain

B ⋅∇P = B ⋅∇Ψ = 0,

i.e., Eq. (56) is satisfied, indicating P = P(Ψ) is a sufficient condition for the force balance in the parallel (to the magnetic field) direction.

5.2 Toroidal force balance

Consider the force balance in the toroidal direction. The ϕ component of Eq. (55) is written

(57)

Since P = P(Ψ), which implies ∂P∕∂ϕ = 0, equation (57) reduces to

(58)

Using the expressions of the poloidal current density (49) and (50) in the force balance equation (58) yields

(59)

which can be further written

(60)

According to the same reasoning for the pressure, we conclude that g = g(Ψ) is a sufficient condition for the toroidal force balance. (The function g defined here is usually called the “poloidal current function” in tokamak literature. The reason for this name is discussed in Sec. A.3.)

5.3 Force balance along the major radius

Consider the force balance in direction. The component of Eq. (55) is written

(61)

Using the expressions of the current density and magnetic field [Eqs. (6) and (50)], equation (61) is written

(62)

Assuming the sufficient condition discussed above, i.e., P and g are a function of only Ψ, i.e., P = P(Ψ) and g = g(Ψ), Eq. (62) is written

(63)

which can be simplified to

(64)

which is the requirement of force-balance along the major radius. On the other hand, we know J_ϕ can be expressed in Ψ via Eq. (53). Combining this with Eq. (64) yields

(65)

i.e.,

(66)

Equation (66) is known as Grad-Shafranov (GS) equation.

[Note that the Z component of the force balance equation is written

JRBϕ − JϕBR = ⇒− −△∗Ψ = μ0 ⇒− −△∗Ψ = μ0 ⇒− −△∗Ψ = μ0 ⇒△∗Ψ = −μ0R2 −g

which turns out to be identical with the Grad-Shafranov equation. This is not a coincidence. The reason is that the force balance equation has been satisfied in three different directions (namely, , , and B direction) and thus it must be satisfied in all the directions.]

5.4 Axisymmetric equilibrium magnetic field

A general axisymmetric magnetic field (which does not necessarily satisfy the force balance), is given by Eq. (10), i.e.,

(67)

For the above axisymmetric magnetic field to be consistent with the force balance equation (55), there are additional requirements for Ψ and g. Specifically, Ψ is restricted by the GS equation and g should be a function of only Ψ. Therefore an axisymmetric equilibrium magnetic field is fully determined by two functions, Ψ = Ψ(R,Z) and g = g(Ψ). The Ψ is determined by solving the GS equation with specified RHS source terms and boundary conditions.

The RHS source terms in the GS equation (66) are P(Ψ) and g(Ψ), both of which must be specified before the GS equation can be solved. For most cases, the source terms are nonlinear about Ψ and thus the GS equation is a two-dimensional (in R and Z) nonlinear partial differential equation for Ψ.

For most choices of P(Ψ) and g(Ψ), the GS equation (66) has to be solved numerically. For some particular choices of P and g profiles, analytical solutions are available, one of which is the Solovév equilibrium and is discussed in Appendix A.1.

Note that we solve the GS equation in order to obtain the poloidal magnetic flux Ψ and thus the poloidal magnetic field. The toroidal magnetic field must be specified in some way before we can solve the GS equation. There are several ways of specifying the toroidal magnetic field: (1) given g(Ψ), (2) given ⟨j_∥⟩, (3) given the safety factor q(Ψ). There are simple relations between g, ⟨j_∥⟩, and q, which allows translation form one to another (discussed later). In transport simulations, ⟨j_∥⟩ is obtained from current drive models and neoclassical bootstrap current models. Note that the specification of the source terms (P, g, q, and ⟨j_∥⟩) usually involve the unknown Ψ (via not only the explicit presence of Ψ, but also the flux-surface averaging which implicit involves Ψ). This indicates that iterations are needed when numerically solving the GS equation.

5.5 Axisymmetric equilibrium current density

Since J = μ₀⁻¹∇× B, the current density J can be inferred from a given magnetic field. The components of J (expressed in terms of g and Ψ) are given by Eqs. (49), (50) and (53) and these expressions can be further simplified by using the equilibrium constraints, such as g = g(Ψ) and △^∗Ψ = −μ₀R² −g. The simplified J expressions are given in A.8.

On the other hand, in the kinetic equilibrium reconstruction[19], it is the current that is first (partially) specified (e.g., by summing sources of current drive and bootstrap current) and then the current is used as constraints for the GS equation, i.e, constraints for the magnetic field.

5.6 Equilibrium scaling

The GS equation is given by Eq. (66), i.e.,

(68)

If a solution to the GS equation is obtained, the solution can be scaled to obtain a family of solutions. Given an equilibrium with Ψ(R,Z), P(Ψ), g(Ψ), then it is ready to prove that Ψ₂ = sΨ, P₂ = s²P(Ψ), and g₂ = ±sg(Ψ) is also a solution to the GS equation, where s is a constant. In this case, both the poloidal and toroidal magnetic fields are increased by a factor of s, and thus the safety factor remains unchanged. Also note that the pressure is increased by s² factor and thus the value of β (the ratio of the therm pressure to magnetic pressure) remains unchanged. Note that g₂ = ±sg(Ψ), which indicates that the direction of the toroidal magnetic field can be reversed without breaking the force balance. Also note that Ψ₂ = sΨ and s can be negative, which indicates that the direction of the toroidal current can also be reversed without breaking the force balance.

The second kind of scaling is to set Ψ₂ = Ψ, P₂ = P(Ψ), and g₂² = g²(Ψ) + c. It is ready to prove that the scaled expression is still a solution to the GS equation because g₂g₂′ = gg′. This scaling keep the pressure and the poloidal field unchanged and thus the poloidal beta β_p remains unchanged. This scaling scales the toroidal field and thus can be used to generate a family of equilibria with different profiles of safety factor.

Another scaling, which is trivial, is to set Ψ₂ = Ψ, P₂ = P(Ψ) + c, and g₂ = g(Ψ). This scaling can be used to test the effects of the pressure (not the pressure gradient) on various physical processes.

When a numerical equilibrium is obtained, one can use these scaling methods together to generate new equilibria that satisfy particular global conditions. Note that the shape of magnetic surfaces of the scaled equilibrium remains the same as the original one.

The above scaling is made under the constraint that the the GS equation (68), is satisfied. In practice, we may scale Ψ by a factor while fixing g and P. This does not satisfy the GS equation, but allows more flexibility in changing the the safety factor profile. Scaling Ψ by a factor corresponds to scaling the plasma current and the poloidal magnetic field.

6 Numerical methods of reconstructing equilibrium magnetic field

6.1 By using integral form and Picard iteration

Given all the toroidal currents, Ψ can be calculated using a formula similar to the Biot-Savart Law (see (667)):

(69)

where J_ϕ is the plasma toroidal current density, I_i is the toroidal curent in No. i PF coil located at (R = R_i,Z = Z_i), G is a function given by

(71)

Substituting this into Eq. (69) gives an implicit formula for Ψ, which can be iterated (assuuming an initial guess of Ψ(R,Z) on the rhs and assuming the function forms, P(Ψ) and g(Ψ), are known). Will the iteration turn out to converge? We do not know for sure, and numerical experiments indicate it does in some cases[17]. Let us discuss these cases. The 1D functions, P(Ψ) and g(Ψ), are modeled as polynomials in Ψ with coefficients α_n and β₀ to be determined:

(72)

(73)

where x = Ψ = (Ψ − Ψ_M)∕(Ψ_B − Ψ_M), Ψ_M and Ψ_B are values of Ψ at the magnetic axis and LCFS (values of Ψ_M and Ψ_B are figured out from Ψ(R,Z)).

Choosing an initial guess of Ψ(R,Z) on gridpoints, we can get the values of α_n, β₀, and I_i by solving a least square problem that minimises the difference between quantities computed (from Ψ and I_i) and the corresponding quantities measured in actual experiments. After this, the rhs of Eq. (69) is fully known, and we can use this equation to update Ψ(R,Z). Then, using the latest Ψ, we solve the least square problem again. Then we obtain new Ψ using Eq. (69) again. We repeat the above procedure until convergecne in Ψ(R,Z). This is called Picard iteration.

I use the above method in developing a numerical code called HEQ. The following subsection discusses the least square problem in details.

(74)

(75)

(76)

(77)

(78)

6.2 Response matrix Γ

Plugging expressions (72) and (73) into (71), and organizing the result according to the coefficients (α₀,α₁,α₂,β₀), we obtain

then Eq. (79) is written as

(80)

Plugging expression (80) into Eq. (69), we obtain

(82)

for (j = 0,…,3), and

(83)

for (j = 4,…,N − 1) . Here (R_i^(FL),Z_i^(FL)) is the location where the ith flux measurement is made.

The matrix Γ is called “response matrix”, which is of shape (M,N), The value of M is equal to number of measurements included (measurements are merged to Γ by row stacking). We will include M_Ψ measurements of the poloidal flux, M_B measurements of the poloidal magnetic field, and 1 measurement of plasma current. Then M = M_Ψ + M_B + 1. (For EAST, M_Ψ = 35,M_B = 38.)

Next, we consider the poloidal field measurement:

(84)

where _i is the positive normal direction of the No. i magnetic probe, _i,R = _i ⋅e_R, _i,Z = _i ⋅e_Z. On the other hand, B_p(R_i,Z_i) can be inferred from Ψ via

(85)

where

(86)

and G_BR and G_BZ is given by Eq. (671) and (672). The rhs of Eq. (85) indicates that the matrix elements of Γ for (i = M_Ψ,…,M_Ψ + M_B − 1) are given by:

(87)

for (j = 0,…,3), and

(88)

for (j = 4,…,N − 1) .

Next, we consider the constraint of plasma current. The plasma current I_p can be inferred from Ψ via

(89)

The rhs of Eq. (89) indicates that the matrix elements of Γ for i = M_Ψ + M_B + 1 are given by:

(90)

for (j = 0,…,3), and zero for (j = 4,…,n − 1) .

6.3 Algorithm summary

Step 0. Initial guess of Ψ: Ψ^(k) with k = 0.

Step 1. Compute elements of response matrix Γ using Ψ^(k).

Step 2. Solve linear least-square problem to get (c₁,c₂,c₂,c₄,I₁,…,I_Nc).

Step 3. Update Ψ using the Picard iteration:

(91)

or in its discrete form:

(92)

Step 4: k = k + 1 and goto Step 1.

In Step 1, we first search for the magnetic axis and boundary surface of Ψ^(k) so that we can get the values of Ψ^(k) at those locations. These values are needed in computing x = (Ψ^(k) − Ψ_M)∕(Ψ_B − Ψ_M).

A formula for computing q value at the magnetic axis [ref. P. M. Bellan’s paper]:

qaxis = ,

where

(93)

Jϕ axis = ∑ j=03c jbj(x) = c0b0(0) + c3b3(0) = c0Raxis + c3

(94)

(95)

6.4 Solving by inverting the Laplace operator–to be deleted

If we want to invert the Laplace operator △^∗ of the GS equation to solve for Ψ, we encounter two issues: (1) the RHS of the GS equation involves functions of Ψ, where both the functions and Ψ are unknown; (2) the values of Ψ on a given boundy should be prescribed, which in turn invovles two issues: (a) how to choose the shape of the bounday? (b) how to know the value of Ψ on the boundary? The answer to isssues (a) and (b) is that we can choose any boundy we like (e.g., rectangular region) and guess the values of Ψ on it (which will be updated in iterations discussed later). The answer to issue (1) is that the function forms are chosen by users, usually as polymials in Ψ, e.g. [17]

where x = (Ψ − Ψ_M)∕(Ψ_B − Ψ_M), Ψ_M and Ψ_B are values of Ψ at the magnetic axis and LCFS (but the locations of the axis and LCFS are unknown), the free parameters α_n and β_n are assumed known (which will be later adjusted in the least-square method to make the resulting solution approximately match some measurements and constraints imposed). To make the RHS fully known, we need to guess values of Ψ in the inner region. Then inverting  △^∗ to get new Ψ. The values we obtained via inerting △^∗ usually differ from what we used in the RHS, so we need to interate until convergence. We call this interation as the first-level iteration.

After the first-level iteration coverges, we can calculate J_ϕ through Eq. (64) or Eq. (53). After this, all the current (current in the plasma, and current in the external coils, which are assumed to be known and will be adjusted in the least-square method) perpendicular to the poloidal plane is known, we can calculate the value of Ψ on the boundary, Ψ_b, by using the Green’s function method. If Ψ_b calculated this way differs from the initial guess of the value of Ψ on the boundary, we need to iterate: use the Ψ_b calculated as a new guess value of Ψ on the boundary and repeat the above procedures until convergence. We call this iteration as the second-level iteration.[?].

The thrid level iteration is to iterate over the free parameters α_n, β_n, the PF coil currents, to make the solution match match some measurements and constraints imposed in the least square sence. This iteration is an optimization problem.

7 Fixed boundary equilibrium and choices of coordinates

Plasma physists love to consider a kind of simplified problem: the fixed boundary equilibrium problem, where the shape of the boundary flux surface is given (the value of Ψ is a constant on this boundary). In dealing with the fixed boundary problem, the curvilinear coordinate system is useful. Specifically, the convenience is that the coordinates can be adjusted to make one of the coordinate surfaces coincide with the given boundary flux surface, so that the boundary condition becomes trivial.

More generally, curvilinear coordinates can be ajusted to make coordinate surfaces coincide with magnetic surfaces (dicussed later). This often simplifies analysis of wave and transport problems, especially when we properly choose the poloidal/toroidal coordinate to make magnetic field lines look like straight lines in terms of these coordinates.

Next section discusses the basic theory of curvilinear coordinates system[4].

In many studies of tokamak plasmas, one need construct a curvilinear coordinate system based on a given magnetic cofiguration in order to make the problem amenable to analytical methods or numerical methods. Specifically, many theories and numerical codes use the curvilinear coordinate systems that are constructed with one coordinate surface coinciding with magnetic surfaces. In these coordinate systems, we need to choose a poloidal coordinate 𝜃 and a toroidal coordinate ζ. As metioned above, a particular choice for 𝜃 and ζ is one that makes the magnetic field lines be straight lines in (𝜃,ζ) plane. These kinds of coordinates are often called magnetic coordinates. That is, “magnetic coordinates  are defined so they conform to the shape of the magnetic surfaces and trivialize the equations for the field lines.”

A further tuned magnetic coordinate system is the so-called field aligned (or filed-line following) coordinate system, in which changing one of the three coordinates with the other two fixed would correspond to following a magnetic field line. The field aligned coordinates are discussed in Sec. 15.

Next, let us discuss some general properties about coordinates transformation.

8 Curvilinear coordinate system

8.1 Coordinates transformation

In the Cartesian coordinates, a point is described by its coordinates (x,y,z), which, in the vector form, is written as

(96)

where r is the location vector of the point; , , and are the basis vectors of the Cartesian coordinates, which are constant, independent of spactial location. The transformation between the Cartesian coordinates system, (x,y,z), and a general coordinates system, (x₁,x₂,x₃), can be expressed as

(97)

For example, cylindrical coordinates (R,ϕ,Z) can be considered as a general coordinate systems, which are defined by  r = R cosϕ + R sinϕ + Z.

The transformation function in Eq. (97) can be written as

(98)

8.2 Jacobian

A useful quality characterizing coordinate transformation is the Jacobian determinant (or simply called Jacobian), which, for the transformation in Eq. (98), is defined by

(99)

which can also be written as

(100)

It is easy to prove that the Jacobian 𝒥 in Eq. (100) can also be written (the derivation is given in my notes on Jacobian)

(101)

Conventionally, the Jacobian of the transformation from the Cartesian coordinates to a particular coordinate system σ is called the Jacobian of σ, without explitly mentioning that this transformation is with respect to the Cartesian coordinates.

Using the defintion in Eq. (99), the Jacobian 𝒥 of the Cartesian coordinates can be calculated, yielding 1. Likewise, the Jacobian of the cylindrical coordinates (R,ϕ,Z) can be calculated as follows:

8.3 Orthogonality relation between two sets of basis vectors

In a curvilinear coordinate system (x₁,x₂,x₃), there are two kinds of basis vectors: ∇x_i and ∂r∕∂x_i, with i = 1,2,3. These two kinds of basis vectors satisfy the following orthogonality relation:

(102)

where δ_ij is the Kronical delta function. [Proof: Working in a Cartesian coordinate system (x,y,z) with the corresponding basis vectors denoted by (,,), then the left-hand side of Eq. (102) can be written as

[The cylindrical coordinate system (R,ϕ,Z) is an example of general coordinates. As an exercise, we can verify that the cylindrical coordinates have the property given in Eq. (102). In this case, x = x₁ cosx₂, y = x₁ sinx₂, z = x₃, where x₁ ≡ R, x₂ ≡ ϕ, x₃ ≡ Z.]

It can be proved that ∇xⁱ is a contravariant vector while ∂r∕∂xⁱ is a covariant vector (I do not prove this and do not bother with the meaning of these names, just using this as a naming scheme for easy reference).

The orthogonality relation in Eq. (102) is fundamental to the theory of general coordinates. The orthogonality relation allows one to write the covariant basis vectors in terms of contravariant basis vectors and vice versa. For example, the orthogonality relation tells that ∂r∕∂x₁ is orthogonal to ∇x₂ and ∇x₃, thus, ∂r∕∂x₁ can be written as

(105)

where A is a unknown variable to be determined. To determine A, dotting Eq. (105) by ∇x₁, and using the orthogonality relation again, we obtain

(106)

which gives

(108)

Similarly, we obtain

(109)

and

(110)

Equations (108)-(110) can be generally written

(111)

where (i,j,k) represents the cyclic order in the variables (x₁,x₂,x₃). Equation (111) expresses the covariant basis vectors in terms of the contravariant basis vectors. On the other hand, from Eq. (108)-(110), we obtain

(112)

which expresses the contravariant basis vectors in terms of the covariant basis vectors.

8.4 An example: (ψ,𝜃,ζ) coordinates

Suppose (ψ,𝜃,ζ) is an arbitrary general coordinate system. Following Einstein’s notation, contravariant basis vectors are denoted with upper indices as

(115)

In term of the contravairant basis vectors, A is written

(116)

where the components are easily obtained by taking scalar product with e_ψ,e_𝜃,ande_ζ, yielding A_ψ = A ⋅ e_ψ, A_𝜃 = A ⋅ e_𝜃, and A_ζ = A ⋅ e_ζ. Similarly, in term of the covariant basis vectors, A is written

(117)

where A^ψ = A ⋅ e^ψ, A^𝜃 = A ⋅ e^𝜃, and A^ζ = A ⋅ e^ζ.

Using the above notation, the relation in Eq. (111) is written as

(118)

(119)

(120)

where  𝒥 = [(∇ψ ×∇𝜃) ⋅∇ζ]⁻¹. Similarly, the relation in Eq. (112) is written as

(121)

(122)

(123)

8.5 Gradient and directional derivative in general coordinates (ψ,𝜃,ζ)

The gradient of a scalar function f(ψ,𝜃,ζ) is readily calculated from the chain rule,

(124)

Note that the gradient of a scalar function is in the covariant representation. The inverse form of this expression is obtained by dotting the above equation respectively by the three contravariant basis vectors, yielding

(125)

(126)

(127)

Using Eq. (124), the directional derivative in the direction of ∇ψ is written as

(128)

8.6 Divergence operator in general coordinates (ψ,𝜃,ζ)

To calculate the divergence of a vector, it is desired that the vector should be in the contravariant form because we can make use of the fact:

(129)

for any scalar quantities α and β. Therefore we write vector A as

(130)

where A^(ψ) = A ⋅∇ψ, A^(𝜃) = A ⋅∇𝜃, A^(ζ) = A ⋅∇ζ.   Then the divergence of A is readily calculated as

8.7 Laplacian operator in general coordinates (ψ,𝜃,ζ)

The Laplacian operator is defined by ∇² ≡∇⋅∇. Then ∇²f is written as (f is an arbitrary function)

8.8 Curl operator in general coordinates (ψ,𝜃,ζ)

To take the curl of a vector, it should be in the covariant representation since we can make use of the fact that ∇×∇α = 0. Thus the curl of A is written as

8.9 Metric tensor for general coordinate system

Consider a general coordinate system (ψ,𝜃,ζ). I define the metric tensor as the transformation matrix between the covariant basis vectors and the contravariant ones. Equations (111) and (112) express the relation between the two sets of basis vectors using cross product. Next, let us express the relation in matrix from. To obtain the metric matrix, we write the contrariant basis vectors in terms of the covariant ones, such as

(139)

Taking the scalar product respectively with ∇ψ, ∇𝜃, and ∇ζ, Eq. (139) is written as

(140)

(141)

(142)

Similarly, we write

(143)

Taking the scalar product with ∇ψ, ∇𝜃, and ∇ζ, respectively, the above becomes

(144)

(145)

(146)

The same situation applies for the ∇ζ basis vector,

(147)

Taking the scalar product with ∇ψ, ∇𝜃, and ∇ζ, respectively, the above equation becomes

(148)

(149)

(150)

Summarizing the above results in matrix form, we obtain

(151)

Similarly, to convert contravariant basis vector to covariant one, we write

(152)

Taking the scalar product respectively with ∇𝜃 ×∇ζ𝒥 , ∇ζ ×∇ψ𝒥 , and ∇ψ ×∇𝜃𝒥 , the above equation becomes

(153)

(154)

(155)

For the second contravariant basis vector

(156)

(157)

(158)

(159)

For the third contravariant basis vector

(160)

(161)

(162)

(163)

Summarizing these results, we obtain

(164)

where

M = ,

This matrix and the matrix in Eqs. (151) should be the inverse of each other. It is ready to prove this by directly calculating the product of the two matrix.

Special case: metric tensor for (ψ,𝜃,ϕ) coordinate system

Suppose that (ψ,𝜃,ϕ) are arbitrary general coordinates except that ϕ is the usual toroidal angle in cylindrical coordinates. Then ∇ϕ = 1∕R is perpendicular to both ∇ψ and ∇𝜃. Using this, Eq. (151) is simplified to

(165)

Similarly, Eq. (164) is simplified to

(166)

[Note that the matrix in Eqs. (165) and (166) should be the inverse of each other. The product of the two matrix,

(167)

can be calculated to give

,

where

By using the definition of the Jacobian in Eq. (101), it is easy to verify that A = 1, i.e.,

(168)

]

9 Covariant/contravariant representation of equilibrium magnetic field

The axisymmetric equilibrium magnetic field is given by Eq. (67), i.e.,

(169)

In a general coordinate system (ψ,𝜃,ϕ) (not necessarily magnetic surface coordinates), the above expression can be written as

(170)

where the subscripts denote the partial derivatives with the corresponding subscripts. Note that Eq. (170) is a mixed representation, which involves both covariant and contravariant basis vectors. Equation (170) can be converted to the contravariant form by using the metric tensor, giving

(171)

Similarly, Eq. (170) can also be transformed to the covariant form, giving

(172)

For the convenience of notation, define

(173)

then Eq. (172) is written as

(174)

10 Magnetic surface coordinates (ψ,𝜃,ϕ)

A coordinate system (ψ,𝜃,ϕ), where ϕ is the usual cylindrical toroidal angle, is called a magnetic surface coordinate system if Ψ is a function of only ψ, i.e., ∂Ψ∕∂𝜃 = 0 (we also have ∂Ψ∕∂ϕ = 0 since we are considering axially symmetrical case). In terms of (ψ,𝜃,ϕ) coordinates, the contravariant form of the magnetic field, Eq. (171), is written as

(175)

where Ψ′≡ dΨ∕dψ. The covariant form of the magnetic field, Eq. (172), is written as

(176)

10.1 Local safety factor

The local safety factor is defined by

(177)

which characterizes the local pitch angle of a magnetic field line in (𝜃,ϕ) plane of a magnetic surface. Substituting the contravariant representation of the magnetic field, Eq. (175), into the above equation, the local safety factor is written

(178)

Note that the expression in Eq. (178) depends on the Jacobian 𝒥 . This is because the definition of depends on the definition of 𝜃, which in turn depends on the the Jacobian 𝒥 .

In terms of , the contravariant form of the magnetic field, Eq. (175), is written

(179)

and the parallel differential operator B₀ ⋅∇ is written as

(180)

If happens to be independent of 𝜃 (i.e., field lines are straight in (𝜃,ϕ) plane), then the above operator becomes a constant coefficient differential oprator (after divided by 𝒥⁻¹). This simplification is useful because different poloidal harmonics are decoupled in this case. We will discuss this issue futher in Sec. 14.

10.2 Global safety factor

The global safety factor defined in Eq. (36) is actually the poloidal average of the local safety factor, i.e.,

Next, let us transform the 𝜃 integration in expression (182) to a curve integral in the poloidal plane. Using the relation dℓ_p and d𝜃 [Eq. (190)], expression (182) is further written

(184)

We note that expression (184) involves magnetic surface averaging, which is unknown before we know Ψ. Therefore iteration is usually needed in solving the fixed-q equilibrium (i.e., we guess the unknown Ψ, so that the magnetic surface averaging in expression (184) can be performed, yielding the values of g.)

Using B_p = |∇Ψ|∕R and B_ϕ = g∕R,  the absolute value of q in expression (183) is written

10.3 Relation between Jacobian and poloidal angle 𝜃

Given the definition of a magnetic surface coordinate system (ψ,𝜃,ϕ), the Jacobian of this system is fully determined. On the other hand, given the definition of ψ, ϕ, and the Jacobian, the definition of 𝜃 is fully determined (can have some trivial shifting freedoms). Next, let us discuss how to calculate 𝜃 in this case. In (ψ,𝜃,ϕ) coordinates, a line element is written

(187)

The line element that lies on a magnetic surface (i.e., dψ = 0) and in a poloidal plane (i.e., dϕ = 0) is then written

(189)

Using the fact that ∇ψ and ∇ϕ are orthogonal and ∇ϕ = ∕R, the above equation is written as

(190)

Given |𝒥∇ψ|, Eq. (190) can be integrated to determine the 𝜃 coordinate of points on a magnetic surface.

—–

(191)

10.4 Calculating poloidal angle

Once |𝒥∇ψ| is known, the value of 𝜃 of a point can be obtained by integrating expression (190), i.e.,

(192)

where the curve integration is along the contour Ψ = Ψ_j, x_ref,j is a reference point on the contour, where value of the poloidal angle is chosen as 𝜃_ref,j. The choice of the positive direction of 𝜃 is up to users. Depending on the positive direction chosen, the sign of the Jacobian of the constructed coordinates can have a sign difference from the 𝒥 appearing in Eq. (192). Denote the Jacobian of the constructed coordinates by 𝒥′, then

(193)

This sign can be determined after the radial coordinate and the positive direction of the poloidal angle are chosen. In GTAW code, I choose the positive direction of 𝜃 to be in anticlockwise direction when observers look along the direction of . To achieve this, the line integration in Eq. (192) should be along the anticlockwise direction. (I use the determination of the direction matrix, a well known method in graphic theory, to determine the direction from a given set of discrete points on a magnetic surface.)

Normalized poloidal angle

The span of 𝜃 defined by Eq. (192) is usually not 2π in one poloidal loop. This poloidal angle can be scaled by s(ψ) to define a new poloidal coordinate 𝜃, whose span is 2π in one poloidal loop, where s(ψ) is a magnetic surface function given by

(194)

Then 𝜃 is written as

[The reference points x_ref,j and the values of poloidal angle at these points can be chosen by users. One choice of the reference points x_ref,j are those points on the horizontal ray in the midplane that starts from the magnetic axis and points to the low filed side of the device and 𝜃_ref,j at these points is chosen as zero (this is my choice in the GTAW code). In the TEK code, the reference points are chosen at the high-field side of the midplane and 𝜃_ref,j = −π at the reference points.]

Jacobian for equal-arc-length poloidal angle

If the Jacobian 𝒥 is chosen to be of the following form

(197)

Then 𝜃_i,j in Eq. (195) is written

(198)

and the Jacobian of new coordinates (ψ,𝜃,ϕ), 𝒥_new, which is given by Eq. (196), now takes the form

(199)

Equation (198) indicates a set of poloidal points with equal arc intervals corresponds to a set of uniform 𝜃_i points. Therefore this choice of the Jacobian is called the equal-arc-length Jacobian. Note that Eq. (198) does not involve the radial coordinate ψ. Therefore the values of 𝜃 of points on any magnetic surface can be determined before the radial coordinate is chosen.

Jacobian for equal-volume poloidal angle (Hamada coordinate)

The volume element in (ψ,𝜃,ϕ) coordinates is given by dV = |𝒥|d𝜃dϕdψ. If we choose a Jacobian that is independent of 𝜃, then uniform 𝜃 grids will correspond to grids with uniform volume interval. In this case, 𝒥 is written as

(200)

where h(ψ) is a function independent of 𝜃. Then 𝜃_i,j in Eq. (195) is written

(201)

and the Jacobian of the new coordinates (ψ,𝜃,ϕ), 𝒥_new,  is given by Eq. (196), which now takes the following form:

(202)

Note that both 𝜃_i,j and 𝒥_new are independent of the function h(ψ) introduced in Eq. (200). (h(ψ) is eliminated by the normalization procedure specified in Sec. 10.4 due to the fact that h(ψ) is constant on a magnetic surface.) The equal-volume poloidal angle is also called Hamada poloidal angle.

The equal-volume poloidal angle is useful in achieving loading balance for parallel particle simulations. Assume that markers are loaded uniform in space and the poloidal angle is domain decomposed and assigned to different MPI process. Then the equal-volume poloidal angle can make marker number in each MPI process be equal to each other and thus work loading to each process be equal. (**check**If the domain decomposition is also applied to the radial direction, to achieve loading balance, then the radial coordinate ψ should be chosen in a way that makes ∮ (R∕∇ψ)dl_p be independent of ψ, so that 𝒥_new in Eq. (202) is constant in space.**)

Jacobian of Boozer’s poloidal angle

If the Jacobian 𝒥 is chosen to be of the Boozer form:

(203)

then the poloidal angle in Eq. (195) is written as

(204)

The final Jacobian is given by

(205)

The usefullness of Boozer poloidal angle will be further discussed in Sec. 14.8 after we introduce a gneralized toroidal angle.

Jacobian for PEST poloidal angle (straight-field-line poloidal angle)

If the Jacobian 𝒥 is chosen to be of the following form

then Eq. (178) implies that the local safety factor, (ψ,𝜃) = −g∕Ψ′, is a magnetic surface function, i.e., the magnetic field lines are straight in (ψ,𝜃) plane. Then the poloidal angle in Eq. (195) is written

(206)

The Jacobian 𝒥_new given by Eq. (196) now takes the form

(207)

Let us denote an arbitrary poloidal angle by 𝜃 and the above straight-field-line poloidal angle by 𝜗, then it is ready to find the following relation between 𝜃 and 𝜗:

(208)

where is the local safety factor corresponding to the arbitrary poloidal angle 𝜃, i.e., = B⋅∇ϕ∕(B⋅∇𝜃). [Proof: Using d𝜃 = dl_p, the poloidal angle 𝜗 given in Eq. (206) is written as

Note that Boozer poloidal angle is very close to the poloidal angel disccused here because the two Jacobians are very similar:

(210)

Comparison between different types of poloidal angle

All Jacobians introduced above can be written in a general form:

(211)

The choice of (i = 2,j = k = 0) gives the PEST coordinate, (i = j = 0,k = 2) give the Boozer coordinate, (i = j = 1,k = 0) gives the equal-arc coordinate, (i = j = k = 0) gives the Hammada coordinate.

Figure 7 compares the equal-arc-poloidal angle and the straight-line poloidal angle, which shows that the resolution of the straight-line poloidal angle is not good near the low-field-side midplane. Since ballooning modes take larger amplitude near the low-field-side midplane, better resolution is desired there. This is one reason that I often avoid using the straight-line poloidal angle in my numerical codes.

Fig. 7: The equal-arc poloidal angle (left) and the straight-line poloidal angle (right) for EAST equilibrium #38300@3.9s. The blue lines correspond to the equal-𝜃 lines, with uniform 𝜃 interval between neighbour lines. The red lines correspond to the magnetic surfaces, which start from = 0.01714 (the innermost magnetic surface) and end at = 0.9851 (the boundary magnetic surface), and are equally spaced in .

Verification of Jacobian

After the magnetic coordinates are constructed, we can evaluate the Jacobian 𝒥_new by using directly the definition of the Jacobian, i.e.,

(212)

which can be further written as

(213)

where the partial differential can be evaluated by using numerical differential schemes. The results obtained by this way should agree with results obtained from the analytical form of the Jacobian. This consistency check provide a verification for the correctness of the theory derivation and numerical implementation. In evaluating the Jacobian by using the analytical form, we may need to evaluate ∇ψ, which finally reduces to evaluating ∇Ψ. The value of |∇Ψ| is obtained numerically based on the numerical data of Ψ given in cylindrical coordinate grids. Then the cubic spline interpolating formula is used to obtain the value of |∇Ψ| at desired points. (𝒥_new calculated by the second method (i.e. using analytic form) is used in the GTAW code; the first methods are also implemented in the code for the benchmark purpose.) In the following sections, for notation ease, the Jacobiban of the constructed coordinate system will be denoted by 𝒥 , rather than 𝒥_new.

Radial coordinate

The radial coordinate ψ can be chosen to be various surface function, e.g., volume, poloidal or toroidal magnetic flux within a magnetic surface. The frequently used radial coordinates include Ψ, and , where Ψ is defined by

(214)

where Ψ₀ and Ψ_a are the values of Ψ at the magnetic axis and LCFS, respectively. Other choices of the radial coordinates: the toroidal magnetic flux and its square root, Ψ_t, and , where Ψ_t and Ψ_t are defined by

(215)

and

(216)

respectively, where Ψ_t(0) and Ψ_t(1) are the values of Ψ_t at the magnetic axis and LCFS, respectively.

If ψ = , then

(217)

The cylindrical coordinates (R,ϕ,Z) is a right-hand system, with the positive direction of Z pointing vertically up. In GTAW code, the positive direction of 𝜃 is chosen in the anticlockwise direction when observers look along the direction of . Then the definition 𝒥⁻¹ = ∇ψ ×∇𝜃 ⋅∇ϕ indicates that (1) 𝒥 is negative if ∇ψ points from the magnetic axis to LCFS; (2) 𝒥 is positive if ∇ψ points from the LCFS to the magnetic axis. This can be used to determine the sign of Jacobian after using the analytical formula to obtain the absolute value of Jacobian.

If ψ = Ψ or ψ = , ∇ψ points from the magnetic axis to LCFS.

The volume between magnetic surfaces can also be used as a radial coordinate. The differential volume element is written as

(218)

Integrating over the toroidal angle, we obtain

(219)

Further integrating over the poloidal angle, we obtain

(220)

i.e.,

(221)

In codes I wrote, I stick to using Ψ as the radial coordinate when doing computation, and transform to other radial coordinates when presenting the results if needed.

11 Constructing magnetic surface coordinate system from discrete Ψ(R,Z) data

Given an axisymmetric tokamak equilibrium in (R,ϕ,Z) coordinates (e.g., 2D data Ψ(R,Z) on a rectangular grids (R,Z) in G-file), we can construct a magnetic surface coordinates (ψ,𝜃,ϕ) by the following two steps. (1) Find out a series of magnetic surfaces on (R,Z) plane and select radial coordinates for each magnetic surface (e.g. the poloidal flux within each magnetic surface). (2) Specify the Jacobian or some property that we want the poloidal angle to have. Then calculate the poloidal angle of each point on each flux surface (on the ϕ = const plane) by using Eq. (190) (if the Jacobian is specified) or some method specified by us to achieve some property we prefer for the poloidal angle (if a Jacobian is not directly specified). Then we obtain the magnetic surface coordinates system (ψ,𝜃,ϕ).

11.1 Finding magnetic surfaces

Two-dimensional data Ψ(R,Z) on a rectangular grids (R,Z) is read from the G_EQDSK file (G-file) of EFIT code. Based on the 2D array data, I use 2D cubic spline interpolation to construct a interpolating function Ψ = Ψ(R,Z). To construct a magnetic surface coordinate system, I need to find the contours of Ψ, i.e., magnetic surfaces. The values of Ψ on the magnetic axis, Ψ₀, and the value of Ψ on the last closed flux surface (LCFS), Ψ_b, are given in G-file. Using these two values, I construct a 1D array “psival” with value of elements changing uniform from Ψ₀ to Ψ_b. Then I try to find the contours of Ψ with contour level value ranging from Ψ₀ to Ψ_b. This is done in the following way: construct a series of straight line (in the poloidal plane) that starts from the location of the magnetic axis and ends at one of the points on the LCFS. Combine the straight line equation, Z = Z(R), with the interpolating function Ψ(R,Z), we obtain a one variable function h = Ψ(R,Z(R)). Then finding the location where Ψ is equal to a specified value Ψ_i, is reduced to finding the root of the equation Ψ(R,Z(R)) − Ψ_i = 0. Since this is a one variable equation, the root can be easily found by using simple root finding scheme, such as bisection method (bisection method is used in GTAW code). After finding the roots for each value in the array “psival” on each straight lines, the process of finding the contours of Ψ is finished. The contours of Ψ found this way are plotted in Fig. 8.

In the above, we mentioned that the point of magnetic axis and points on the LCFS are needed to construct the straight lines. In G-file, points on LCFS are given explicitly in an array. The location of magnetic axis is also explicitly given in G-file. It is obvious that some of the straight lines Z = Z(R) that pass through the location of magnetic axis and points on the LCFS will have very large or even infinite slope. On these lines, finding the accurate root of the equation Ψ(R,Z(R)) −Ψ_i = 0 is difficult or even impossible. The way to avoid this situation is obvious: switch to use function R = R(Z) instead of Z = Z(R) when the slope of Z = Z(R) is large (the switch condition I used is |dZ∕dR| > 1).

In constructing the flux surface coordinate with desired Jacobian, we will need the absolute value of the gradient of Ψ, |∇Ψ|, on some specified spatial points. To achieve this, we need to construct a interpolating function for |∇Ψ|. The |∇Ψ| can be written as

(222)

By using the center difference scheme to evaluate the partial derivatives with respect to R and Z in the above equation (using one side difference scheme for the points on the rectangular boundary), we can obtain an 2D array for the value of |∇Ψ| on the rectangular (R,Z) grids. Using this 2D array, we can construct an interpolating function for ∇Ψ by using the cubic spline interpolation scheme.

11.2 Expression of metric elements of magnetic coordinates (ψ,𝜃,ϕ)

Metric elements of the (ψ,𝜃,ϕ) coordinates, e.g., ∇ψ ⋅∇𝜃, are often needed in practical calculations. Next, we express these metric elements in terms of the cylindrical coordinates (R,Z) and their partial derivatives with respect to ψ and 𝜃. Note that, in this case, the coordinate system is (ψ,𝜃,ϕ) while R and Z are functions of ψ and 𝜃, i.e.,

(223)

(224)

Then ∇R and ∇Z are written as

(225)

(226)

wehre R_ψ ≡ ∂R∕∂ψ, etc. Equations (225) and (226) can be solved to give

(227)

(228)

Using the above expressions, the Jacobian of (ψ,𝜃,ϕ) coordinates, 𝒥 , is written as

(230)

Using this, Expressions (227) and (228) are written as

(231)

and

(232)

Then the elements of the metric matrix are written as

(233)

(234)

and

(235)

Equations (233), (234), and (235) are the expressions of the metric elements in terms of R, R_ψ, R_𝜃, Z_ψ, and Z_𝜃. [Combining the above results, we obtain

(236)

Equation (235) is used in GTAW code. Using the above results, h^αβ = ∇α ⋅∇β are written as

(237)

(238)

(239)

As a side product of the above results, we can calculate the arc length in the poloidal plane along a constant ψ surface, dℓ_p, which is expressed as

12 Constructing model tokamak magnetic field

In some cases (e.g., turbulence simulation), model tokamak magnetic field, which is not an exact solution to the GS equation, is often used. In the model, the safety factor profile q(ψ), toroidal field function g(ψ), and the magnetic surface shape (R(ψ,𝜃),Z(ψ,𝜃)) are given. To use this model in a simulation, we need to calculate its poloidal magnetic field, which is determined by the poloidal magnetic flux function Ψ. To determine Ψ, we need to use Eq. (182), i.e.,

(241)

and re-organize the formula as

(242)

which can be integrated to obtain Ψ and thus the poloidal magnetic field B_p = ∇Ψ ×∇ϕ, where 𝒥 = [(∇ψ ×∇𝜃) ⋅∇ϕ]⁻¹. The toroidal magnetic field can be obtained from g(ψ) by B_ϕ = g∕R. In most papers, g(ψ) is chosen to be a constant.

A typical magnetic surface shape used in simulations is the Miller shape, which is given by

(243)

where r = ψ is the radial coordinate, R₀(r), δ(r), and κ_r(r) are the Shafronov shift, triangularity, and elongation profiles, which can be arbitrarily specified. For the special case of R₀ being a constant, δ(r) = 0, and κ(r) = 1, this shape reduces to a concentric-circular magnetic field. An example of calculating the poloidal field of the model field is given in Sec. 16. This kind of model field can be called theoretical physicists’ tokamak. Computational physicists often use this model for code benchmarking purpose. The famous DIII-D cyclone base case is an example, which was extensively used for benchmarking gyrokinetic simulation of ion temperature driven turbulence.

13 Magnetic surface averaging

13.1 Definition

The magnetic surface average of a physical quantity G(ψ,𝜃,ϕ) is defined by

(244)

where the volume integration is over a small volume between two adjacent flux surfaces with Ψ differing by △Ψ. [This formula (with finite ΔΨ) is used in TEK code to calculate the radial heat flux.]

The above 3D volume integration can also be written as a 2D surface integration. The differential volume element is given by d³V = |𝒥|dψd𝜃dϕ, where 𝒥 is the Jacobian of (ψ,𝜃,ϕ) coordinates. Using this, equation (244) is written as

(246)

where G₀(𝜃) is defined by the following Fourier expansion:

(247)

_________________________________________________________________________________________________________________________

“Zonal” and “mean” components

⟨G⟩ is sometimes called the “zonal” component of G if the radial wavelength of ⟨G⟩ is much smaller than the equilibrium scale length. If the radial wavelength of ⟨G⟩ is comparable to the equilibrium scale length, ⟨G⟩ is usually called “mean” component in tokamak literature. For example, mean flows are of system space scale and thus are easy to be observed in experiments. On the other hand, the “zonal” flow, which usually refers to the turbulence generated secondary flow, is of much smaller radial scale (the radial wavelength of zonal flow is of several Larmor radius) and thus is difficult to observe in experiments.

_________________________________________________________________________________________________________________________

Sometimes, we do not want the Jacobian to explicitly appear in the formula. This can be achieved by writing the differential volume element as

(248)

Using B_p = |∇Ψ|∕R, the volume element is further written as

(249)

Using this, the averaging defined in Eq. (244) is written as

(251)

(Equation (251) is used in the GTAW code to calculate the magnetic surface averaging.) Using Eq. (190) and B_p = |∇Ψ|∕R, equation (251) can also be written as

(252)

Using the expression of the volume element dτ = |𝒥|d𝜃dϕdψ, the volume within a magnetic surface is written

(253)

Using this, the differential of V with respect to ψ is written as

(254)

Using this, Eq. (252) is written as

⟨G⟩ = ∫ 02πG|𝒥|d𝜃

13.2 Flux Surface Functions—to be deleted

Next, examine the meaning of the following volume integral

(255)

where the volume V = V (ψ) is the volume within the magnetic surface labeled by ψ. Using ∇⋅ B = 0, the quantity D can be further written as

(256)

Note that 𝜃 is not a single-value function of the spacial points. In order to evaluate the integration in Eq. (256), we need to select one branch of 𝜃, which can be chosen to be 0 ≤ 𝜃 < 2π. Note that function 𝜃 = 𝜃(R,Z) is not continuous in the vicinity of the contour of 𝜃 = 0. Next, we want to use the Gauss’s theorem to convert the above volume integration to surface integration. Noting the discontinuity of the integrand 𝜃B in the vicinity of the contour of 𝜃 = 0, the volume should be cut along the contour, thus, generating two surfaces. Denote these two surfaces by S₁ and S₂, then equation (256) is written as

(258)

Using the expression of the volume element dτ = |𝒥|d𝜃dϕdψ, Ψ_p can be further written in terms of flux surface averaged quantities.

Similarly, the toroidal flux within a flux surface is written as

(260)

the poloidal current within a flux surface is written as

(261)

and toroidal current within a flux surface is written as

(262)

(**check**)The toroidal magnetic flux is written as

⇒ Ψt′ = g

⇒ = g

(264)

Next, calculate the derivative of the toroidal flux with respect to the poloidal flux.

(266)

By using the contravariant representation of current density (497), the poloidal current within a magnetic surface is written as

−∇ψ ×∇𝜃 − g′∇ϕ ×∇ψ,

The toroidal current is written as

(269)

Next, calculate another useful surface-averaged quantity,

14 Magnetic coordinates (ψ,𝜃,ζ) with general toroidal angle ζ

14.1 General toroidal angle ζ

In Sec. 10.1, we introduced the local safety factor (ψ,𝜃). Equation (178) indicates that if the Jacobian is chosen to be of the particular form 𝒥 = h(ψ)R², then the local safety factor is independent of 𝜃, i.e., magnetic line is straight in (𝜃,ϕ) plane. On the other hand, if we want to make field line straight in (𝜃,ϕ) plane, the Jacobian must be chosen to be of the specific form 𝒥 = h(ψ)R². We note that, as mentioned in Sec. 10.3, the poloidal angle is fully determined by the choice of the Jacobian. The specific choice of 𝒥 = α(ψ)R² is usually too restrictive for achieve a desired poloidal resolution (for example, the equal-arc poloidal angle can not be achieved by this choice of Jacobian). Is there any way that we can make the field line straight in a coordinate system at the same time ensure that the Jacobian can be freely adjusted to obtain desired poloidal angle? The answer is yes. The obvious way to achieve this is to define a new toroidal angle ζ that generalizes the usual toroidal angle ϕ. Define a new toroidal angle ζ by[10]

(272)

where δ = δ(ψ,𝜃) is a unknown function to be determined by the constraint of field line being straight in (𝜃,ζ) plane. Using Eq. (179), the new local safety factor in (ψ,𝜃,ζ) coordinates is written as

(274)

where c(ψ) can be an arbitrary function of ψ. A convenient choice for c(ψ) is c(ψ) = q, i.e., making the new local safety factor be equal to the original global safety factor, i.e., _new = q. In this case, equation (274) is written as

(275)

which, on a magnetic surface labed by ψ, can be integrated over 𝜃 to give

(276)

where 𝜃_ref is an starting poloidal angle arbitrarily chosen for the integration, and δ(ψ,𝜃_ref) is the constant of integration. In the following, both 𝜃_ref and δ(ψ,𝜃_ref) will be chosen to be zero. Then the above equations is written

(277)

Substituting the above expression into the definition of ζ (Eq. 272), we obtain

(278)

which is the formula for calculating the general toroidal angle. If 𝜃 is a straight-field line poloidal angle, then ζ in Eq. (278) reduces to the usual toroidal angle ϕ.

In summary, magnetic field line is straight in (𝜃,ζ) plane with slope being q if ζ is defined by Eq. (278). In this method, we make the field line straight by defining a new toroidal angle, instead of requiring the Jacobian to take particular forms. Thus, the freedom of choosing the form of the Jacobian is still available to be used later to choose a good poloidal angle coordinate. Note that the Jacobian of the new coordinates (ψ,𝜃,ζ) is equal to that of (ψ,𝜃,ϕ). [Proof:

[In numerical implementation, the term ∫ ₀^𝜃d𝜃 appearing in δ is computed by using

For later use, from Eq. (277), we obtain

14.2 Contravariant form of magnetic field in (ψ,𝜃,ζ) coordinates

Recall that the contravariant form of the magnetic field in (ψ,𝜃,ϕ) coordinates is given by Eq. (179), i.e.,

(282)

Next, let us derive the corresponding form in (ψ,𝜃,ζ) coordinates. Using the definition of ζ, equation (282) is written as

The expression of the magnetic field in Eq. (284) can be rewritten in terms of the flux function Ψ_p and Ψ_t discussed in Sec. 13.2. Equation (284) is

(285)

which, by using Eq. (259), i.e., ∇Ψ = ∇Ψ_p∕(2π), is rewritten as

(286)

which, by using Eq. (266), i.e., q = dΨ_t∕dΨ_p, is further written as

(287)

14.3 Relation between the partial derivatives in (ψ,𝜃,ϕ) and (ψ,𝜃,ζ) coordinates

Noting the simple fact that

(288)

where c is a constant, we conclude that

(289)

(since ϕ = ζ + q(ψ)δ(ψ,𝜃), where the part q(ψ)δ(ψ,𝜃) acts as a constant when we hold ψ and 𝜃 constant), i.e., the symmetry property with respect to the new toroidal angle ζ is identical with the one with respect to the old toroidal angle ϕ. On the other hand, generally,

(290)

and

(291)

In the special case that f is axisymmetric (i.e., f is independent of ϕ in (ψ,𝜃,ϕ) coordinates), then two sides of Eqs. (290) and (291) are equal to each other. Note that the partial derivatives ∂∕∂ψ and ∂∕∂𝜃 in Sec. 14.1 and 14.2 are taken in (ψ,𝜃,ϕ) coordinates. Because the quantities involved in Sec. 14.1 and 14.2 are axisymmetric, these partial derivatives are equal to their counterparts in (ψ,𝜃,ζ) coordinates.

14.4 Steps to construct a straight-line magnetic coordinate system

In Sec. 11, we have provided the steps to construct the magnetic surface coordinate system (ψ,𝜃,ϕ). Only one additional step is needed to construct the straight-line flux coordinate system (ψ,𝜃,ζ). The additional step is to calculate the generalized toroidal angle ζ according to Eq. (272), where δ is obtained from Eq. (277). Also note that the Jacobian of (ψ,𝜃,ζ) coordinates happens to be equal to that of (ψ,𝜃,ϕ) coordinates.

14.5 Form of operator B ⋅∇ in (ψ,𝜃,ζ) coordinates

The usefulness of the contravariant form [Eq. (284] of the magnetic field lies in that it allows a simple form of B ⋅∇ operator in a coordinate system. (The operator B₀ ⋅∇ is usually called magnetic differential operator.) In (ψ,𝜃,ζ) coordinate system, by using the contravariant form Eq. (284), the operator is written as

(293)

where h = h(ψ,𝜃,ζ) is some known function. Using Eq. (292), the magnetic differential equation is written as

(294)

Note that the coefficients before the two partial derivatives of the above equation are all independent of 𝜃 and ζ. This indicates that different Fourier harmonics in 𝜃 and ζ are decoupled. As a result of this fact, if f is Fourier expanded as

(295)

(note that, following the convention adopted in tokamak literature[10], the Fourier harmonics are chosen to be e^{i(m𝜃−nζ)}, instead of e^i(m𝜃+nζ)), and the right-hand side is expanded as

(296)

then Eq. (294) can be readily solved to give

(297)

The usefulness of the straight line magnetic coordinates (ψ,𝜃,ζ) lies in that, as mentioned previously, it makes the coefficients before the two partial derivatives both independent of 𝜃 and ζ, thus, allowing a simple solution to the magnetic differential equation.

14.6 Resonant surface of a perturbation

Equation (297) indicates that, for the differential equation (293), there is a resonant response to a perturbation e^{i(m𝜃−nζ)} on a magnetic surface with m − nq = 0. Therefore, the magnetic surface with q = m∕n is called the “resonant surface” for the perturbation e^{i(m𝜃−nζ)}.

The phase change of the perturbation e^{i(m𝜃−nζ)} along a magnetic field is given by mΔ𝜃 − nΔζ, which can be written as Δ𝜃(m−nq). Since m−nq = 0 on a resonant surface, this indicates that there is no phase change along a magnetic field line on a resonant surface, i.e., the parallel wavenumber k_∥ is zero on a resonant surface.

14.7 Helical angle used in tearing mode theory

Next, we discuss a special poloidal angle, which is useful in describling a perturnbation of single harmonic (m,n). This poloidal angle is defined by

(298)

where (m,n) are the mode numbers of the perturbation. The poloidal angle χ is often called helical angle and is special in that its definition is associated with a perturbation (the mode numbers of the perturbation appear in the definition) while the definition of the poloidal angles discussed previously only involve the equilibrium quantities.

The poloidal angle χ is designed to make 3D perturbations of the form ∼ f(ψ,m𝜃 −nζ) reduce to 2D perturbations, i.e.,

(299)

It is ready to verify that the Jacobian of coordinates (ψ,χ,ζ) is equal to that of coordinates (ψ,𝜃,ζ) [proof: (𝒥′)⁻¹ = ∇ψ ×∇χ ⋅∇ζ = ∇ψ ×∇(𝜃 − nζ∕m) ⋅∇ζ = ∇ψ ×∇𝜃 ⋅∇ζ = 𝒥⁻¹].

The component of B along ∇χ direction (i.e., the covariant component) is written

On the other hand, the component of B along ∇𝜃 direction is written

(302)

14.8 Covariant form of magnetic field in (ψ,𝜃,ζ) coordinate system

In the above, we have obtained the covariant form of the magnetic field in (ψ,𝜃,ϕ) coordinates (i.e., Eq. (176)). Next, we derive the corresponding form in (ψ,𝜃,ζ) coordinate. In order to do this, we need to express the ∇ϕ basis vector in terms of ∇ψ, ∇𝜃, and ∇ζ basis vectors. Using the definition of the generalized toroidal angle, we obtain

(304)

This expression can be further simplified by using equation (275) to eliminate ∂δ∕∂𝜃, which gives

(307)

with I(ψ) = h(ψ) − gq. The magnetic field expression in Eq. (307) frequently appears in tokamak literature[29]. In this form, the coefficients before both ∇𝜃 and ∇ζ depends on only the radial coordinate. In terms of I(ψ), the Jabobian can also be written as

(308)

14.9 Form of operator (B ×∇ψ∕B²) ⋅∇ in (ψ,𝜃,ζ) coordinates

In solving the MHD eigenmode equations in toroidal geometries, besides the B ⋅∇ operator, we will also encounter another surface operator (B×∇ψ∕B²) ⋅∇. Next, we derive the form of the this operator in (ψ,𝜃,ζ) coordinate system. Using the covariant form of the equilibrium magnetic field [Eq.  (306)], we obtain

(309)

Using this, the (B ×∇ψ∕B²) ⋅∇ operator is written as

Examining Eq. (311), we find that the coefficients before the two partial derivatives will be independent of 𝜃 and ζ if the Jacobian 𝒥 is chosen to be of the form 𝒥 = h(ψ)∕B², where h is some magnetic surface function. It is obvious that the independence of the coefficients on 𝜃 and ζ will be advantageous to some applications. The coordinate system (ψ,𝜃,ζ) with the particular choice of 𝒥 = h(ψ)∕B² is called the Boozer coordinates. The usefulness of the new toroidal angle ζ is highlighted in Boozer’s choice of the Jacobian, which makes both B ⋅∇ and (B ×∇ψ∕B²) ⋅∇ be a constant-coefficient differential operator. For other choices of the Jacobian, only the B⋅∇ operator is a constant-coefficient differential operator.

14.10 Radial differential operator

In solving the MHD eigenmode equations in toroidal geometry, we also need the radial differential operator ∇ψ ⋅∇. Next, we derive the form of the operator in (ψ,𝜃,ζ) coordinates. Using

∇f = ∇ψ + ∇𝜃 + ∇ζ,

the radial differential operator is written as

(313)

This formula is used in GTAW code.

15 Field-line-following coordinates

15.1 Definition of the field-line-following coordinates (ψ,𝜃,α)

In (ψ,𝜃,ζ) coordinates, a magnetic field line is straight in (𝜃,ζ) plane with slope being q. Then the equation for a magnetic field line is written as

(314)

where α is a constant in (𝜃,ζ) plane and can be used to label magnetic field lines on a magnetic surface. This motivates us to use α, i.e.,

(315)

to replace ζ. Then the magnetic field in Eq. (284) is written as

(316)

which is called the Clebsch form. The direction

(317)

is parallel (or anti-parallel) to the magnetic field direction. Due to this fact, (ψ,𝜃,α) coordinates are usually called “field-line-following coordinates” or “field-aligned coordinates” [2, 6].

Equation (316) implies that

(318)

and

(319)

i.e., both α and ψ are constant along a magnetic field line. Taking scalar product of Eq. (316) with ∇𝜃, we obtain

(320)

which is nonzero, i.e., only 𝜃 among (ψ,𝜃,α) is changing along a magnetic field line. (Here 𝒥 = (∇ψ ×∇𝜃 ⋅∇α)⁻¹ is the Jacobian of the coordinate system (ψ,𝜃,α), which happens to be equal to the Jacobian of (ψ,𝜃,ζ) coordinates.)

Using Eq. (316), the magnetic differential operator B ⋅∇ in the new coordinate system (ψ,𝜃,α) is written

(321)

which is just a partial derivative over 𝜃, as is expected, since only 𝜃 is changing along a magnetic field line.

By the way, note that (B ⋅∇α)∕(B ⋅∇𝜃) = 0, i.e., the magnetic field lines are straight with zero slope on (𝜃,α) plane.

Using Eqs. (315) and (278), α can be written as

It is widely believed that turbulence responsible for energy transport in tokamak plasmas usually has k_∥≪ k_⊥, where k_∥ and k_⊥ are the parallel and perpendicular wavenumbers, respectively. Due to this elongated structure along the parallel direction, less grids can be used in the parallel direction than that in the perpendicular direction in turbulence simulation. In this case, the field-aligned coordinates (ψ,𝜃,α) provide suitable coordinates to be used, where less gridpoints can be used for 𝜃 coordinate in simulations and even some ∂∕∂𝜃 derivatives can be neglected (high-n approximation), which simplifies the equations that need to be solved. There is another reason why almost all gyrokinetic codes use field-aligned coordinates: the  stability of numerical algorthims is improved when we use coarse grids in the parallel direction because the parallel Courant condition (for explicit schemes) Δt ≤ ΔL_∥∕v_∥ can be more easily satisfied (especially for the cases with kinetic electrons), where ΔL_∥ is the parallel grid spacing, which is larger when coarse grids are used in the parallel direction. This is also mentioned in Ref. [22] and it seems to be right from my experiences of testing several algorithm but a strict test is needed to verify this. This can also be understood in the following way: the coarse parallel grid automatically filters out physically irrelevant but numerically problematic high-k_∥ modes, permitting much longer time steps for explicit time stepping, in both particle and fluid codes[12].

15.2 Some discussions

The fact B ⋅∇α = 0 implies that α is constant along a magnetic field line. At first glance, a magnetic line on an irrational surface seems to sample all the points on the surface. This seems to indicate that α is a flux surface label for irrational surface. However, α must be a non-flux-surface-function so that it can provide a suitable toroidal coordinate. I had once been confused by this conflict for a long time. The key point to resolve this confusion is to realize that it is wrong to say there is only one magnetic line on an irrational surface, i.e. it is wrong to say a magnetic line on an irrational surface samples all the points on the surface. There are still infinite number of magnetic field lines that can not be connected with each other on an irrational surface. Then the fact B ⋅∇α = 0 does not imply that α must be the same on these different magnetic field lines. In fact, although B ⋅∇α = 0, the gradient of α on a flux-surface along the perpendicular (to B) direction is nonzero, i.e., B ×∇Ψ ⋅∇α≠0. [Proof:

In (ψ,𝜃,ϕ) coordinates, ∇ϕ is perpendicular to ∇ψ. However, in field-line-following coordinates (ψ,𝜃,α), ∇α is not perpendicular to ∇ψ. Therefore ∇α is not along the binormal direction B ×∇ψ.

15.3 Expression of ∇α, ∇ψ ⋅∇α, and ∇α ⋅∇𝜃

The generalized toroidal angle α is defined by Eq. (322), i.e., α = ϕ −δ, where δ = ∫ ₀^𝜃 d𝜃. The gradient of α is then written as

(325)

(326)

Another method of computing the above quantities: Using Eqs. (231) and (232), expression (324) is written as

15.4 Field-aligned coordinates in GEM[7] and GENE[13] codes

In GEM[7] and GENE[13] codes, the field-aligned coordinates (x,y,z) are defined by

(330)

(331)

(332)

where r is an arbitrary flux surface label with length dimension, which is often chosen in GEM to be the minor radius of a magnetic surface in the midplane. Here r₀ and R₀ are constant quantities of length dimension, r₀ is the minor radius of a reference magnetic surface (usually corresponding to the center of the radial simulation box), R₀ is the major radius of the magnetic axis, q₀ is the safety factor value on the r = r₀ surface. The constant length q₀R₀ introduced in the definition of z is to make z approximately correspond to the length along the field line in the large-aspect ratio limit. The constant length r₀∕q₀ introduced in the definition of y is to make y corresponds to the arc-length in the poloidal plane traced by a field line when its usual toroidal angle increment Δϕ is α. This explanation makes y look like a poloidal coordinate whereas y is actually a toroidal coordinate.

Next, let us calculate the wave number along the y direction, k_y, for a mode with toroidal wave number n. The wavelength along the y direction, λ_y, is given by

(333)

Then the wavenumber k_y is written as

(334)

On the other hand, the poloidal wavenumber k_𝜃 is given by

(335)

where m is the poloidal mode number of the mode in (ψ,𝜃,ϕ) coordinates. If the mode has the property k_∥≈ 0, i.e., m ≈ nq₀, then k_𝜃 is equal to the k_y. The motivation of introducing the constant length r₀∕q₀ in the definition of y is to make k_y ≈ k_𝜃 for a mode with k_∥≈ 0.

Some authors call k_y or k_𝜃 by the name “binormal wavenumber”, which is not an appropriate name in my opinion. Some authors call y the binormal direction, which is also an inappropriate name since neither ∇y nor ∂r∕∂y is along the binormal direction B₀ ×∇ψ.

15.5 Visualization of gridpoints in field aligned coordinate system

In this section, I try to visualize gridpoints in the field aligned coordinates. The directions of the covariant basis vectors of (ψ,𝜃,α) coordinates are as follows:

(336)

(337)

(338)

Here ∂r∕∂ψ|_𝜃,α is a combination of the usual radial and toroidal direction, which needs some clarification. Note that, ϕ is related to α by Eq. (322), i.e.,

(339)

(where the second equality becomes exact if 𝜃 is the straight-field-line poloidal angle defined in Sec. 10.4.), which indicates that, for q′(ψ)≠0 and 𝜃≠0, the usual toroidal angle ϕ is changing when changing ψ and holding 𝜃 and α fixed. Figure 17b shows how the usual toroidal angle ϕ changes when we change ψ and hold 𝜃 and α fixed.

Fig. 17: (a): A 𝜃 contour on ϕ = 0 plane. Here 𝜃 = 9×2π∕63.  This is the direction of ∂r∕∂ψ|_𝜃,ϕ. (b): ψ coordinate lines in (ψ,𝜃,α) coordinates ( ∂r∕∂ψ|_𝜃,α are the tangent lines to these curves) on the isosurface of 𝜃 = 9 × 2π∕63. Here different lines correspond to different values of α. (c): α coordinate lines (∂r∕∂α|_ψ,𝜃 are tangent lines to these curves), which are along the usual toroidal direction . (d): Grid on the isosurface of 𝜃 = 9×2π∕63, where the red lines are α coordinate lines while the blue lines are ψ coordinate lines. Magnetic field from EAST discharge #59954@3.03s.

The relation ϕ ≈ α + q(ψ)𝜃 given by Eq. (339) indicates that the toroidal shift along ∂r∕∂ψ|_α,𝜃 for a radial change form ψ₁ to ψ₂ is given by (q(ψ₂) − q(ψ₂))𝜃, which is larger on 𝜃 isosurface with larger value of 𝜃. An example for this is shown in Fig. 18 for 𝜃 = 19 × 2π∕63, where ∂r∕∂ψ|_α,𝜃 has larger toroidal shift than that in Fig. 17 for 𝜃 = 9 × 2π∕63.

The ∂r∕∂ψ|_𝜃,α curves can be understood from another perspective. Examine a family of magnetic field lines that start from 𝜃 = 0 and ϕ = ϕ₁ but different radial coordinates. These starting points all have the same value of α, which is equal to ϕ₁. When following these field lines to another isosurfce of 𝜃, the intersecting points of these field lines with the 𝜃 isosurface will trace out a ∂r∕∂ψ|_𝜃,α line with α = ϕ₁. Examine another family of magnetic field lines similar to the above but with the starting toroidal angle ϕ = ϕ₂. They will trace out another ∂r∕∂ψ|_𝜃,α line (with α = ϕ₂) on the 𝜃 isosurface. Continue the process, we finally get those curves in Fig. 17b and Fig. 18b.

Radial wavenumber in field-aligned coordinates (ψ,𝜃,α)

For a harmonic in (ψ,𝜃,ϕ) coordinates given by A(ψ,𝜃,ϕ) = exp(ik_ψψ + im𝜃 + inϕ), the radial wave number is k_ψ. Let us calculate the corresponding radial wavenumber k_ψ^⋆ in the new coordinates (ψ,𝜃,α), which is defined by

(340)

where the phase is given by phase = k_ψψ + m𝜃 + nϕ. Then the above expression is written as

Let us examine how many ψ grid points are needed to resolve the ψ dependence in (ψ,𝜃,α) coordinates on the high-field side (𝜃 = π). Assume k_ψ ≈ 0, then k_ψ^⋆ at 𝜃 = π is given by k_ψ^⋆ = nπq′. The corresponding wave-length is given by λ_ψ^⋆ = 2π∕k_ψ^⋆. The grid spacing Δψ should be less than half of this wave-length (sampling theorem). Then the grid number should satisfy that

(343)

where L_ψ is the radial width of the computational domain.

The number of Fourier harmonics that need to be included is given by

(344)

For DIII-D cyclone base case, choose the radial coordinate ψ as r. At the radial location ψ = r₀ = 0.24m, q₀ = 1.4, ŝ₀ = 0.78, then q′₀ = s₀q₀∕r₀ = 4.5m⁻¹. Then N_r^⋆ = n× 0.45 for the radial width L_r = 0.10m.

Figure 19 plots ∂r∕∂ψ|_𝜃,α lines on the 𝜃 = 0,2π isosurfaces, which are chosen to be on the low-field-side midplane. On 𝜃 = 0 surface, ∂r∕∂ψ|_𝜃,α lines are identical to ∂r∕∂ψ|_𝜃,ϕ lines. On 𝜃 = 2π surface, each ∂r∕∂ψ|_𝜃,α line has large ϕ shift. In old version of my code, 𝜃 = 0,2π surfaces are chosen as the 𝜃 cuts (in the new version 𝜃 = [−π,π]). A connection condition for the perturbations is needed between these two surfaces. This connection condition is discussed in Sec. 15.5.

Periodic conditions of physical quantity along 𝜃 and α in field-line-following coordinates (ψ,𝜃,α)

Since (ψ,𝜃,ϕ) and (ψ,𝜃 + 2π,ϕ) correspond to the same spatial point, a real space continuous quantity f expressed in terms of coordinates (ψ,𝜃,ϕ), i.e., f = f(ψ,𝜃,ϕ), must satisfy the following periodic conditions along 𝜃:

(345)

Since (ψ,𝜃,ϕ) and (ψ,𝜃,ϕ + 2π) correspond to the same spatial point, f must satisfy the following periodic conditions along ϕ:

(346)

Since (ψ,𝜃,α) and (ψ,𝜃,α + 2π) correspond to the same spatial point, a real space continuous quantity g expressed in terms of field-line-following coordinates (ψ,𝜃,α), i.e., g = g(ψ,𝜃,α), must satisfy the following periodic condition along α:

(347)

However, generally there is no periodic condition along 𝜃,

(348)

because P₁ = (ψ,𝜃,α) and P₂ = (ψ,𝜃 + 2π,α) are generally not the same spatial point. In fact, equation (322) implies, for point P₁, its toroidal angle ϕ₁ is given by

(349)

while for point P₂, its toroidal angle ϕ₂ is given by

(350)

i.e., ϕ₁ and ϕ₂ are different by 2πq. From this, we know that (ψ,𝜃,α) and (ψ,𝜃 + 2π,α − 2πq) correspond to the same spatial point. Therefore we have the following periodic condition:

(351)

or equivalently

(352)

Numerical implementation of periodic condition (351)

For the fully kinetic ion module of GEM code that I am developing, 𝜃 is chosen in the range [0 : 2π]. The condition (352) imposes the following boundary condition:

(353)

If α is on a grid, α + 2πq is usually not on a grid. Therefore, to get the value of g(ψ,0,α + 2πq), an interpolation of the discrete date over the generalized toroidal angle α (or equivalently ϕ) is needed, as is shown in Fig. 20.

α contours on a magnetic surface

Figure 21 compares a small number of ϕ contours and α contours on a magnetic surface.

As is shown in the left panel of Fig. 21, with ϕ fixed, an 𝜃 curve reaches its starting point when 𝜃 changes from zero to 2π. However, as shown in the right panel of Fig. 21, with α fixed, an 𝜃 curve (i.e. a magnetic field line) does not necessarily reach its starting point when 𝜃 changes from zero to 2π. There is a toroidal shift, 2πq, between the starting point and ending point. Therefore there is generally no periodic condition along 𝜃 since q is not always an integer. A mixed periodic condition involves both 𝜃 and α is given in (351).

In field-line-following coordinates (ψ,𝜃,α), a toroidal harmonic of a physical perturbation can be written as

(354)

where n is the toroidal mode number, m′, which is not necessary an integer, is introduced to describe the variation along a field line. The periodic condition given by Eq. (351) requires that

(355)

To satisfy the above condition, we can choose

(356)

where N is an arbitrary integer, i.e.,

(357)

We are interested in perturbation with a slow variation along the field line direction (i.e., along ∂r∕∂𝜃|_ψ,α) and thus we want the value of m′ to be small. One of the possible small values given by expression (357) is to choose N = −n × NINT((q_max + q_min)∕2), so that m′ is given by

(358)

where NINT is a function that return the nearest integer of its argument, q_max and q_min is the maximal and minimal value of the safety factor in the radial region in which we are interested. [In the past, I choose m′(ψ) = nq − NINT(nq). However, m′(ψ) in this case is not a continuous function of ψ and thus is not physical.] This form is used to set the initial density perturbation in the fully kinetic code I am developing. Note that m′ in Eq. (358) depends on the radial coordinate ψ through q(ψ). Also note that m′ here is different from the poloidal mode number m in (ψ,𝜃,ϕ) coordinate system. It is ready to show that the perturbation given by Eq. (354) with m′∼ 1 and n ≫ 1 has large poloidal mode number m when expressed in (ψ,𝜃,ϕ) coordinates. [Proof: Expression (354) can be written as

(359)

Assume that 𝜃 is the straight-field-line poloidal angle in (ψ,𝜃,ϕ) coordinate system, then δ(ψ,𝜃) = q𝜃 and the above equation is written as

(360)

which indicates the poloidal mode number m in (ψ,𝜃,ϕ) coordinates is given by m = m′−nq. For the case with m′∼ 1 and n ≫ 1, m is much larger than one.]

Since α contours on a magnetic surface are magnetic field lines, they span out the 3D shape of magnetic surface when there are many α contours on a magnetic surface, as is shown by the right-panel of Fig. 22.