Theory of multipole solutions
to the sourceless
Grad-Shafranov equation
in plasma physics

A. Ferreira

R. Goiás, 1021, Jardim Santa Cruz

18700-140  Avaré, São Paulo, Brazil


The rules to write out any one of the linearly independent functions belonging to the infinite set of those in polynomial form that satisfy the sourceless Grad-Shafranov equation as stated in the toroidal-polar coordinate system are established. It is found that a polynomial solution even in the poloidal angle is given by the product of an integral power of the radial coordinate variable by a complete polynomial of equal degree in this same variable with angular-dependent coefficient functions that are linear combinations of a finite number of Chebyshev polynomials in the cosine of the poloidal angle, the numerical coefficients of these being expressed in terms of the binomial numbers of Pascal’s arithmetic triangle. Tables of the ten polynomial solutions of the lowest degrees are provided in variables of the toroidal-polar and of the cylindrical coordinate systems.


In the literature of Plasma Physics devoted to the equilibrium of the toroidal pinch, multipole solutions are generally understood as solutions to the Grad-Shafranov equation with no source terms as represented by the gradient of the plasma pressure and the gradient of half the squared toroidal field function (also called poloidal current function) in flux space, some examples of which have been known since the early days of thermonuclear fusion research [1], [2]. In the present paper we shall designate by multipole fields the magnetic fields in vacuum that are invariant under rotation about a fixed axis in space and whose field lines are entirely contained in the planes passing by the symmetry axis; the fluxes of such fields coincide then with the multipole solutions in the sense this expression is costumarily employed in Plasma Physics. 111It should be noted that an azimuthal magnetic field that falls with the inverse of the distance away from the axis of rotational symmetry and that the presence of a conducting fluid with a uniform pressure distribution in space still do not provide the Grad-Shafranov equation with a source term. We shall exclude such magnetic field and material medium from our definition of a multipole field.

Multipole solutions have found theoretical application in the solution of the free boundary problem of tokamak plasmas, and, by extension, in the calculation of external magnetic field configurations required by such confining devices and design of the system of coils intended to generate them [2]. As it is our aim to show in the article that follows this in the current issue of this journal [3], the utility of the multipole solutions goes beyond this scope, standing their use on the basis of a method of solution to the Grad-Shafranov equation for which the sources are constant in flux space. This multiplicity of utilizations then justifies the study of the solutions to the partial differential equation whose derivation we pass now to outline [4].

For frame of reference we choose the cylindrical coordinate system as illustrated in Fig. 1, which is, among all rotational coordinate systems, the one with the simplest metrical properties222The cyclic order of the unity vectors in the cylindrical system represented in Fig. 1, with the left hand rule being presumed, is . Since this order is the same as that of the unity vectors of same names in the conventional right-handed cylindrical coordinate system, the results for all vector operations in both systems have equal forms. In particular, the terms in the expressions for the curl of a vector have the same respective signals in one and the other systems..

z axis of rotational symmetry
FIG. 1 The left-handed Cartesian coordinate system , the left-handed cylindrical coordinate system and the right-handed toroidal-polar coordinate system .

The multipole fields satisfy the equations:




The solenoidal law, which in the coordinate system of our choice takes the form:


is automatically satisfied if the radial and the axial components of the magnetic field are written in terms of a scalar “stream function” as:


For a rotationally symmetric vector field with components lying on meridian planes the curl admits a nonnull projection only on the azimuthal -direction, which, in the left-handed cylindrical coordinate system represented in Fig. 1, is expressed as:


Replacing and as given by Eqs. (\the@equationgroup@ID) and (\the@equationgroup@ID) in Eq. (id1), we can cast Ampère’s law into the form:


which is, in vector notation, the equation governing the multipole fields. The stream function can be shown [4] to be physically the flux associated with the meridian (or poloidal) magnetic field divided by .

Equation (id1), which we shall call the sourceless Grad-Shafranov equation or the multipole equation, is a partial differential equation of the elliptic type which, referring to the two coordinate systems most widely used in equilibrium studies, the cylindrical and the toroidal-polar (see Fig. 1), admits of infinitely many solutions in the form of polynomials in one coordinate variable with the meaning of length ( or ), the coefficients of which are themselves polynomials, if the first of those two systems is used, also in a linear coordinate variable (), or, if the second one is used, in the cosine of the angular coordinate (). With reference to an equatorial plane in space defined by the coordinate , these solutions can be grouped in even ones and in odd ones, and, inside each parity group, linearly independent solutions can be singled out according to the degree of the polynomial they are in a linear coordinate variable (say). From a mathematical viewpoint it is these independent solutions that we shall recognize as the multipole solutions.

It is the object of the present paper to derive the infinite set of the even multipole solutions to the sourceless Grad-Shafranov equation according to the foregoing definition here given to them. Of the two previously cited coordinate systems, our choice will fall upon the toroidal-polar as the ground for the analytical work, since the developments in this system lead naturally to the representation of the angular dependence of the solutions in terms of orthogonal polynomials that have simple and well known properties.

We shall commence by assigning the multipole solutions the general form of polynomials in the radial coordinate variable with coefficients that are unknown functions of the cosine of the angular coordinate variable.333This is true for solutions of both parities; the odd ones are made so by effect of an overall multiplying sine factor. From the construction of a few examples of the lowest degrees, the precise form of polynomial dependence the linearly independent solutions must bear on the radial coordinate variable is inferred. The ensuing problem to tackle is that of the dependence of the solutions on the angular variable, and accordingly it takes the guise, no longer of a partial differential equation in two variables, but that of ordinary differential equations in the angular one for the coefficient functions of the polynomials whose dependence on the radial coordinate variable has already been established. It is at this stage of the analysis that the representation of the angular dependences of the multipole solutions in terms of orthogonal polynomials appears forcefully as the analytical recourse to be adopted since they are advanced by the very form of the differential equations governing those dependences. Two families of orthogonal polynomials present themselves as natural candidates to play the role of basis for the (finite) expansion of the angular-dependent coefficients: that of the Chebyshev polynomials and that of the associated Legendre polynomials of order unity, of which we have given preference to the former, on account both of the simpler defining properties of its members and of the form we should lend to the solution as suggested by the isolated examples worked out.

The finite set of Chebyshev polynomials that are to enter a multipole solution of a given degree in the radial variable is found by a blend of general arguments and induction. The next task in the progression is then to determine the numerical coefficients by which they must be multiplied in a combination that satisfies the ordinary differential equations for the coefficient functions of the assumed polynomial solution.

At this level the problem will appear as transmuted from that of differential equations, partial and ordinary, into that of a restricted set of difference equations for the numerical coefficients of nonseparable solutions in two variables. Three are these difference equations, two of them ordinary and one partial, the solutions of the ordinary ones serving as boundary conditions to the partial one.

Besides the recurrent construction of the set of numerical coefficients entering the multipole solution of a particular degree they allow for, we show that these difference equations can be given solutions under the form of general expressions by which capability the step-by-step procedure of evaluation of the coefficients can be circumvented. Moreover, since the end values are ipso facto contained in the solution to the partial difference equation, this one comprises also the solutions to the ordinary difference equations and a single expression which attends to all the coefficients belonging to any and all multipole solutions can be written.

Up to the present moment uses seem not to have been ever made of the odd multipole solutions, which however would find application in systems of confinement that would not exhibit up-down symmetry with respect to the equatorial plane of the torus. In the present paper we shall not pursue the full characterization of such odd solutions, limiting ourselves to showing how they arise along a process of systematic construction of polynomial solutions to the multipole equation.



As stated in Section I, we shall adopt the toroidal-polar coordinate system in which to express the multipole equation and to work out its solutions. This system, which is pictorially represented in Fig. 1, is formed by rotating the usual polar system of the plane about a fixed axis in space, which we identify with the -axis. In the planar system, the polar axis is placed perpendicularly to the axis of rotation; the pole is located on the polar axis at a distance from that of rotational symmetry. As the plane of the system is revolved, the polar axis generates a plane, which we refer to as the equatorial plane, perpendicular to the -axis. We shall call meridian plane any plane containing the -axis. The coordinates of a point belonging to the meridian plane determined by the -axis and a given position of the polar axis are then the radius , which measures the distance from the pole to , the polar angle , the angle between the positive direction of the polar axis (the direction of growing distances from the -axis) and the radius, and the azimuthal angle , the angle of rotation given to a reference plane passing by the axis of rotation to bring that plane from an arbitrarily chosen conventional position, to which we assign the azimuthal angle , to coinciding with the meridian plane containing the -axis, the polar axis and the point , such that form a positive set according to the right-hand rule. In this coordinate system Eq. (id1) takes the form:


where ,


is the radial coordinate normalized to and is the flux function normalized to an arbitrary flux . From now on it is Eq. (id1) that we shall refer to as the sourceless Grad-Shafranov equation or the multipole equation.

We shall recognize that Eq. (id1) has infinitely many linearly independent solutions of polynomial form in the variable with coefficients that are polynomials in the variable . Considered as functions of the angle , these independent solutions can be grouped into two sets, of the even ones and of the odd ones; we shall designate the solutions belonging to the first set by the Greek letter and those belonging to the second set by the Greek letter .

Since in Eq. (id1) only derivatives of the unknown function appear, a constant is a solution, and we write:


which we refer to as the even multipole solution of order zero.

A polynomial solution of the first degree in writes in general as:


To determine the forms of the functions and we substitute this expression in the multipole equation. Collecting terms of equal power in and equating the resulting coefficient of each power to zero, we obtain from the zeroth power:


The solution of this differential equation is:


where and are two arbitrary constants. Since our interest is restricted to solutions that are periodic in the variable , we demand to be zero, and thus reduces to:

We next consider the equation resulting from putting the coefficient of the first power in equal to zero. We have:


the second equality coming from the use of equation (2.6a) for in the right hand side of the first equality. This equation can be recognized as an instance of Chebyshev’s differential equation [5], whose solution can be written as:




is the Chebyshev polynomial type I of the first order,


is the Chebyshev polynomial type II of the zeroth order, and and are arbitrary constants.

The terms now remaining in the equation obtained by substituting Eq. (id1) in the multipole equation are all of the second power in . At this last level in the succession of equations for the -dependence of the assumed polynomial solution we meet a condition again regarding , to know:


The solution of this differential equation can be easily obtained and is:


where and are arbitrary constants.

To make equal the two expressions we have for we require that


in Eq. (2.8) and that


in Eq. (2.12). Note that this choice for the free constants and suppresses the terms that are even in in each of the two distinct expressions for .

With and so determined, the expression that results for the solution of the form proposed in Eq. (id1) is:


The constant term is just a multiple of the multipole solution , and the new, independent solution brought about by Eq. (2.16) is:


where we have adopted, to designate it, the notation appropriate to an odd multipole solution in the angular variable, which it is. We shall call this the odd multipole solution of the zeroth order.

We next consider a polynomial solution of the second degree in for the multipole equation:


The procedure of substitution in Eq. (id1) and of equating the coefficients of the powers of to zero leads again, from the lowest power, to Eq. (id1) for , which is solved by Eq. (2.6a) as before. From the equality of the term in to zero we obtain Eq. (2.7) for , the solution of which is given by Eq. (2.8). Since a condition like that of Eq. (2.11) is no longer found when we go to the next power of , we keep both terms in Eq. (2.8) (that is to say, we take ).

The term in in the multipole equation gives rise to the equation:


where the second equality follows from the use of Eq. (2.8) for in the right hand side of the first equality. We are again in the presence of a Chebyshev equation, this time inhomogeneous, the complete solution of which is:


where and are respectively the Chebyshev polynomial type I of the second order, defined by


and the Chebyshev polynomial type II of the first order, defined by


and and are two arbitrary constants.

The constraint that results from equating the term in in the multipole equation to zero closes the set, and, the same as in the previous case of a polynomial solution of the first degree in , takes on the form of a homogeneous differential equation for the -dependent coefficient of the term of the highest power in in the proposed solution, in this way imposing a demand on of its own. We have:


In general, for an assumed polynomial solution of any degree, by following the two opposite leads in the chain of equations for the coefficient functions, one corresponding to the equation associated with the lowest power, and the other corresponding to the equation associated with the highest power of in the multipole equation, we are able to generate two parallel sets of solutions for the -dependent functions. All the consistency of the method of construction of polynomial solutions and in fact the very existence of such solutions hinge on the possibility of making equal these two sets by appropriately choosing the arbitrary constants appearing in both.

A possible manner of bringing , as given by Eq. (2.20), to satisfy Eq. (2.23) is simply by substituting the former into the latter and then choosing for the constants , and the values that solve the ensuing system of algebraic equations. An alternative manner consists in comparing the solution of Eq. (2.23) with the expression for stated in Eq. (2.20). Adopting this second one, we write down the solution of Eq. (2.23):


where and are constants, and then, by equating it to the solution for as given by Eq. (2.20), we are led to determine the constants which appear in both as:


Note that these choices for and have the effect of suppressing the odd functions in that enter the expressions of the two solutions obtained for .

With all these results at hand we are in position to state the polynomial solution of the second degree in for the multipole equation as:


Now, the first two terms on the right hand side are just a combination of the two previously found multipole solutions and , and we recognize the third, new term, as a multiple of the even multipole solution of the first order, which we shall denote by :


If we next consider a polynomial solution of the third degree in , by following the same steps of the procedure we employed to derive the even and odd multipole solutions of lower degrees, we arrive at the odd multipole solution of the first order, which is:




is the Chebyshev polynomial type II of the second order.

The examination of these and of higher order multipole solutions in this way obtained shows that they have definite parity, being either even or odd in the polar angle , the general form of the even ones of order being:


and that of the odd ones of order being:


where the angular functions and should be expected to express themselves naturally as combinations of Chebyshev polynomials type I and type II respectively. The problem then reduces to finding the coefficients of these combinations given the order of the multipole solution. In this paper we shall treat only the case of the solutions that are even with respect to the polar angle .






We start by introducing Eq. (2.32) in Eq. (id1). Collecting terms of equal powers of and equating the resulting coefficient for each power to zero, we obtain three differential equations for the angular functions according to the value of the exponent of . From the term in we obtain an equation for , which is:


and, from the term in , an equation for , namely:


both of which are homogeneous. The coefficient of , for , gives rise to a differential relation involving two angular functions of subscripts differing by unity, which can be written as:


We now proceed by considering each of these equations.

(1) The equation for

Equation (3.1) can be immediately recognized as being Chebyshev’s equation, whose even solution is the Chebyshev polynomial type I of order , [5]. We write thus:


where is a constant.

(2) The homogeneous equation for

Although our ultimate aim is to express as a combination of Chebyshev polynomials, it is nonetheless of interest to consider also other representations which bear a more direct reference to the general pattern the differential equation we found to govern this function obeys. Indeed, the even solution in of Eq. (3.2) finds a concise representation as:


where is the associated Legendre function of degree and order 1 of the first kind [5]. (The second independent solution to Eq. (3.2), which writes as:

being the associated Legendre function of degree and order 1 of the second kind, though it is regular in the whole interval of variation of , , is however odd in the angle and must thence be rejected.) Another representation of that we shall find to be useful is obtained by noting that, upon the transformation:


Eq. (3.2) becomes:


which can be recognized as a form of the hypergeometric equation [5]. The solution of interest for can be written down immediately as:


Since the first argument of the hypergeometric function vanishes for and is a negative integer for , this solution is actually a polynomial of degree in , for which the constant term is missing. (To obtain a convenient representation of the second solution to Eq. (3.2), the transformation of the independent variable proves to be more advantageous than that of Eq. (3.6). By the sole consideration of the form into which the equation is converted when the free variable passes to be , it becomes apparent that one of the independent solution for is [6]:

Because of the sine term multiplying the hypergeometric function, this solution is odd in the variable and therefore to be abandoned.)

We now turn ourselves to our main objective with regard to the angular function , which is to find for it a representation in terms of the Chebyshev polynomials. For this purpose we shall take as starting point the differential equation for itself rather than the explicit definitions we have found by solving it. Equation (3.2) can be concisely written as:


where we have made use of the notation:


being an integer. For the attainment of our aim we have to know which are the effects of the differential operator above introduced as it acts on the Chebyshev polynomials in general.

By having recourse to the frequently quoted relations for the Chebyshev polynomials [5], [7]:


the following properties of the operator can be easily established:



Next, from the considerations on the form of brought forth in the sequel of Eq. (3.8), it is apparent that this function can be suitably expressed as a finite combination of Chebyshev polynomials, restricted these to the ones of the orders smaller than and equal to , and even, to know:


where the ’s are constants.

If now we insert Eq. (3.16) in Eq. (3.9), by making use of the formulae stated in Eqs. (3.14) and (3.15), we are able to transform the differential equation for into a null identity for a combination of Chebyshev polynomials of odd orders, whose coefficients are combinations two by two respectively of the constants ’s. By appeal to the orthogonality property of the Chebyshev polynomials it follows that the coefficient of each polynomial in the identity must in its turn be null, and this gives a recursive pair of formulae for the constants. From the coefficient of we get:


and from the coefficient of we get:


Once is specified, Eqs. (3.17) and (3.18) provide us with the means to evaluate the remaining coefficients that enter the making up of according to Eq. (3.16). We shall defer the accomplishment of this task to the next Section, leaving the matter as it is at this point, and proceed by considering the remaining differential relation for the angular dependences of the even multipole solutions we derived at the beginning of the present Section.

(3) The inhomogeneous equation for ,

For the treatment of Eq. (3.3) it is useful to introduce the differential operator:


which has the property:


as it can be easily proved with the help of Eq. (3.11). Recalling the definition of the operator from Eq. (3.10), we may restate Eq. (3.3) as:


The first use we shall give to Eq. (3.21) is to find the general form the angular functions must obey. For this purpose we shall view it as an inhomogeneous equation for the function , therefore assuming that the right hand side is known. We shall proceed by induction, starting by considering the equation for the index , which is:


the second equality coming from the property of the operator stated in Eq. (3.15) when the object function is as given by Eq. (3.4).

Equation (3.22) is an inhomogeneous Chebyshev differential equation, which admits as complementary solution of even parity in the Chebyshev polynomial type I of order , and whose particular solution should exhibit the same dependence on as the driving term, being therefore proportional to the Chebyshev polynomial of order . We write then for the complete solution:


where and are constants.

We next consider Eq. (3.21) with the index taken to be equal to , namely:


where we have used the form just derived for in the right hand side and again the property of Eq. (3.15) for the differential operator . The even component of the complementary solution to Eq. (3.24) is a multiple of the Chebyshev polynomial type I of order , while the particular solution is a combination of the two Chebyshev polynomials that appear on the right hand side. The complete solution of even parity in the angle then writes as:


where the ’s are constants.

From a consideration of the form of , given by Eq. (3.4), and of those of and , given respectively by Eqs. (3.23) and (3.25), we infer that the form of the function in general must be:


where the ’s are constants. Given that, according to Eq. (2.32), the even multipole solution of order , , depends on angular functions () and since by Eq. (3.26) each of these comprises constants (), the total number of constants needed to specify completely equals  , one of them remaining arbitrary, as it should be proper to the solution of a homogeneous equation of the second order with a definite parity; we shall take this constant, in terms of which all others will be expressed, as . The problem of finding the multipole solution of order to the sourceless Grad-Shafranov equation is then the problem of determining the set of these constants.

To deduce the recursion formulae connecting the coefficients that enter the representation of , we start by rewriting Eq. (3.21) as:

By inserting Eq. (3.26) in the left hand side of this equation and by using the property of Eq. (3.20) for the differential operator , we obtain after a brief calculation:


The reduction of the right hand side of Eq. requires a somewhat lengthier manipulation than that of the left hand side, but it is otherwise straightforward. With the help of the formula stated in Eq. (3.15) we obtain:


The process of equating the left and the right hand sides of Eq. , as given by Eqs. (3.27) and (3.28) respectively, yields two recursion formulae for the coefficients . The first one comes from the terms proportional to on both sides of the equation and is:


(The order in which the values of are written above corresponds to the sequence in which the values of the coefficients are generated by systematic application of the formula, knowing the initial value .)

The second recursion formula stems from the equality of the coefficients multiplying on both sides of Eq. , and is:


The set of difference equations constituted by Eq. (3.17) together with Eq. (3.18), Eq. (3.29) and Eq. (3.30) encompasses all the coefficients that are needed to specify the whole of the angular functions that enter the composition of a multipole solution of a given order and forms thus a complete formulation of the problem of finding the totality of even polynomial solutions to Eq. (id1). In Section IV we shall show how the difference equations can be used to generate the coefficients by a step-by-step calculation procedure; next, partly in Section IV and partly in Section V, we shall derive a solution in closed form for each of the three difference equations, such that a coefficient knowingly belonging to the scope of one of them can be evaluated by the appropriate solving formula from the sole knowledge of its indices; finally, also in Section V, we shall show that the solutions found can be merged into a single formula, from which the whole of the coefficients pertaining to a multipole solution can be obtained in terms of , being enough for that to specify the order of the multipole.


An insight into the mathematical structure of the problem posed by the difference equations to which the original problem of the partial differential equation for the multipole fields was reduced in the last Section can be gained by displaying the Chebyshev coefficients for the solution of a given order in a Cartesian array of columns and rows in which the number attached to a column specifies the first suffix of the coefficients keeping position in that column and the number attached to a row specifies the second suffix of the coefficients belonging to that row. According to Eq. (3.26) the first suffix of the numerical coefficients belonging to the set of those characterized by possessing a specified as the common first suffix is defined to be the subscript of the angular function into whose composition they enter as factors multiplying Chebyshev polynomials, and, by Eq. (2.32), the subscripts of the angular functions that are summoned to participate in the combination that builds up the multipole solution rank from to . Thus the label of the columns in the array must run from to . With regard to the second suffix of the coefficient , which is equal to the order of the Chebyshev polynomial this coefficient multiplies in Eq. (3.26), it must be clear from this equation that its range of variation, considered in its wholeness the constellation of coefficients involved in the construction of the multipole solution , is determined by the angular function , the general expression of which is given by Eq. (3.16). From this we see that the label of the rows in the array must run from to .

Figure 2 exhibits such an array, in which is equal to 4 and the coefficient is taken to be unity. As this example illustrates, not all of the positions are to be filled, since, by virtue of the conventions we have adopted regarding notation, there are no coefficients associated with every pair of indices within the wideness of range of the matrix . In general, the field of indices of the coefficients for a fixed is defined by:


The area enclosing the coefficients resembles that of a triangle, and for this reason we shall refer to this array as the triangle of coefficients. In due time we shall show that it bears a close kinship with the array of binomial coefficients known as Pascal’s triangle.

Each of the difference equations we have derived in Section III applies to the coefficients occupying a different region of the triangle according to a pattern of correspondence which we pass to expound.

4 5 6 7 8
0 1

FIG. 2 The triangle of the coefficients for . The coefficient is taken to be unity. The long arrows on the sides indicate the sequence in which the peripheral coefficients are generated by recursion starting with the coefficient . The short interior arrows intend to signify that the element combines with the element to generate the coefficient .

(1) The equation for the coefficients in the column on the right hand side of the triangle of coefficients

This is Eq. (3.18) together with Eq. (3.17), here reproduced as:


where we have made use of the Kronecker symbol in order to unify the expressions for and for . We find it useful to define a new discrete variable relating to through:


in terms of which Eq. (4.2) is restated as:


Note that, with this transformation, the free variable that appears in the recursion formula for the coefficients is now denoted by the same symbol that, as a suffix, indicates their positions in the column . By choosing any value for (unity, for example, as in the case illustrated by Fig. 2), all the coefficients belonging to the column on the right hand side of the triangle can be evaluated by setting , up to in succession. The arrow on the side of the column in Fig. 2 indicates the sequence in which the positions are filled according to this procedure.

A closed form solution to the homogeneous difference equation for can be easily obtained. Writing down Eq. (4.4) for up to a generic (even) value and then multiplying all the relations so obtained one by the other in succession, we arrive at: