Efficient magnetic fields for supporting toroidal plasmas
The magnetic field that supports tokamak and stellarator plasmas must be produced by coils well separated from the plasma. However, the larger the separation, the more difficult it is to produce a given magnetic field in the plasma region, so plasma configurations should be chosen that can be supported as efficiently as possible by distant coils. The efficiency of an externally-generated magnetic field is a measure of the field’s shaping component magnitude at the plasma compared to the magnitude near the coils; the efficiency of a plasma equilibrium can be measured using the efficiency of the required external shaping field. Counterintuitively, plasma shapes with low curvature and spectral width may have low efficiency, whereas plasma shapes with sharp edges may have high efficiency. Two precise measures of magnetic field efficiency, which correctly identify such differences in difficulty, will be examined. These measures, which can be expressed as matrices, relate the externally-produced normal magnetic field on the plasma surface to the either the normal field or current on a distant control surface. A singular value decomposition (SVD) of either matrix yields an efficiency ordered basis for the magnetic field distributions. Calculations are carried out for both tokamak and stellarator cases. For axisymmetric surfaces with circular cross-section, the SVD is calculated analytically, and the range of poloidal and toroidal mode numbers that can be controlled to a given desired level is determined. If formulated properly, these efficiency measures are independent of the coordinates used to parameterize the surfaces.
In schemes for magnetic plasma confinement, including both tokamaks and stellarators, one must create and control magnetic fields in the plasma region using electromagnetic coils at a distance from the plasma. In tokamaks, not only is it typically beneficial to provide axisymmetric shaping (elongation, triangularity, an X-point, etc), but also to manipulate nonaxisymmetric magnetic field components in order to control error fields and edge localized modes (ELMs). In stellarators, coils shapes are optimized to achieve a plasma shape that has good physics properties. In both tokamaks and stellarators, the distance between the coils and plasma is a critical quantity. Nearby coils cannot be used in a nuclear environment, and a small coil-plasma separation can complicate design and construction, increasing costs. (Indeed, much of the challenge in assembling the W7-X stellarator stemmed from the need to fit many device components in the small space between plasma and coils, “lesson 1” in Klinger et al. (2013).) The distance between coils and plasma can sometimes be increased at the cost of greater current in the coils and/or greater coil curvature, but this increase leads to other engineering challenges such as greater electromagnetic forces in the coils. However, beyond a certain plasma-coil separation, it is typically impossible to find practical coils to support a given plasma shape within a reasonable tolerance Merkel (1987). Thus, magnetic confinement fusion program would greatly benefit from improved and systematic understanding of the efficient control of magnetic fields using distant coils.
Efficiency is a property of both external magnetic fields (those generated by currents outside the plasma) and of plasma equilibria. For an external magnetic field, efficiency is a measure of the field’s shaping component magnitude at the plasma divided by some measure of the magnitude far from the plasma. The efficiency of a plasma equilibrium can then be measured by the efficiency of the external field required to support the plasma (in the sense of a free-boundary equilibrium). These rough definitions will be made precise below.
One might naively expect that the efficiency of a given plasma equilibrium is directly related to the minimum spatial length scale of the plasma boundary shape, with structures of length scale being inefficient to produce using currents displaced from the plasma by a length . However, this reasoning is too simplistic, as the following counterexample shows. A diverted tokamak has a corner at the X-point, meaning the plasma shape contains arbitrarily small length scales. Nonetheless, it has been feasible to produce the X-point shape in many tokamak facilities, and in the design of ITER and many reactor concepts. The divertor coil current required to produce an X-point, which for ITER is the plasma current, has not been prohibitive in practice, so the infinitesmal spatial scales evidently do not imply an infinitesmal efficiency. Another quantity which fails to measure the efficiency of a plasma equilibrium is the spectral width 111 Note that the condition of small curvature is not simply related to the condition of small spectral width of and , and vice versa. The circle , has low curvature but large spectral width when parameterized in terms of the alternative coordinate . Conversely, the shape given by , has small spectral width but infinite curvature. of the boundary shape, which we define loosely as the number of Fourier modes needed to represent the shape via and for some poloidal and toroidal angles . Many Fourier modes are required to represent the boundaries of tokamak plasmas with divertors, while other shapes with very small spectral width, such as shapes with concavities, can be difficult to create with external coils. Furthermore, the spectral width depends on the specific choice of coordinates. Better measures of efficiency are needed.
Stellarator optimizations to date have typically been based upon using a small number of Fourier modes to describe the plasma shape Nührenberg and Zille (1988), This space of shape parameters may exclude some desirable plasma configurations with large spectral width and/or sharp corners, although this is not certain - one can drop the outer surfaces of such an equilibrium to find a smaller similar equilibrium within the parameter space. This space of shape parameters will also include regions where it is difficult to find reasonable coils. To address this issue, recent optimizations have included measures of coil complexity in the overall target functionPomphrey et al. (2001). In such a procedure, measures should be chosen that do not necessarily penalize sharp edges, since such shapes are not necessarily inefficient to produce, and sharp edges may indeed be advantageous for a divertor.
Here, we will examine and compare two superior possible definitions of magnetic field efficiency, the “transfer operator” efficiency
and the “inductance operator” efficiency
where denotes the magnetic field component normal to the surface, and the current potential is the stream function for a divergence-free skin current. The control surface in these definitions is a toroidal surface offset by some distance (uniform or nonuniform) from the plasma. Depending on the application, the control surface may be either a specific surface on which we wish to design coils, or any surface that is relatively far from the plasma.
Our definitions emphasize because this is the ‘shaping component’, the relevant component for controlling the plasma surface shape. For any fixed-boundary MHD equilibrium solution, the same equilibrium could be realized in the realistic scenario of a free boundary if coils could be designed to achieve on the desired plasma boundary Merkel (1987). The linearity of Ampère’s Law means we can write as a sum of contributions from current in the plasma, from fixed coil currents outside the plasma, and from other external shaping currents which are being determined. The term may represent the field due to standard planar toroidal field coils, and/or the net poloidal current on a modular coil winding surface which creates the main toroidal field. The term may represent the field due to saddle coils, nonplanar excursions of modular coils, poloidal shaping coils, etc. (Our division of the externally generated into and corresponds to Merkel’s division into and terms in Merkel (1987).) The problem of designing shaping coils is thus to find shaping currents that achieve , hence our focus on controlling .
The two notions of efficiency (1)-(2) arise from two linear operators, the transfer and inductance operators. The inductance operator relates skin currents on the control surface to on the plasma surface. This approach is closely related to the NESCOIL method of stellarator coil design Merkel (1987). The transfer operator instead relates on the control surface to on the plasma surface. For the part of the magnetic field driven by external coils, for some scalar potential satisfying . A given on the control surface (with normal vector ) represents a Neumann boundary condition, so Laplace’s equation inside the control surface has a unique solution up to a constant, and so a unique on the plasma surface can be determined. The transfer operator encodes this relationship between on the two surfaces. Finite-dimensional inductance and transfer matrices arise from discretizing the quantities of interest on the two surfaces by expanding in a finite number of basis functions.
For both the inductance and transfer matrices, a quantitative measure of magnetic field efficiency can be obtained using a singular value decomposition (SVD). The SVD results in an orthonormal set of right singular vectors, each representing a pattern of or the current potential on the control surface, paired with a left singular vector representing the resulting pattern of on the plasma surface. Each pair of left and right singular vectors is associated with a singular value, which measures the efficiency by which currents or on the control surface drive on the plasma surface. Each singular value, with the associated right and left singular vectors, represents a vacuum solution of throughout the volume enclosed by the control surface, and hence a pattern of everywhere in this volume, not just on the two surfaces. We call each of these spatial patterns of (as well as the associated on the surfaces) a ‘magnetic field distribution’.
Given the transfer or inductance matrix SVD and , the efficiency of a plasma equilibrium shape can be measured using either of two sequences, which we call the efficiency and feasibility sequences. The efficiency sequence is the projection of onto the left singular vectors. It has an exponentially decreasing trend, with rapid decrease for plasma shapes that can be efficiently produced. For example, we will show this decrease is rapid for a diverted tokamak shape, despite the sharp corner at the X point. This sequence is insensitive to the particular choice of control surface, so its rate of decrease measures the intrinsic efficiency of the plasma shape. The feasibility sequence is the efficiency sequence divided by the singular values, representing the amplitudes of or current on the control surface required to achieve on the plasma surface. If this sequence is (exponentially) decreasing, the required amplitudes on the control surface are bounded, so it is feasible to support this plasma shape with coils on this control surface. A plasma shape could be optimized by targeting the rate of decrease of either the efficiency or feasibility sequence.
We examine both the inductance and transfer matrices since each has advantages and disadvantages. The transfer matrix relates a to a , so both the matrix and its singular values are naturally dimensionless. Its singular values have the satisfying property that they approach 1 as the plasma and control surfaces approach each other. In contrast, the inductance matrix is not dimensionless, and its singular values do not have a simple limiting value when the surfaces coincide. However, ultimately we are interested in generating the magnetic field with currents on the control surface, and the connection to the generating current is clearer for the inductance matrix than for the transfer matrix. The inductance matrix is easier and faster to compute, since computation of the transfer matrix requires inverting a dense linear system, and one must carefully integrate over certain singularities (as is common in boundary integral equation methods.) We will find that the inductance and transfer matrices generally give a similar description of external control of the magnetic field. Some long-wavelength structures are found to be more efficient using the transfer matrix than using the inductance matrix. The reason is that magnetic field on the plasma surface is generated by gradients of the current potential rather than by the current potential itself; long-wavelength magnetic field structures can be inefficiently generated even though they propagate efficiently.
Other measures of magnetic field efficiency could be considered besides the ones we consider here. One example is the radius of curvature of the conductors generating the shaping fields, related to curvature of contours of the current potential. Another example is the magnetic energy in the region enclosed by the plasma bounding surface divided by the magnetic energy throughout space.
The present work builds upon the work of several previous authors. SVD techniques have been applied to an inductance matrix in the extension of NESCOIL by Pomphrey et al Pomphrey et al. (2001), and are discussed in Boozer (2000a, b). The concept of the transfer matrix was introduced by Pomphrey in unpublished work, and his numerical calculation appears in figure 2 of Ref Boozer (2015a). The concept of the transfer matrix and the importance of its singular value decomposition were also discussed elsewhere in Ref Boozer (2015a) and in Boozer (2015b).
In the following two sections we introduce preliminary concepts, then review the definitions of the inductance matrix and transfer matrix. Some details of our numerical implementation are given in section IV. In section V we build some intuition for the two matrices and their SVDs by considering the case of axisymmetric surfaces with circular cross-section, since in this limit the SVDs can be computed analytically using a large-aspect-ratio approximation. We also demonstrate that these analytical results agree extremely well with finite-aspect-ratio numerical calculations. Shaped axisymmetric surfaces are examined in section VI, and here we demonstrate (figure 10) that both the inductance and transfer matrix find the X-point divertor shape to be efficient, despite the sharp corner. A numerical example for the W7-X stellarator is presented in section VII. Here we show that our efficiency measures are consistent with the choice of coil winding surface for W7-X, that it is as far from the plasma as possible while still allowing the desired plasma shape to be supported. In section VIII, we prove that for a suitable definition of orthogonality and within discretization error, the SVD is independent of the coordinates chosen to parameterize the surfaces, and also independent of the specific choice of basis functions. We conclude in section IX. An alternative method for visualizing the relationship between the inductance and transfer matrices is presented in the appendix.
Ii Current potential and basis functions
Divergence-free currents on a surface have the form
where denotes a surface current density, is the outward normal vector associated with the surface, and the quantity is known as the current potential Merkel (1987). Contours of the current potential represent streamlines of the current.
In general the current potential may be multiple-valued, because there may be net current that topologically links the surface. For example, the standard toroidal field coils in a tokamak correspond to a current potential that increases as one moves toroidally around the surface, where is the toroidal coordinate. In general we can write
where and are the total currents linking the surface poloidally and toroidally, and are any poloidal and toroidal angle coordinates, and is single-valued. Since and are fixed given the enclosed toroidal flux and the choice of current topology (corresponding to saddle, modular, or helical coils) Merkel (1987, 1988), we can consider the magnetic field driven by to be part of the ‘known’ , whereas generates and is considered ‘unknown’, adopting the same philosophy as in NESCOIL. Hence, for the rest of this paper we will be concerned only with the single-valued part , dropping the subscript to simplify notation.
To describe scalar quantities on surfaces, such as the current potential and the normal component of , we will expand in some set of basis functions . Each surface generally may have its own set of basis functions. It is convenient to require that the basis functions satisfy an orthogonality relation
where represents an area integral over the surface, and is some weight, normalized so that
As we will see below, it can sometimes be useful to consider non-constant . We write the normal component of the magnetic field as a weighted linear combination of the basis functions:
for some components . (Repeated indices imply summation.) Using (5),
The components may be thought of as magnetic fluxes, as they are area integrals of the magnetic field. Similarly, the single-valued part of the current potential may be written
for some components , satisfying
Both and have units of Ampères. Notice that there is some freedom of definition in whether to include the weight in (7) or (8), and in (9) or (10). The choices used in (7)-(10) are made for consistency with Boozer (2015a), where they enable several analogies to be drawn between toroidal magnetic fields and electric circuits.
One natural choice of basis functions is the Fourier basis, consisting of sines and cosines in the poloidal and toroidal angles. This basis requires a specific choice of weight . Observe that area integrals of any quantity may be written
where is the number of identical toroidal periods (e.g. for W7-X.) Also, where
for integers and are valid basis functions satisfying our requirements. (We exclude the constant basis function with . The component associated this basis function must have zero amplitude due to and (8). We can also require the amplitude of the constant basis function to be zero for the expansion of the current potential , since we can shift by an arbitrary constant without altering the physical surface current.) Stellarator symmetry (the 3D generalization of tokamak up-down symmetry) can be imposed conveniently by considering only the sine functions in (II).
As we will discuss in section VIII, the Fourier basis (II) has an unsatisfying property that some important quantities we will compute depend on the choice of coordinates and used to parameterize the surfaces. For example, the plasma surface could be parameterized by the coordinates of the widely-used VMEC code Hirshman and Whitson (1983) (in which magnetic field lines are not straight), or any of the popular straight-field-line coordinates. In section VIII we will show this coordinate-dependence is a result of the fact that the weight (13) depends on the choice of coordinates. Therefore, it is advantageous to consider a different set of basis functions orthogonal under a weight that is constant on the surface, and therefore independent of the choice of coordinates. From (6), such a weight must be where is the area of the surface. There are multiple ways to generate suitable basis functions. One approach is to use the functions
for integers and . This is the basis set used for the rest of this paper, except where noted otherwise. A different set of basis functions orthogonal under the same constant weight is discussed in section VIII.
Iii Inductance and transfer matrices
Consider a toroidal surface on which a skin current flows. We can compute the resulting magnetic field at any point using the Biot-Savart Law
Integrating by parts in and to remove the derivatives on , after some algebra we obtain
where . This magnetic field is precisely that generated by magnetic dipoles of density covering the surface and oriented normal to it.
Let us call the surface considered so far the “control” surface outside the plasma, and consider a second toroidal surface inside this control surface, representing the boundary of a plasma. We evaluate the component of (18) normal to this inner surface, giving a linear relation between and the magnetic field normal to the plasma surface, denoted . The operator in this linear relation is the inductance operator. A discretized linear relation is obtained by applying the operation , where are the basis functions on the plasma surface. Recalling (8), (9), and (10), and letting denote the basis functions on the control surface, we obtain
(or in alternate notation ) where
Quantities without primes refer to the plasma surface. We call the (mutual) inductance matrix between the two surfaces, since it relates currents to fluxes.
Next, we define the transfer operator and matrix. As noted in the introduction, given on the control surface, denoted , Laplace’s equation gives a unique solution for on the plasma surface, denoted . Due to the linearity of Laplace’s equation, is a linear functional of : for some linear operator . We call the transfer operator. A discrete version of this linear relation also holds: the vector of fluxes on the plasma surface must be a linear function of the vector of fluxes on the control surface . The transfer matrix is precisely this linear relationship:
For any real matrix , including the inductance and transfer matrix, one can find a singular value decomposition
where and are orthogonal ( and ), is diagonal, and the diagonal elements of are real and non-negative. These diagonal elements are the singular values, and the columns of and are the left singular vectors and right singular vectors respectively.
The SVD provides an intuitive way to understand the action of the matrix on any vector . First, is decomposed in the orthonormal basis provided by the columns of . Each component is then scaled by the corresponding singular value. The results then give the amplitudes of the response in the basis provided by the columns of . Hence, the singular values of the transfer and inductance matrix represent the efficiency by which quantities on the control surface (either or respectively) propagate to the plasma surface. The left- and right-singular values associated with each singular value represent the modes of or which propagate with the given efficiency. The singular vectors of the transfer and inductance matrices give our efficiency-ordered magnetic field distributions.
Let and denote the flux vectors associated with the normal magnetic fields on the plasma surface and defined in the introduction. If the number of basis functions retained on the plasma and control surfaces is the same, then the singular-valued part of the current potential required to achieve total can be computed in principle from . We say ‘in principle’ because is very poorly conditioned, for a large current in a sufficiently small loop on the control surface will cause little change to the plasma surface . If the number of basis functions on the surfaces is not the same, one can still apply aforementioned relation if is interpreted as the Moore-Penrose pseudoinverse, which is equivalent to solving for in a least-squares sense. This procedure for computing a current potential is essentially the one used in the NESCOIL code. Using the SVD of we can write
To understand why large shaping currents result from this procedure for some plasma shapes but not for other shapes, we will consider various subsets of the products on the right-hand side of (23). Let denote the rightmost product, , which is the projection of the externally generated plasma surface onto the left singular vectors (which represent distributions on the plasma surface.) We call the elements of the efficiency sequence. The elements of this sequence are relatively insensitive to the distance between the plasma and control surfaces, since this separation mainly affects the singular values rather than the unitary matrix . (Recall the eigenvalues of a unitary matrix all have unit norm.) Including the next term in the product (23), we define the feasibility sequence , which is the efficiency sequence divided by the singular values. Roughly speaking, the sum of the feasibility sequence determines the magnitude of the required shaping currents , since and is unitary. An increasing feasibility sequence has a diverging sum, so the magnitude of the required shaping currents diverges, hence it is likely infeasible to find suitable coils. Conversely, a rapidly decreasing feasibility sequence means the amplitudes of shaping currents are small, hence it should be feasible to find suitable coils. As the feasibility sequence depends on the singular values, it will be sensitive to the separation between plasma and control surfaces. Efficiency and feasibility sequences can be defined for the transfer matrix just as for the inductance matrix.
Iv Numerical solution
The inductance matrix can be computed from the definition (20), discretizing the two surface integrals using uniform grids in and on each surface. A convenient expression for numerical evaluation of the inductance matrix is
where the matrix element is given by basis function on the plasma surface evaluated at a grid point indexed by , an analogous definition holds for (representing the basis functions on the control surface),
and denote the spacing of grid points, and are evaluated at a grid point on the plasma surface indexed by , and and are evaluated at a grid point on the control surface indexed by . Corresponding to (11), the and sums represent integrals over the interval (i.e. over one of the identical toroidal periods) and indexes the toroidal periods. Due to symmetry the corresponding sum is reduced to the factor of in (25). Computation times are typically dominated by the assembly of and by the matrix multiplications in (24). Calculations can be significantly accelerated with little programming effort by using OpenMP parallelism for the former and a multi-threaded BLAS subroutine for the latter.
The basic resolution parameters of the numerical calculation are the number of grid points and on each surface, determining both dimensions of and determining the number of rows of and , and the maximum poloidal and toroidal mode numbers and for the set of basis functions on each surface, determining the number of columns of and and determining both dimensions of . Typically it is appropriate to choose and . In principle, different numbers of grid points could be used on the two surfaces, and/or different numbers of basis functions could be used on the two surfaces. However for all results presented here, the resolution parameters are identical on the two surfaces.
Now consider computation of the transfer matrix, which amounts to solving Laplace’s equation in the region between the plasma and control surfaces. This problem is a classic case in which boundary integral equation methods excel Hess (1975); Jaswon and Symm (1977); Merkel (1986, 1988); Becker (1992), since these methods avoid the need to discretize the volume inside the control surface; only the surfaces need to be discretized. Within boundary integral equation methods, there are many options. Methods can be divided in two classes, those which use Green’s theorem as in Merkel (1986), and those which use a surface source distribution Hess (1975) in which we write a surface integral
Here, is some surface distribution which is solved for in terms of the Neumann boundary condition on the control surface, and is a solution of Laplace’s equation with a singularity at , such as a monopole or dipole. The sources are outside the control volume (by a finite or infinitesimal distance) and are not physical, introduced to span the space of possible solutions to Laplace’s equation inside the control surface, so can be the monopole function even for this magnetic field calculation. For surface source distribution methods, one can choose the surface on which the sources are located (i.e. the integration surface in (26)) to either be infinitesimally outside the control surface, or it can be a distinct third surface which is outside the control surface by a finite distance. In the latter case, the transfer matrix must be independent of the specific choice of outermost surface if the calculation is well resolved numerically, since the transfer matrix is originally defined without reference to this third surface. For good numerical convergence, this outer surface must be far from the control surface compared to the spacing of integration points, but not so far that the coupling between outer and control surfaces is poorly conditioned. The 3-surface method has the disadvantage that selecting a third surface shape which is an appropriate distance from the control surface sometimes requires manual adjustments; however the 2-surface source distribution method and Green’s theorem method have the disadvantage that one must carefully integrate over a singularity at the point , as in (26). These singularities may be handled using specialized quadrature rules, or using the remarkable regularization technique of Merkel Merkel (1986, 1988) in which one adds and subtracts a singular function that can be integrated analytically. (Note that the appendix of Merkel (1986) contains many typographic errors. Readers should refer instead to Merkel (1988) which contains the correct equations.)
For any of these schemes, the transfer matrix is computed using a similar sequence of steps. First, a dense matrix is assembled which relates either (in the Green’s theorem method) or (in the surface source method) on the outermost surface to on the control surface. The inverse of this matrix is left-multiplied by another matrix that relates or on the outermost surface to on the plasma surface. For example, the method described in Boozer (2015a) for computing the transfer matrix is a 3-surface source distribution method with dipole sources. One computes the outer-to-plasma-surface and outer-to-control-surface inductance matrices and , and then the transfer matrix follows from . This method amounts to spanning the set of solutions to Laplace’s equation using dipoles on the outer surface, where the dipole amplitudes are chosen to satisfy the Neumann boundary condition on the control surface.
For results shown here, we use two different surface source distribution methods, both the 2-surface scheme with monopole source distribution and Merkel’s regularization, and the 3-surface scheme with dipole source distribution, and we have verified the two approaches give identical results as a check of correctness. For the 2-surface method, we use the same uniform grid for and . For every calculation shown in the figures below, computation of the transfer and inductance matrices and their SVDs took in total 160 s on one node (24 cores) of the NERSC Edison computer, and computations could be accelerated with further code optimization if desired.
V Circular axisymmetric surfaces
To gain familiarity with the transfer and inductance matrices and their SVDs, let us first consider the case in which both the plasma and control surfaces are axisymmetric, though the magnetic field distributions we consider may be nonaxisymmetric. We begin by analyzing the cylindrical limit, corresponding to large aspect ratio and circular cross section of the surfaces. In this limit, we can compute the singular values and vectors of the transfer and inductance matrices analytically, because the singular vectors correspond to separable solutions of Laplace’s equation for the scalar magnetic potential in cylindrical geometry. This correspondence can be understood as follows. Laplace’s equation becomes
where constant. Some of the separable solutions are
for integers and , where is the modified Bessel function, and and are arbitrary constants. There are also “axisymmetric” solutions
Forming and noting the weight is constant, these fields represent orthogonal flux vectors on both the plasma and control surfaces. Thus, a singular value decomposition for the transfer matrix can be constructed as follows: the columns of represent the (normalized) distributions (V)-(V) for various evaluated on the plasma surface, the columns of represent the (normalized) distributions (V)-(V) for various evaluated on the control surface, and the singular values are the ratios of fluxes between the plasma and control surfaces. Since the SVD is unique up to the signs of the singular vectors under suitable assumptions (singular values in decreasing order, matrix is square, etc.), then the SVD constructed above must be “the” SVD of the transfer matrix. The singular values are thus
for the fields (V), where and are the minor radii of the plasma and control surfaces, and
for the fields (V). One can see from these formulae that in the limit that the plasma and control surfaces approach each other, , the singular values become 1.
If is not as large as the aspect ratio of the surfaces, the Bessel functions in (30) can be expanded for small argument, giving for , or for . In the latter case, it is interesting that the singular value becomes independent of the physical wavelength , so the distributions are actually less efficient (i.e. have singular values that are smaller by ) than the shorter wavelength distributions. The distributions are comparable in efficiency to the distributions.
A similar analysis can be done in this cylindrical limit to compute the singular values of the inductance matrix. In this case, we must also consider the Laplace equation solutions which extend outside the control surface to infinity, replacing in (V) and (the modified Bessel function of the second kind) in (V). Continuity of at the control surface is imposed, as is a jump condition on the tangential magnetic field associated with the surface current. Supposing the current potential has the form
for some , we thereby obtain the radial magnetic field inside the control surface to be
when , and
when . Computing the magnetic flux on the plasma surface and the current on the control surface, the singular values then follow from their ratio:
when , and
for . These formulae, unlike those for the transfer matrix, remain dependent on , , and the aspect ratio even when the plasma and control surfaces coincide, .
Some insight into (V) can be gained by expanding the Bessel functions for small argument, which is valid when is below the aspect ratio. In this case, when is nonzero, we obtain again (36), indicating the singular values are nearly independent of . However, when , the expansion of (V) yields instead
This expression indicates the singular values actually grow like , so the dominant singular vectors correspond to short wavelengths rather than long wavelengths. This behavior is one instance of a phenomenon mentioned in the introduction: the inductance matrix sometimes treats long-wavelength field distributions as inefficient, because the field is driven by the gradient of the current potential, and this gradient is small for long wavelengths. At sufficiently large however the approximation leading to (37) breaks down, and the singular values begin to decrease with .
For the geometries in figure 1, the singular values of both transfer and inductance matrices computed using these approximate analytic formulae are shown in figure 2. The singular values of the inductance matrix have been normalized to the characteristic value so the largest singular values are . Overall, the distributions of singular values for the transfer and inductance matrices are quite similar, aside from two factors. First, the behavior is different, as already discussed, with the singular values of the transfer matrix decreasing monotonically with whereas those of the inductance matrix peak at intermediate . Second, the singular values of the transfer matrix generally decrease faster with and compared to the singular values of the inductance matrix, and so a different color scale is used for the two matrices in the figure. The physical reason for this different rate of decrease is the same effect described in the previous paragraph.
One practical consequence of the analysis in this section is that we can analytically determine the toroidal mode numbers which can be controlled with comparable efficiency to given poloidal mode numbers. In other words, given some poloidal mode number , what is the toroidal mode number such that the and modes have comparable efficiency, as determined by the singular values? To proceed, we compute from (30), expanding the Bessel functions for large argument, . Equating the result to from (31), we find
In figure 2 we show ellipses in -space, where the extent in is determined by (38) for several choices of . It can be seen that (38) indeed gives quite an accurate estimate for the range of and which have comparable efficiency. While (38) was derived based on the transfer matrix singular values, it does not appear possible to form such an explicit expression for the inductance matrix singular values due to their more complicated dependence on and . However the transfer and inductance matrix contours in figure 2 are quite similar.
Next, we compare these analytic results to direct numerical calculations which fully account for toroidal effects. First considering only axisymmetric basis functions for simplicity, results are shown in figure 3 for three values of the minor radius of the plasma surface, taking the aspect ratio of the control surface to be . The agreement with (31) and (36) is extremely close. The distributions given by the singular vectors are found to correspond nearly exactly to Fourier modes of increasing , as expected, coming in pairs with or symmetry. The factor of in (36) which is not present in (31) causes the singular values of the inductance matrix to decrease less rapidly than those of the transfer matrix, particularly at low .
Next, we consider the set of basis functions given by various for , with results shown in figure 4. Excellent agreement is seen with the analytical results (30) and (V). As expected, the distributions given by the singular vectors are found to correspond nearly exactly to Fourier modes of increasing for the transfer matrix, and to Fourier modes of decreasing for the inductance matrix. For these figures, numerical parameters used were and .
Figure 5 mirrors figure 2, but showing data from the fully toroidal numerical calculation. In this calculation each singular vector contains a mixture of and values, so the figure shows the singular value associated with the left singular vector that has maximal amplitude of each given . Figures 2 and 5 are not exactly identical, due to toroidal effects which are not accounted for in the analytic expressions plotted in the former figure. However, the agreement is generally quite good, demonstrating that the analytic expressions (30), (31), (V), and (36) are remarkably accurate even at realistic values of aspect ratio. The predicted non-monotonic dependence of the singular values of the inductance matrix can be clearly seen.
Equations (30)-(31) and (V)-(36) may actually be useful in applications that do not directly involve the transfer or inductance matrices, since they provide an efficiency ordering to pairs of poloidal and toroidal mode numbers . For example, considering the Dommaschk potentials Dommaschk (1986) for vacuum magnetic fields, the potentials for various numbers can be sorted by to order them by efficiency.
Finally, figure 6 illustrates the behavior of the transfer matrix SVD at finite aspect ratio. Figure 6.a shows the geometry, with the plasma aspect ratio taken to be 3, and chosen to be 1.7, representing a reasonable distance from the plasma at which control coils might be placed. Considering both sine and cosine symmetries, figure 6.b shows there are approximately 500 distributions of with singular values above 0.01. For comparison, the results of (30)-(31) are also plotted, although the large-aspect-ratio approximation used to derive these analytic results is not very well satisfied. Despite this questionable approximation, the analytic method does a good job of predicting the distribution of the ‘true’ (numerical) singular values. Numerical parameters used were , , , and . We verified that none of the points in figure 6 moved visibly when any of the resolution parameters was doubled.
Red crosses in figure 6.b show the numerically computed singular values if only the axisymmetric distributions are retained. These singular values correspond to translation, elongation, triangularity, squareness, etc., and they come in pairs (barely visible in the plot) corresponding to sine and cosine symmetry. It can be clearly seen that to control to any given tolerance, there are many more nonaxisymmetric degrees of freedom than there are axisymmetric degrees of freedom.
Vi Shaped axisymmetric surfaces
Now that the singular vectors and values of the inductance and transfer matrices are known for the case of circular surfaces, we proceed to consider plasma and control surfaces which are axisymmetric but non-circular. At this point we will be able to demonstrate one of our main results: both the inductance and transfer matrix methods can correctly identify the X-point divertor shape as being efficient, despite its sharp corner. We will also give examples of plasma shapes that have low curvature and low spectral width, yet are inefficient to produce.
The four shapes we will consider are shown in figure 7. The ‘divertor’ shape is taken from a numerical Grad-Shafranov solution for the ITER 15 MA baseline scenario. Since our numerical methods presume differentiable surface shapes, the surface at normalized poloidal flux = 0.995 rather than 1 is chosen so the corner at the X-point is very slightly rounded off. The other three shapes are chosen to have small curvature and small spectral width, in the sense that all can be represented by Fourier series with 4 or fewer terms:
(For contrast, we retain 100 terms in the Fourier series for the divertor shape.) The ‘ellipse’ shape is chosen to have the same width, height, and median major radius as the divertor shape. The ‘peanut’ and ‘H’ shapes are chosen just large enough to encircle the divertor shape, so any difficulty generating these two shapes cannot be due to larger distance from a common control surface or coils. For the ellipse, peanut, and H shapes, we have generated MHD equilibrium solutions just as for the divertor shape, this time using the VMEC code Hirshman and Whitson (1983) run in fixed-boundary mode. To emphasize that the divertor shape has a significantly sharper corner than the other shapes, figure 8 shows the curvature (i.e. inverse radius of curvature) of the poloidal cross-section
where primes indicate , and is independent of the poloidal coordinate used.
Intuitively, we expect the ellipse and divertor shapes to be efficient and the peanut and H shapes to be inefficient in some sense, for plasmas of roughly the ellipse and divertor shapes are common in experiments whereas the peanut and H shapes are uncommon. One calculation supporting this notion of efficiency is presented in figure 9. Here, we first compute , the magnetic field normal to the desired plasma surface driven by plasma current inside it, using the virtual casing principle Shafranov and Zakharov (1972); Drevlak et al. (2005). Next, we solve the linear least-squares problem to determine the currents in the actual ITER shaping coils (6 poloidal field coils and 6 central solenoid coils) which minimize , where (in this case ) and is the contribution from the shaping coils. Contours are then plotted of the total poloidal flux , summing the plasma and external contributions. It is clear that the ellipse and divertor shapes can be produced to high accuracy by the ITER coilset, whereas the least-squares approximations to the peanut and H shapes lack the desired concave regions. This failure comes despite the fact, mentioned above, that the peanut and H target shapes are everywhere closer to the coils compared to the divertor shape. (The maximum coil currents for these least-squares fits are 29 MA, 13 MA, 45 MA, and 19 MA for the ellipse, peanut, H, and divertor shapes respectively.) If additional coils are included between the actual ITER coils at similar distances from the plasma, the concave plasma shapes can be better approximated, but only with an unrealistic number of coils and large coil currents.
For calculations using the inductance and transfer matrices, we will explore the effect of the choice of control surface by considering two options. One option, the ‘ITER control surface’, is obtained by fitting the locations of the ITER poloidal field and central solenoid coils with a 4-term Fourier series as in (39). This control surface thus represents a realistic surface on which shaping coils may be placed in an experiment. For contrast, the other option is an ‘offset’ control surface obtained by expanding the plasma surface by a uniform 1.5 m distance. This choice will enable us to demonstrate that tailoring the control surface to the plasma shape has little effect on our results. The control surfaces are shown in figure 7. For these calculations we consider only axisymmetric basis functions and singular vectors, and we allow up-down asymmetry.
Our main result is figure 10, which shows the efficiency sequences for the four shapes, i.e. the projection of onto the left singular vectors. The toroidal field can be ignored as it gives no , and we can thus take . The four equilibria have the same total plasma current (15 MA) and hence comparable poloidal field, so it is meaningful to compare the absolute magnitudes in figure 10 in addition to the rates of decrease. The figure shows that to cancel to a relative accuracy of , only the first 20-25 efficiency-ordered distributions are required for the ellipse and divertor shapes. The divertor shape is not significantly harder to produce than the ellipse in this regard. In contrast, the ‘long tails’ for the peanut and H shapes in figure 10 indicate that many more distributions are required to create these shapes. For the peanut shape, even with the first 80 efficiency-ordered distributions, it is not possible to cancel to a relative accuracy of . The trends are nearly identical for the transfer and inductance matrices.
Comparing the top row of subplots in figure 10 to the bottom row, it can be seen that changing the control surface may change the precise values of the efficiency sequence, but the overall rate of decrease is not significantly changed. We have also repeated the analysis with other control surfaces, including surfaces with elliptical or circular cross-section, and the results are always quite similar to figure 10. This robustness is expected, for as pointed out previously, the efficiency sequences are independent of the singular values, and hence largely independent of the distance between plasma and control surfaces. Thus, the rates of decrease of the data in figure 10 are a good measure of the intrinsic efficiency of producing a plasma shape.
To better understand these results, figure 11 displays additional details of the calculation for the divertor plasma shape with ITER control surface. Figure 11.a shows the field lines of the part of the magnetic field driven by the plasma current. It can be seen that this magnetic field points into the plasma surface on one side of the X point and out of the plasma surface on the other side of the X point. Hence, has a sharp transition from one sign to the other in this region (figure 11.b). We use high numerical resolution () so this transition is well resolved numerically. Figure 11.c then displays the first few right singular vectors of the transfer matrix, corresponding to the magnetic field component normal to the control surface for the most efficient magnetic field distributions. These functions resemble , , , , etc., corresponding to vacuum fields with slow spatial variation, as expected from section V. The first right singular vectors of the inductance matrix (11.d) have slightly faster spatial variation, consistent with the trend we observed in section V. The components of these distributions normal to the plasma surface are plotted in figures 11.e-f, corresponding to the left singular vectors. These singular vectors have fine structure near corresponding to the X point, for exactly the same reason does: a with slow variation is projected normal to a surface with a sharp corner. Hence, the basis of left singular vectors is able to resolve fine structure near the X point, precisely where has fine structure.
Next, in figure 12 we examine the singular values for the four shapes and two control surfaces. The overall pattern of singular values is extremely similar to the case of circular surfaces, figure 3. The singular values for circular surfaces behaved like (in pairs with sin or cos symmetry), and we can interpret the singular values for shaped surfaces as following the same pattern for an effective ratio of surface radii . Particularly for the ellipse and divertor shapes, for which the 1.5 m offset control surface is everywhere closer to the plasma than the ITER control surface, the rate of decrease in the singular values in figure 12 is significantly faster for the more distant surface, as expected. For the ITER control surface, the singular values for the peanut and H shapes are larger than for the ellipse and divertor shapes, which merely reflects the fact that the ellipse and divertor shapes are smaller. For a given effective , the singular values are not very sensitive to the details of the surface shapes. Hence the singular values themselves are not an ideal target for optimization. The projections shown in figure 10 are more sensitive to the plasma shape, making them a better target for optimization.
Finally, we divide the efficiency sequences (figure 10) by the singular values (figure 12) to get the feasibility sequences in figure 13. Depending on whether the efficiency sequences or singular values decrease more rapidly, the feasibility sequences may be (exponentially) increasing or decreasing. Indeed, the sequence for the ellipse shape is decreasing in the figure, while the sequences for the peanut and H shapes are increasing. The sequences for the divertor shape have a significantly lower slope than the sequences for the peanut and H shapes, reflecting the fact that the divertor shape is easier to generate despite its larger curvature. (We only show the first 24 sequence elements for the divertor shape, since after this point the projection becomes polluted by discretization error from the numerical ITER Grad-Shafranov solution we have available. )
We have repeated the analysis of figures 10 and 13 taking to be generated by an infinitesimally thin ring of current at the magnetic axis, and find only minor changes to the results. The details of the current distribution inside the surface are unlikely to have a strong effect on generally since the plasma current is usually concentrated near the magnetic axis, relatively far from the boundary where is evaluated.
Vii Nonaxisymmetric surfaces
Now that some understanding of the inductance and transfer matrix SVD properties has been developed in the previous sections, let us present one example from the important but more complicated case of nonaxisymmetric surfaces. Some features from the previous axisymmetric cases will persist. Application of the inductance and transfer matrix SVD for nonaxisymmetric shape optimization will be considered in future work.
Figures 14-18 present the SVD of the transfer and inductance matrices for a W7-X case. The plasma surface is taken from a fixed-boundary VMEC equilibrium () from prior to the selection of the W7-X coil shapes, so there is no ripple due to the discrete coils of the actual experiment. We consider two control surfaces. The ‘real’ control surface is chosen to be a surface on which the W7-X coils lie. The distance between the plasma and this control surfaces is not exactly uniform. (Not surprisingly, based on the difficulty of generating concave shapes as discussed in the previous section, the designers of W7-X chose to make the coil surface closer to the plasma in the concave regions.) We also consider an ‘expanded’ control surface, obtained by adding to and to of the real control surface. For computing the inductance and transfer matrices, in light of the 5-fold toroidal periodicity of W7-X, we enforce 5-fold toroidal periodicity in the basis functions. We also force to have stellarator symmetry by including only stellarator-symmetric basis functions in the calculation.
Figure 14.c shows the singular values of the transfer and inductance matrices. Just as in the circular axisymmetric case of section V, the singular values decrease roughly exponentially, and the rate of decrease is a bit faster for the transfer matrix. The rate of decrease is also faster for the expanded control surface since it is farther from the plasma.
Figures 15-16 show the left and right singular vectors of the transfer matrix for the real control surface, transformed to -space. Comparing these two figures, the patterns on the plasma and control surfaces are similar. As one would expect, there is a general trend for the spatial frequencies to become higher as one proceeds through the sequence of singular vectors.
Figures 17-18 show the left and right singular vectors of the inductance matrix for the real control surface, again transformed to -space. The spatial frequencies of these singular vectors are noticeably higher than those of the transfer matrix singular vectors, particularly on the control surface (figure 18). This phenomenon is yet another instance of the inductance matrix favoring shorter wavelengths, due to the gradient in the relation between current potential and magnetic field.
Finally, in figure 19 we examine the difficulty of achieving a total on the W7-X plasma surface, showing the efficiency sequences given by projecting onto the left singular vectors of the transfer and inductance matrices. Since the surface shapes are nonaxisymmetric, the net poloidal current on the control surface (which we take to flow at constant , as in NESCOIL) produces a nonzero . (The trends in figure 19 are unchanged if is considered to arise instead from an axisymmetric toroidal magnetic field .) The plasma current in W7-X is small, and in this example is negligible compared to . Figures 19.a-b illustrate that the general rates of decrease in the efficiency sequences are the same for the two control surfaces, as expected. To achieve a total on the plasma surface by adding non-planarity to modular coils on the control surface, the required magnitudes of the or current distributions on the control surface are obtained by dividing the figure 19.a-b sequences by the singular values, yielding the feasibility sequences. Results are shown in figures 19.c-d. For the real control surface, the rate of decrease in figure 19.a-b is just barely faster than the rate of decrease in the singular values, so the resulting sequences in figures 19.c-d show a trend that is nearly flat but very slightly decreasing. This trend is in fact expected, for it indicates that the W7-X coils are just close enough to the plasma to achieve the necessary control of , consistent with the aim of the W7-X designers to push the coils as far from the plasma as possible. If coils were further from the plasma so as to lie on the expanded control surface, the singular values of the inductance and transfer matrices decrease more rapidly (figure 14.c,) yielding feasibility sequences in figures 19.c-d that are generally increasing. The amplitude of or current needed on the control surface is thus unbounded, and one would effectively need to control an infinite number of the distributions on the control surface. Thus, the nearly flat trends in figures 19.c-d for the real control surface are consistent with the choice of this winding surface for the actual W7-X coils: it is as far from the plasma as possible, at limit of where it is practical to achieve total on the desired plasma surface.
Viii Universality and uniqueness: choice of angles and basis functions
We now analyze several technical issues related to the uniqueness of the transfer and inductance matrices: are the SVDs of these matrices independent of the specific basis functions chosen? And are the singular values and vectors independent of the choice of poloidal and toroidal angles and ? For example, one could parameterize the plasma surface using the coordinates of the VMEC code (in which magnetic field lines are not straight), or the straight-field-line coordinates associated with the normal cylindrical toroidal angle, or Boozer coordinates, etc. We will prove and demonstrate numerically that the SVDs are indeed independent of these choices of basis functions and angles, under reasonable assumptions. Specifically, we loosely prove the following two results for both the inductance and transfer matrices: (1) Two different sets of basis functions which are orthogonal under the same weight yield the same numerically converged singular values, and the associated singular vectors represent the same patterns in real space; and (2) as long as the weight is specified in a coordinate-independent manner, changing the choice of angles does not alter the numerically converged singular values or real-space patterns represented by the singular vectors.
viii.1 Change of basis
To carry out the rough proofs, we begin by considering a single surface with two different sets of basis functions, and . Let the two basis sets be associated with weights and which may or may not be equal. We introduce a matrix with components which relates the bases:
(For finite basis sets the basis functions in one set may not be exact linear combinations of the functions in the other set, but once a sufficient number of basis functions are included in each set to resolve a given number of singular vectors to a given level, (41) holds in an approximate sense. Due to this issue, our conclusions will apply only to the singular vectors and values that are well resolved numerically, rather than to all of them.) Applying the operation to (41) gives
We now form
where we have applied (42) to the first factor of and then used (41). If , then the right-hand side of (VIII.1) becomes , so . We thus arrive at the important conclusion that if the two bases have the same weight, then is an orthogonal matrix.
viii.2 Uniqueness of inductance matrix
Now suppose we have two inductance matrices, both relating the current potential on the control surface to on the plasma surface, but using two different sets of basis functions on the plasma surface, and :
where the integrals are over the two surfaces, superscripts indicate = plasma surface and = control surface, and
Now if we have the SVD for basis 1,
Let us now suppose the two sets of basis functions use the same weight so the matrix is orthogonal, as proved above. Then since both and are orthogonal, so is . Hence the SVD for basis 2 is
We have thus proved that the singular values and the right singular vectors are unchanged under the change of basis. The th left singular vector in basis 1 represents the following