Vacuum Solution for Solov'ev's Equilibrium Configuration in Tokamaks

In this work, we have revisited the Solov'ev analytical solution for a tokamak equilibrium. We point out that the vacuum solution in the Solov'ev formulation is inapplicable, since a distributed current density is assumed to fill the vacuum. The realistic vacuum should be current-free. To amend this vacuum solution problem, we use Green's function method to compute the plasma current contribution, together with a homogeneous solution to the Grad-Shafranov equation, to construct the full solution. Matching with the Solov'ev solution on the last closed flux surface is performed to determine the homogeneous solution. The total solution is then extended into the vacuum region to get a realistic vacuum solution. We find that the actual vacuum solution is different from the Solov'ev solution in the vacuum region, especially the X-point structure. The X-point obtained at the last closed flux surface is not like the letter"X", and the expanded angle in the vacuum is larger than corresponding angle in the plasma at the null point. The results are important for understanding the X-point and separatrix structure. At the end of the paper, we have extended the classic Solovev's configuration to an ITER-like configuration, and obtained the full solution.


Introduction
Calculation of plasma equilibria is a crucial component in the design and operation of the magnetic confinement devices employed in nuclear fusion research. The equilibrium poloidal flux in axisymmetric toroidal magnetic confinement devices is determined by solving the Grad-Shafranov equation [1][2][3], which is a nonlinear, elliptic, partial differential equation. The solution of the Grad-Shafranov equation is usually effected via an iterative numerical scheme.
In 1968, Solov'ev [4] obtained a family of exact analytic solutions of the Grad-Shafranov equation. Solov'ev's solutions are useful for benchmarking plasma equilibrium codes [5,6], as well as for magnetohydrodynamical stability analysis of tokamak plasmas [7]. Solov'ev's solutions have been employed to construct model up-down-symmetric tokamak equilibria [8], as well as equilibria with magnetic divertors.
In this work, we construct a tokamak plasma equilibrium generated by a combination of currents flowing within the plasma and in distant external magnetic field coils. The plasma toroidal current density takes the simple form j ϕ (r, z) = −a r − b R 2 /r inside the plasma (where r and z are cylindrical coordinates, and a, b, and R are constants), and is zero in the surrounding vacuum. The Grad-Shafranov equation possesses a well-known exact analytic solution within the plasma due to Solov'ev. We use a Green's function method to compute the poloidal magnetic flux generated by plasma currents, together with a parameterized homogeneous solution to the Grad-Shafranov equation (that is well-behaved at small r and z), to construct the vacuum solution. The vacuum solution is matched to the analytic solution on the last closed magnetic flux surface (LCFS) to determine the parameters in the homogeneous solution. This procedure is performed for both up-down-symmetric double-null and up-down-asymmetric single-null equilibria. We find that any magnetic X-points on the LCFS are distorted due to the fact that one quadrant is filled by a current-carrying plasma, whereas the other three are filled by a vacuum in which no current flows. In particular, the vacuum quadrant opposite the plasma-filled quadrant expands at the expense of the other three quadrants.

Letter
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Solov'ev's analytic solutions of the Grad-Shafranov equation (as well as the other aforementioned analytic solutions) are derived on the assumption that there is distributed toroidal current filling all space. In reality, the plasma current density is zero in the vacuum region that lies beyond the last closed (magnetic) flux surface (LCFS). Now, the solution of the Grad-Shafranov equation in the vacuum region consists of two components. The first component is the contribution from the toroidal currents flowing within the plasma. The second component is the contribution from currents flowing in external magnetic field coils. In this study, for the sake of simplicity, the external field coils are assumed to be very distant from the plasma. In order to construct the correct solution in the vacuum region, we use a Green's function method to compute the contribution from the plasma currents, and then combine this with a parameterized homogeneous solution of the Grad-Shafranov equation that is well behaved on the toroidal symmetry axis. The parameters that characterize the homogeneous solution are determined by demanding that, when the homogeneous solution is added to the plasma solution, the net solution is constant on the LCFS. In this manner, we obtain a realistic vacuum solution that is consistent with Solov'ev's analytic solution within the plasma.
This paper is organized as follows. The Grad-Shafranov equation is introduced in section 2. Solov'ev's analytic solution, and its problem, is discussed in section 3. Our method for obtaining the vacuum field for up-down-symmetric double-null equilibria is presented in section 4. Our method is extended to deal with up-down-asymmetric single-null equilibria in section 5. The paper is summarized in section 6.

Grad-Shafranov equation
Let r, ϕ, z be conventional cylindrical coordinates. In the following, all lengths are normalized to some convenient scale length, R 0 , all magnetic field strengths to some convenient scale field, B 0 , all poloidal magnetic fluxes to R 2 0 B 0 , all toroidal current densities to B 0 /(µ 0 R 0 ), and all plasma pressures to B 2 0 /µ 0 . (In fact, B 0 is the vacuum toroidal magnetic field strength at r = R 0 .) Consider an axisymmetric toroidal plasma equilibrium whose magnetic field is written Here, ψ(r, z) is the poloidal magnetic flux, and I A (ψ) is an arbitrary flux function that takes the value unity in the vacuum region surrounding the plasma. As is well known, the poloidal flux satisfies the Grad-Shafranov equation [1][2][3], where j ϕ (r, z) = r p + is the toroidal current density. Here p(ψ) is the plasma pressure, which is an arbitrary flux function that takes the value zero in the vacuum region surrounding the plasma. Moreover, denotes a derivative with respect to argument.

Solov'ev's solution and its problem
Following Solov'ev [4], let us assume that where a, b, and R are constants. This implies that In this case, equation (2) possesses the exact analytic solution [4] ψ(r, z) where c 0 is a constant, and As illustrated in figure 1, for a > c 0 > −b, the contours of ψ(r, z) map out a family of up-down-symmetric plasma equilibria with positive triangularity. The magnetic axis is located at r = R, z = 0. Note that ψ = 0 on the axis. The magnetic flux surfaces in the immediate vicinity of the magnetic axis are nested ellipses whose major radii are aligned with the r and z axes. The constant ratio of the major radii of these ellipses is l z /l r = [(a − c 0 )/(b + c 0 )] 1/2 . The magnetic configuration possesses two magnetic X-points, located at ζ = ζ X and z = ±z X , where Both X-points lie on the LCFS, ψ = ψ X , where Generally speaking, we expect the region beyond the LCFS to be a vacuum (or near vacuum) in which the plasma current density is zero. However, according to equation (6), Solov'ev's analytic solution of the Grad-Shafranov equation assumes a toroidal current density distribution that extends to infinity, and is, therefore, not zero in the vacuum region. It is precisely this deficiency of the Solvov'ev solution that we aim to remedy in our paper.

Our method for double-null equilibria
Suppose that the toroidal current density within the plasma is the same as that assumed in the Solov'ev calculation, but that there is zero current density in the vacuum region outside the LCFS. Now, the poloidal magnetic flux in the vacuum region consists of two components. The first component, ψ p (r, z), is the flux generated by the known toroidal current distribution within the LCFS. We will use a Green's function method to compute this component. The second component, ψ h (r, z), is the contribution from currents flowing in remote external magnetic field coils. This component is a homogeneous solution to the Grad-Shafranov equation that is well behaved at small r and z.
In order to calculate ψ p (r, z), it is helpful to define the flux coordinates ψ and θ, where At given values of ψ and θ, we can determine ζ by combining equations (7) and (12) to give the cubic equation (13) For 0 < ψ < ψ X (i.e. within the LCFS), this equation possesses three real roots. For cos θ > 0 the largest root is appropriate, whereas the intermediate root is appropriate when cos θ < 0. Once ζ has been determined, z follows from equation (12). Now, where which implies that which yields (21) The poloidal magnetic flux, ψ p (r, z), generated by the toroidal currents flowing within the plasma is written where the integral is over the plasma poloidal cross-section, and j ϕ (r, z) is given by equation (6). The Green's function takes the form [12] and It follows that Here, and η are small regularization parameters. The poloidal magnetic flux, ψ h (r, z), generated by currents flowing in distant magnetic field coils is written [13] where the c N are arbitrary coefficients, and The coefficients, c N , that characterize the homogenous solution, (32), are determined by demanding that In order to realize our method, the ψ − θ integral in expression (28) is implemented using a two-dimensional (2D) trapezium rule on a 1000 × 1000 uniform rectangular grid. The regularization parameters and η take the values 10 −12 and 10 −6 , respectively. We find that λ ≡ |ψ(R, 0)|/ψ X is a good measure of the accuracy of both the integration scheme and the matching procedure. Of course, λ, which is the relative magnitude of the normalized poloidal magnetic flux on the magnetic axis, should be equal to zero. Obviously, if N h is too small then the matching procedure is inaccurate because the truncated series (32) cannot adequately represent the true homogeneous solution of the Grad-Shafranov equation. On the other hand, if N h is too large then the matching procedure becomes inaccurate due to rounding errors. We find that λ is minimized, taking the value 5 × 10 −6 , when N h = 10.
Let ψ s (r, z) denote the Solov'ev solution, (7), and let ψ(r, z) represent our solution. We find that the value of |ψ − ψ s |/ψ X , averaged over the plasma volume, is 6.6 × 10 −4 , whereas the maximum value is 5.6 × 10 −3 . Figure 2 shows the difference between our solution and the Solov'ev solution on the LCFS. Of course, these two solutions should be identical on this fluxsurface. It can be seen that the relative deviation between the two solutions peaks at the X-points, taking the value 6 × 10 −3 . Figures 3 and 4 show ψ(r, z)/ψ X profiles calculated by our method along horizontal and vertical lines in the r, z plane that pass through the magnetic axis. These profiles are compared to the equivalent Solov'ev profiles. It can be seen that our method accurately reconstructs the Solov'ev analytic solution within the plasma, but deviates strongly from this solution in the surrounding vacuum region. Moreover, both ψ and its spatial derivatives are clearly continuous across the LCFS.
To reiterate, inside the plasma, our solution for ψ(r, z) coincides with a known analytic solution of the Grad-Shafranov equation. Moreover, outside the plasma, our solution is a homogeneous solution of the Grad-Shafranov equation that is generated by a combination of the known plasma currents and currents at infinity. Finally, our solution and its spatial derivatives are continuous across the LCFS. These facts conclusively demonstrate the validity of our method. Figure 5 shows contours of ψ(r, z) calculated by our method for the case in hand. It can be seen, by comparison with figure 1, that the contours coincide with those of the Solov'ev solution inside the plasma, but are very different in the vacuum region surrounding the plasma. Note, in particular, the distorted magnetic X-points on the LCFS. There is a simple explanation for this phenomenon. In the Solov'ev solution, it is easily demonstrated that the two ψ = ψ X contours that cross at either the upper or the lower X-point subtend an acute angle with one another. For the case in hand, Θ p = 71.4 • . In our solution, we expect the ψ = ψ X contours bordering the plasma-filled quadrant of the X-point to subtend the same angle. On the other hand, in accordance with a well-known general principle, we expect the ψ = ψ X contours bordering the three vacuum-filled quadrants to cross at right angles. There is an obvious inconsistency here, which is due to the fact that there is a significant toroidal current density present in the plasma-filled quadrant, all the way to the magnetic separatrix, whereas there is zero current density in the vacuum-filled quadrants. The net result is that the vacuum-filled quadrant directly opposite the plasma-filled quadrant becomes distended.

Our method for single-null equilibria
Conventional tokamak equilibria are generally up-downasymmetric, with a single magnetic null on the LCFS that lies below the magnetic axis. We can easily extend our analysis to take such equilibria into account by writing Note that the previous expression is an exact analytic solution of the Grad-Shafranov equation corresponding to the toroidal current distribution (6). The magnetic axis still lies at r = R, z = 0, ψ = 0. The magnetic X-points lie at ζ = ζ 1 , z = z 1 , ψ = ψ 1 and ζ = ζ 2 , z = z 2 , ψ = ψ 2 . Here, z 1 and z 2 are the positive and negative roots, respectively, of the quadratic equation (38) Once z 1 and z 2 have been determined, ζ 1 , ζ 2 , ψ 1 , and ψ 2 are obtained from If c 1 < 0 then 0 < ψ 2 < ψ 1 . Hence, we can take ψ = ψ 2 as the LCFS. In this case, the LCFS only intersects a single magnetic X-point, which lies in the lower half plane (z < 0). From now on, we shall write ψ 2 = ψ X , for the sake of consistency with our previous analysis. Slightly generalizing our previous analysis, the poloidal magnetic flux generated by the toroidal plasma currents flowing within the plasma is The poloidal magnetic flux generated by currents flowing in distant magnetic field coils takes the form [13] where the c N are arbitrary coefficients. When N is even, and On the other hand, when N is odd, The coefficients, c N , that characterize the homogenous solution, (45), are determined by demanding that on the control surface ψ = (1 − η) ψ X . Given that there are N h independent c N values, the matching must be performed at N h different values of θ. These values are evenly distributed in the range 0 to 2π, with a small common offset that is chosen so as to ensure that one of the matching points corresponds to the magnetic X-point on the LCFS. Our second example is an up-down-asymmetric singlenull plasma equilibrium calculated with the following parameters: R = 1.0, a = 1.2, b = −1.0, c 0 = 1.1, and c 1 = −0.005. The corresponding Solov'ev solution is shown in figure 6.
As before, the ψ − θ integral in expression (41) is implemented using a 2D trapezium rule on a 1000 × 1000 uniform rectangular grid. The regularization parameters and η again take the values 10 −12 and 10 −6 , respectively. We find that the previously defined error parameter, λ, attains its minimum value, 9 × 10 −4 , when N h = 18.
We find that the value of |ψ − ψ s |/ψ X , averaged over the plasma volume, is 3.4 × 10 −3 , whereas the maximum value is 7.7 × 10 −3 . Figure 7 shows the difference between our solution and the Solov'ev solution on the LCFS. Of course, these two solutions should be identical on this flux-surface. It can be seen that the relative deviation between the two solutions peaks at the X-point, taking the value 3 × 10 −3 . Figure 8 shows contours of ψ(r, z) calculated by our method for the case in hand. It can be seen, by comparison with figure 6, that the contours coincide with those of the Solov'ev solution inside the plasma, but are very different in the vacuum region surrounding the plasma. As before, the magnetic X-point on the LCFS is distorted in such a manner that the vacuum-filled quadrant directly opposite the plasmafilled quadrant becomes distended. The explanation for this phenomenon is the same as that given in the previous section.

Summary
In this paper, we calculate a tokamak equilibrium generated by a combination of currents flowing within the plasma and in distant external magnetic field coils. The poloidal magnetic flux surfaces inside the plasma are identical to those given by the well-known exact analytic solution of the Grad-Shafranov equation due to Solov'ev. However, the flux surfaces outside the plasma correspond to a properly matched vacuum solution. We use a Green's function method to compute the poloidal magnetic flux generated by the known plasma cur rents, together with a parameterized homogeneous solution to the Grad-Shafranov equation (that is well-behaved at small r and z), to construct the vacuum solution. The vacuum solution is matched to the analytic solution on the last closed magnetic flux surface (LCFS) to determine the parameters in the homogeneous solution. This procedure is performed for both up-down-symmetric double-null and up-down-asymmetric single-null equilibria. We find that any magnetic X-points on the LCFS are distorted due to the fact that one quadrant is filled by a current-carrying plasma, whereas the other three are filled by a vacuum in which no current flows. In par ticular, the vacuum quadrant directly opposite the plasma-filled quadrant expands at the expense of the other three quadrants. We conjecture that this phenomenon could be a feature of any magn etically diverted tokamak equilibrium in which there is a large toroidal current density at edge of the plasma (as is indeed generally the case in H-mode equilibria). Finally, the method described in this paper could easily be modified to deal with the extended Solov'ev-type equilibria described in [9][10][11].  necessarily subtend 90 degrees if the plasma current on the LCFS is non-zero, which is inconsistent with the standard scenario in which the field-lines that cross at a magnetic X-point are assumed to do so at right-angles. The authors also thank F.L. Waelbroeck (IFS, University of Texas at Austin) for helpful discussions. This work was funded by the US Department of Energy under Contract DE-FG02-04ER-54742. T.X. was funded by the China Scholarship Council and the ITER Special Foundation of China under Contract 2017YFE0301202.