%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
This paper was motivated by Shafranov who in 1998 attended the American Physical Society meeting. He was my former supervisor and an expert in equilibrium and stability theory of stellarators (Shafranov & Yurchenko Reference Shafranov and Yurchenko1969; Pustovitov & Shafranov Reference Pustovitov, Shafranov and Kadomtsev1990). I was especially pleased to present to him the new concept of reference magnetic coordinates (RMC) invented just two months earlier. He was highly impressed: the definition of RMC and the fast algorithm for their generation could solve the problem which he had known since the 1950s, i.e. how to generate the nested coordinates with the best alignment with a 3-D ergodic magnetic field. The RMC, published recently (Li & Zakharov Reference Li and Zakharov2013; Zakharov and Li Reference Zakharov and Li2015), exclude the widely used, but highly inefficient, field-line tracing technique for generating the nested coordinates. Instead, RMC rely on Fourier transformation, which can further enhance the numerical performance by using graphic processors (GPU).
After our discussion regarding RMC, he pointed out another very old problem in stellarator equilibria related to their nestedness and the Hamada principle. Countless theoretical papers on 3-D equilibria rely on the assumption of nested surfaces but the nestedness in the 3-D case is exceptional. Elimination of islands (especially with codes which do not use RMC) is a tedious problem. Shafranov suggested to me to investigate how to generate the nested equilibria in a more straightforward way. The solution came in 2004 and the present paper describes it for the first time.
The difficulties in calculation of 3-D equilibria are related to the absence of a special function, which could determine the shape of the magnetic surfaces as in 2-D axisymmetric configurations. Due to the formation of magnetic islands on resonant magnetic surfaces, it is impossible to introduce a simple nested toroidal coordinate system perfectly aligned with the magnetic field.
Throughout this paper, we consider the following Fourier representation of the cylindrical coordinates
$r,{\it\varphi},z$
$$\begin{eqnarray}\displaystyle \displaystyle r=\mathop{\sum }_{m,n}r_{mn}(a)\text{e}^{-\text{i}(m{\it\theta}-n{\it\zeta})},\quad \displaystyle z=\mathop{\sum }_{m,n}z_{mn}(a)\text{e}^{-\text{i}(m{\it\theta}-n{\it\varphi})},\quad \displaystyle {\it\varphi}=\mathop{\sum }_{m,n}{\it\varphi}_{mn}(a)\text{e}^{-\text{i}(m{\it\theta}-n{\it\zeta})} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
in terms of toroidal coordinates
$a,{\it\theta},{\it\zeta}$
, where
$a=\text{const.}$
represents a toroidal coordinate surface,
${\it\theta}$
is a poloidal angle and
${\it\zeta}$
is a toroidal angle. The Fourier representation is used also for all oscillating functions in order to convert the differentiation over angles into multiplications.
In practice, two types of numerical codes are in use. One type (e.g. the VMEC code (Hirshman and Lee Reference Hirshman and Lee1986; Hirshman, van Rij & Merkel Reference Hirshman, van Rij and Merkel1986) as the best example) works with the simple nested toroidal coordinates
$a,{\it\theta},{\it\zeta}$
.
The radial Fourier harmonics
$r_{mn}(a),z_{mn}(a),{\it\varphi}_{mn}(a)$
in (1.1) are adjusted to minimize the Kruskal–Kulsrud (K–K) variational functional
$W$
(Kruskal & Kulsrud Reference Kruskal and Kulsrud1958)
$$\begin{eqnarray}\displaystyle W=\int \left(\frac{|\boldsymbol{B}|^{2}}{2{\it\mu}_{0}}+\frac{p}{{\it\gamma}-1}\right)\text{d}^{3}r, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{B}$
is the magnetic field,
$p$
is the plasma pressure and
${\it\mu}_{0}$
is the magnetic permeability of vacuum.
It is widely believed (see, e.g. the review Helander (Reference Helander2014)) that the minimization of the K–K variational principle leads to equilibrium configurations. In fact, this is not the case.
The nested magnetic fields are described by a vector potential
$\boldsymbol{A}$
in the form
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{A}=-\bar{{\it\Phi}}^{\prime }{\it\eta}\boldsymbol{{\rm\nabla}}a+\bar{{\it\Phi}}\boldsymbol{{\rm\nabla}}{\it\theta}+\bar{{\it\Psi}}\boldsymbol{{\rm\nabla}}{\it\zeta},\quad {\it\eta}={\it\eta}(a,{\it\theta},{\it\zeta}),\bar{{\it\Phi}}=\bar{{\it\Phi}}(a),\bar{{\it\Psi}}=\bar{{\it\Psi}}(a).\quad & \displaystyle\end{eqnarray}$$
Here,
$2{\rm\pi}\bar{{\it\Phi}}(a),2{\rm\pi}\bar{{\it\Psi}}(a)$
are the toroidal and poloidal magnetic flux through the contours
$a=\text{const.}$
,
${\it\varphi}=\text{const.}$
and
$a=\text{const.}$
,
${\it\theta}=\text{const}$
. The nestedness is manifested by the absence of the angle dependence of both fluxes. Indeed, the magnetic field
$\boldsymbol{B}=(\boldsymbol{{\rm\nabla}}\times \boldsymbol{A})$
has the form
$$\begin{eqnarray}\displaystyle \boldsymbol{B}=(\bar{{\it\Phi}}^{\prime }+\bar{{\it\Phi}}^{\prime }{\it\eta}_{{\it\theta}}^{\prime })(\boldsymbol{{\rm\nabla}}a\times \boldsymbol{{\rm\nabla}}{\it\theta})-(\bar{{\it\Psi}}^{\prime }+\bar{{\it\Phi}}^{\prime }{\it\eta}_{{\it\varphi}}^{\prime })(\boldsymbol{{\rm\nabla}}{\it\zeta}\times \boldsymbol{{\rm\nabla}}a) & & \displaystyle\end{eqnarray}$$
and does not contain a component normal to the surface
$a=\text{const.}$
More details on this representation for both nested and ergodic magnetic fields can be found in Zakharov and Li (Reference Zakharov and Li2015).
For the sake of clarity, for consideration of 3-D configurations, we adopt here the straight field-line coordinates
$a,\bar{{\it\theta}},\bar{{\it\zeta}}$
with angles
$\bar{{\it\theta}},\bar{{\it\zeta}}$
adjusted to eliminate the function
${\it\eta}$
in
$\boldsymbol{A}$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{A}=\bar{{\it\Phi}}\boldsymbol{{\rm\nabla}}\bar{{\it\theta}}+\bar{{\it\Psi}}\boldsymbol{{\rm\nabla}}\bar{{\it\zeta}},\quad {\it\eta}=0. & \displaystyle\end{eqnarray}$$
The minimization of the K–K functional, performed with conservation of magnetic fluxes, leads to a global equilibrium but also the generation of delta-functional singular currents
$$\begin{eqnarray}\displaystyle i_{m^{\ast }n^{\ast }}={\it\delta}(a-a_{m^{\ast }n^{\ast }})\text{e}^{\text{i}(m^{\ast }\bar{{\it\theta}}-n^{\ast }\bar{{\it\zeta}})} & & \displaystyle\end{eqnarray}$$
at the resonant magnetic surfaces, where
$$\begin{eqnarray}\displaystyle m^{\ast }\bar{{\it\Psi}}^{\prime }+n^{\ast }\bar{{\it\Phi}}^{\prime }=0. & & \displaystyle\end{eqnarray}$$
The existence of such surface currents is a commonly known in ideal MHD. The simple illustration would be a cylindrical,
${\it\rho},{\it\omega}$
, plasma with a current profile stable to a tearing mode, e.g.
$m^{\ast }/n^{\ast }=2/1$
. If the configuration is deformed by a helical perturbation of the plasma boundary (or by a helical external field perturbation), then there are two possible internal solutions satisfying the boundary conditions. One corresponds to a tearing mode with a helical vector potential
${\it\psi}$
continuous across the resonant surface
${\it\rho}={\it\rho}_{m^{\ast }n^{\ast }}$
. This solution is associated with an island
$m^{\ast }/n^{\ast }=2/1$
. Another, the ideal MHD solution, has the plasma displacement
${\it\xi}_{m^{\ast }n^{\ast }}$
continuous across
${\it\rho}_{m^{\ast }n^{\ast }}$
and a zero value of the helical vector potential
${\it\psi}_{m^{\ast }n^{\ast }}=0$
. The jump in radial derivative of
${\it\psi}_{m^{\ast }n^{\ast }}^{\prime }({\it\rho}-{\it\epsilon})-{\it\psi}_{m^{\ast }n^{\ast }}^{\prime }({\it\rho}+{\it\epsilon})={\it\mu}_{0}i_{m^{\ast }n^{\ast }}$
indicates the existence a surface (sheet) current
$i_{m^{\ast }n^{\ast }}$
.
In 3-D configurations, after minimization of the K–K functional, such surface currents
$i_{m^{\ast }n^{\ast }}$
are in equilibrium along the normal direction to the resonant surfaces: the total magnetic and kinetic pressures at both sides of the resonant surfaces are equal in amplitude but opposite in direction. At the same time, the plasma model of an ideal electric conductivity used in the K–K principle does not resolve the internal structure of these sheet currents. Due to resistive (and collisionless) reconnection effects, they have a finite width. Due to the oscillatory nature of the sheet currents along the resonant surface, the amplitude of the magnetic pressure is also oscillatory. As a result, the plasma is not in equilibrium along the resonant surface. In reality it is squeezed from the regions with a higher magnetic pressure to the regions with a lower pressure by a self-consistent reconnection mechanism, leading to formation of the tearing type of magnetic islands.
This consideration shows that the K–K (Kruskal & Kulsrud Reference Kruskal and Kulsrud1958) variational principle cannot lead to rigorous equilibrium configurations (although it provides practical algorithms for equilibrium calculations if the resonant surface effects can be neglected). Although the generation of the sheet currents by the use of the K–K energy principle is well-known (see, e.g. the recent (Loizu et al. Reference Loizu, Hudson, Bhattacharjee and Helander2015)), the absence of an equilibrium is not well recognized (Helander Reference Helander2014). Accordingly, the problem of the elimination of both the surface currents and islands has not been properly addressed. It is addressed in the present paper by a direct implementation of the Hamada principle.
Other types of numerical code (e.g. PIES (Reiman & Greenside Reference Reiman and Greenside1986)) use laboratory coordinates which are not adjusted to the structure of the magnetic field. These codes can solve the field equations regardless of the splitting of the resonant surfaces by islands and magnetic stochasticity. Nevertheless, during the iterations, the specification of the current distribution, corresponding to equilibrium constraints, requires some at least approximate geometry of magnetic surfaces in the situation where rigorous ones do not exist. This represents a conceptual problem for this kind of code. Then, the freedom in choosing boundary conditions (plasma shape or external field) can be used to generate an equilibrium with negligible islands.
Unfortunately, the only technique available for reconstruction of magnetic surfaces in these codes is field-line tracing of the magnetic field. Computationally, this is very time consuming and line tracing becomes more and more questionable for small size islands, i.e. in the vicinity of the desirable nested surface configurations.
A rigorous and fast (based on Newton iteration scheme) algorithm for calculation of 3-D equilibrium configurations with islands is described in Zakharov and Li (Reference Zakharov and Li2015). The so-called reference magnetic coordinates (RMC) were introduced instead of the non-existent nested flux coordinates. Also, the time consuming field-line tracing is replaced by a fast algorithm for iteratively advancing the RMC. The algorithm is based on a variational principle, which is not restricted by flux conservation.
Considering the problem of calculation of a 3-D equilibrium with islands solved conceptually, the present paper is focused on generating the equilibrium configurations with nested magnetic surfaces (without islands). The necessary condition for their existence was formulated by Hamada (Reference Hamada1962). It requires that either the Jacobian
$J(a,\bar{{\it\theta}},\bar{{\it\zeta}})$
of the straight field line
$a,\bar{{\it\theta}},\bar{{\it\zeta}}$
coordinates
$$\begin{eqnarray}\displaystyle & \displaystyle J=r\frac{D(r,{\it\varphi},z)}{D(a,\bar{{\it\theta}},\bar{{\it\zeta}})} & \displaystyle\end{eqnarray}$$
does not contain resonant harmonics at the resonant surfaces or that the plasma pressure gradient should be zero there.
The fundamental step, made here for an implementation of the Hamada principle, consists in imposing the requirement that the calculations preserve a prescribed Jacobian with no resonant harmonics. It is obvious that, e.g. the Hamada-like Jacobian
$J=J(a)$
automatically implements the Hamada principle by having no Fourier spectrum at all.
Section 2 demonstrates this approach and its specifics using 2-D equilibrium configurations. Section 3 introduces the linearized equilibrium equations (LEE) as a fast algorithm for the calculation of nested equilibria. Section 4 extends the approach to 3-D equilibrium configurations, and Summary 5 outlines the importance of the topic, potential implications for fusion research and necessary advances for practical implementation of the LEE in numerical codes.
2. Axisymmetric equilibrium equations with a prescribed Jacobian
The metric tensor of curvilinear coordinates
$a,\bar{{\it\theta}},\bar{{\it\zeta}}$
is determined from (1.1) as
$$\begin{eqnarray}\displaystyle \text{d}l^{2}=g_{aa}\,\text{d}a^{2}+2g_{a\bar{{\it\theta}}}\,\text{d}a\,\text{d}\bar{{\it\theta}}+2g_{a\bar{{\it\zeta}}}\,\text{d}a\,\text{d}\bar{{\it\zeta}}+g_{\bar{{\it\theta}}\bar{{\it\theta}}}\,\text{d}\bar{{\it\theta}}^{2}+2g_{\bar{{\it\theta}}\bar{{\it\theta}}}\,\text{d}\bar{{\it\theta}}\,\text{d}\bar{{\it\zeta}}+g_{\bar{{\it\zeta}}\bar{{\it\zeta}}}\,\text{d}\bar{{\it\zeta}}^{2}, & & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle g_{aa}=r_{a}^{\prime 2}+z_{a}^{\prime 2}+r^{2}{\it\varphi}_{a}^{\prime },\quad g_{a{\it\theta}}=r_{a}^{\prime }r_{{\it\theta}}^{\prime }+z_{a}^{\prime }z_{{\it\theta}}^{\prime }+r^{2}{\it\varphi}_{a}^{\prime }{\it\varphi}_{{\it\theta}}^{\prime },\quad g_{a{\it\zeta}}=r_{a}^{\prime }r_{{\it\zeta}}^{\prime }+z_{a}^{\prime }z_{{\it\zeta}}^{\prime }+r^{2}{\it\varphi}_{a}^{\prime }{\it\varphi}_{{\it\zeta}}^{\prime }, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle g_{{\it\theta}{\it\theta}}=r_{{\it\theta}}^{\prime 2}+z_{{\it\theta}}^{\prime 2}+r^{2}{\it\varphi}_{{\it\theta}}^{\prime 2},\quad g_{{\it\theta}{\it\zeta}}=r_{{\it\theta}}^{\prime }r_{{\it\zeta}}^{\prime }+z_{{\it\theta}}^{\prime }z_{{\it\zeta}}^{\prime }+r^{2}{\it\varphi}_{{\it\theta}}^{\prime }{\it\varphi}_{{\it\zeta}}^{\prime },\quad g_{{\it\zeta}{\it\zeta}}=r_{{\it\zeta}}^{\prime 2}+z_{{\it\zeta}}^{\prime 2}+r^{2}{\it\varphi}_{{\it\zeta}}^{\prime 2} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
and the Jacobian given by
$$\begin{eqnarray}\displaystyle J=rD,\quad D=z_{a}^{\prime }(r_{{\it\theta}}^{\prime }{\it\varphi}_{{\it\zeta}}^{\prime }-r_{{\it\zeta}}^{\prime }{\it\varphi}_{{\it\theta}}^{\prime })+r_{a}^{\prime }(z_{{\it\zeta}}^{\prime }{\it\varphi}_{{\it\theta}}^{\prime }-z_{{\it\theta}}^{\prime }{\it\varphi}_{{\it\zeta}}^{\prime })+{\it\varphi}_{a}^{\prime }(r_{{\it\zeta}}^{\prime }z_{{\it\theta}}^{\prime }-r_{{\it\theta}}^{\prime }z_{{\it\zeta}}^{\prime }). & & \displaystyle\end{eqnarray}$$
For the 2-D axisymmetric configurations considered in this section, the assumption of straight field-line coordinates is not essential and we consider
$\bar{{\it\zeta}}={\it\varphi}$
and
${\it\theta}$
an arbitrary poloidal angle
$$\begin{eqnarray}\displaystyle D=z_{a}^{\prime }r_{{\it\theta}}^{\prime }-r_{a}^{\prime }z_{{\it\theta}}^{\prime }. & & \displaystyle\end{eqnarray}$$
In curvilinear flux coordinates
$\bar{{\it\Psi}}=\bar{{\it\Psi}}(a)$
the 2-D Grad–Shafranov (GSh) equation has the form
$$\begin{eqnarray}\displaystyle \frac{1}{J}\left(\frac{g_{{\it\theta}{\it\theta}}}{J}\bar{{\it\Psi}}^{\prime }\right)_{a}^{\prime }-\frac{1}{J}\left(\frac{g_{a{\it\theta}}}{J}\bar{{\it\Psi}}^{\prime }\right)_{{\it\theta}}^{\prime }+P+\frac{1}{r^{2}}T=0,\quad P\equiv \frac{\text{d}\bar{p}}{\text{d}\bar{{\it\Psi}}},T\equiv \bar{F}\frac{\text{d}\bar{F}}{\text{d}\bar{{\it\Psi}}} & & \displaystyle\end{eqnarray}$$
with
$\bar{p}\equiv {\it\mu}_{0}p(a)$
representing the plasma pressure and
$\bar{F}=\bar{F}(a)\equiv rB_{{\it\varphi}}$
. This equation should be considered to be a nonlinear equation for determination of the flux coordinates
$r(a,{\it\theta})$
,
$z(a,{\it\theta})$
.
When the coordinates
$a,{\it\theta}$
are given, (2.6) is a second-order differential equation for the Fourier harmonics
$r_{mn}(a)$
and
$z_{mn}(a)$
with one
$a$
-derivative hidden in the Jacobian
$J$
and the second one being explicit. Thus, for fixed boundary problems, the plasma boundary is presented in a parametric form
$$\begin{eqnarray}\displaystyle r^{b}=R({\it\omega}),\quad z^{b}=Z({\it\omega}), & & \displaystyle\end{eqnarray}$$
with
${\it\theta}\equiv {\it\omega}$
.
On the other hand, if
$J$
is a given function of coordinates, the first radial derivative disappears and (2.6) becomes the first-order differential equation for the Fourier harmonics of
$r,z$
. Complemented by a definition of the Jacobian (2.4) and a boundary condition, it is equivalent to the second-order GSh equation. In particular, the resulting functions
$r_{mn}(a)$
and
$z_{mn}(a)$
determine the angle
${\it\theta}$
. The regularity condition at the magnetic axis leaves only one set of harmonics
$r_{mn}(a),z_{mn}(a)$
as independent. The remaining freedom is used to satisfy the boundary conditions at the plasma edge.
In the case with a prescribed Jacobian, the boundary conditions represent a nonlinear problem for conversion of
${\it\omega}$
to
${\it\theta}$
$$\begin{eqnarray}\displaystyle r^{b}({\it\theta})=R({\it\omega}({\it\theta})),\quad z^{b}({\it\theta})=Z({\it\omega}({\it\theta})). & & \displaystyle\end{eqnarray}$$
In the next section, the linearization technique will be proposed for solving this problem.
3. Linearized equilibrium equations (LEE)
The most efficient way to solve a nonlinear problem is to linearize it. By introducing the notation
$$\begin{eqnarray}\displaystyle Q\equiv \frac{\bar{{\it\Psi}}^{\prime }}{J}, & & \displaystyle\end{eqnarray}$$
the system of equilibrium equations with a prescribed Jacobian can be written as
$$\begin{eqnarray}\displaystyle JG\equiv (g_{{\it\theta}{\it\theta}}Q)_{a}^{\prime }-(g_{a{\it\theta}}Q)_{{\it\theta}}^{\prime }+JP+\frac{J}{r^{2}}T=0,\quad J=rD,D\equiv z_{a}^{\prime }r_{{\it\theta}}^{\prime }-r_{a}^{\prime }z_{{\it\theta}}^{\prime }. & & \displaystyle\end{eqnarray}$$
The remarkable property of this equation is that it can be easily linearized with respect to the perturbations of the unknowns
${\it\delta}r(a,{\it\theta}),{\it\delta}z(a,{\it\theta})$
and
${\it\delta}Q$
giving a set of linearized equilibrium equations
$$\begin{eqnarray}\displaystyle & \displaystyle (g_{{\it\theta}{\it\theta}}{\it\delta}Q+Q{\it\delta}g_{{\it\theta}{\it\theta}})_{a}^{\prime }-(g_{a{\it\theta}}{\it\delta}Q+Q{\it\delta}g_{a{\it\theta}})_{{\it\theta}}^{\prime }-2\frac{J{\it\delta}r}{r^{3}}T=-JG,\quad r{\it\delta}D+D{\it\delta}r=J-rD, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where the perturbations of the metric tensor and
$Q$
are simple
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}g_{{\it\theta}{\it\theta}}=2r_{{\it\theta}}^{\prime }{\it\delta}r_{{\it\theta}}^{\prime }+2z_{{\it\theta}}^{\prime }{\it\delta}z_{{\it\theta}}^{\prime },\quad {\it\delta}g_{a{\it\theta}}=r_{a}^{\prime }{\it\delta}r_{{\it\theta}}^{\prime }+r_{{\it\theta}}^{\prime }{\it\delta}r_{a}^{\prime }+z_{a}^{\prime }{\it\delta}z_{{\it\theta}}^{\prime }+z_{{\it\theta}}^{\prime }{\it\delta}z_{a}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}Q=\frac{{\it\delta}\bar{{\it\Psi}}}{J},\quad {\it\delta}D=r_{a}^{\prime }{\it\delta}z_{{\it\theta}}^{\prime }+z_{{\it\theta}}^{\prime }{\it\delta}r_{a}^{\prime }-z_{a}^{\prime }{\it\delta}r_{{\it\theta}}^{\prime }-r_{{\it\theta}}^{\prime }{\it\delta}z_{a}^{\prime }. & \displaystyle\end{eqnarray}$$
The right-hand sides in the linearized equations (3.3) represent the discrepancy. In our linearization, we consider the source functions
$P(a),T(a)$
as fixed. In general, they can also be perturbed in accordance with other constraints, e.g. the prescribed
$p(a)$
and
$q(a)$
, as is explained in Zakharov & Pletzer (Reference Zakharov and Pletzer1999). The radial derivatives of the unknowns can be extracted to the left-hand side as
$$\begin{eqnarray}\displaystyle & & \displaystyle g_{{\it\theta}{\it\theta}}{\it\delta}Q_{a}^{\prime }+Qr_{{\it\theta}}^{\prime }{\it\delta}r_{{\it\theta},a}^{\prime \prime }+Qz_{{\it\theta}}^{\prime }{\it\delta}z_{{\it\theta},a}^{\prime \prime }-(Qr_{{\it\theta}}^{\prime })_{{\it\theta}}^{\prime }{\it\delta}r_{a}^{\prime }-(Qz_{{\it\theta}}^{\prime })_{{\it\theta}}^{\prime }{\it\delta}z_{a}^{\prime }\nonumber\\ \displaystyle & & \displaystyle \quad =-Q_{a}^{\prime }{\it\delta}g_{{\it\theta}{\it\theta}}-2Q(r_{{\it\theta},a}^{\prime \prime }{\it\delta}r_{{\it\theta}}^{\prime }+z_{{\it\theta},a}^{\prime \prime }{\it\delta}z_{{\it\theta}}^{\prime })-g_{{\it\theta}{\it\theta},a}^{\prime }{\it\delta}Q+(g_{a{\it\theta}}{\it\delta}Q)_{{\it\theta}}^{\prime }\nonumber\\ \displaystyle & & \displaystyle \qquad +\,(Qr_{a}^{\prime }{\it\delta}r_{{\it\theta}}^{\prime }+Qz_{a}^{\prime }{\it\delta}z_{{\it\theta}}^{\prime })_{{\it\theta}}^{\prime }+2\frac{J{\it\delta}r}{r^{3}}T-JG,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle r(z_{{\it\theta}}^{\prime }{\it\delta}r_{a}^{\prime }-r_{{\it\theta}}^{\prime }{\it\delta}z_{a}^{\prime })=-(r_{a}^{\prime }{\it\delta}z_{{\it\theta}}^{\prime }-z_{a}^{\prime }{\it\delta}r_{{\it\theta}}^{\prime })-D{\it\delta}r+J-rD. & \displaystyle\end{eqnarray}$$
Represented in terms of the Fourier harmonics of unknowns
${\it\delta}r,{\it\delta}z$
, these equations give a standard form of a system of ordinary differential equations (ODE), resolved with respect to its first derivatives.
We assume that a fundamental solution of this system of equations is known (e.g. the shooting technique from the magnetic axis used in the ESC code (Zakharov & Pletzer Reference Zakharov and Pletzer1999) and the DCON code by A. Glasser can generate it)
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\it\delta}r=\mathop{\sum }_{k=0}^{K}C_{k}{\it\delta}r^{k}(a,{\it\theta})=\mathop{\sum }_{k=0}^{K}C_{k}\sum {\it\delta}r_{m}^{k}\text{e}^{-\text{i}m{\it\theta}},\\ \displaystyle {\it\delta}z=\mathop{\sum }_{k=0}^{K}C_{k}{\it\delta}z^{k}(a,{\it\theta})=\mathop{\sum }_{k=0}^{K}A_{k}\sum {\it\delta}z_{m}^{k}\text{e}^{-\text{i}m{\it\theta}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$C_{k}$
are free parameters and
${\it\delta}r^{k}(a,{\it\theta}),{\it\delta}z^{k}(a,{\it\theta})$
are the independent solutions to LEE equations.
The linearized boundary conditions
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \mathop{\sum }_{k=0}^{K}C_{k}{\it\delta}r^{k}(a,{\it\theta})=R({\it\omega})-r^{b}({\it\theta})-r_{{\it\theta}}^{\prime b}{\it\delta}{\it\theta},\\ \displaystyle \mathop{\sum }_{k=0}^{K}C_{k}{\it\delta}z^{k}(a,{\it\theta})=R({\it\omega})-z^{b}({\it\theta})-z_{{\it\theta}}^{\prime b}{\it\delta}{\it\theta},\quad {\it\theta}={\it\theta}({\it\omega})+{\it\delta}{\it\theta}\end{array}\right\} & & \displaystyle\end{eqnarray}$$
determine the coefficients
$C_{k}$
and the perturbation of the angle
${\it\delta}{\it\theta}$
if the initial guess
$r^{b}({\it\theta}),z^{b}({\it\theta}),{\it\theta}({\it\omega})$
is specified.
4. 3-D equilibrium with nested magnetic surfaces
Written for the case of prescribed Jacobian
$J$
, the 3-D equilibrium equations in the following form have in flux coordinates (see, e.g. Zakharov and Li (Reference Zakharov and Li2015) and earlier (Soloviev & Shafranov Reference Soloviev, Shafranov and Leontovich1970; Zakharov & Shafranov Reference Zakharov, Shafranov and Leontovich1986))
$$\begin{eqnarray}\displaystyle & \displaystyle G_{t}\equiv (-g_{{\it\theta}{\it\theta}}Q+g_{{\it\theta}{\it\zeta}}S)_{a}^{\prime }+(g_{a{\it\theta}}Q-g_{a{\it\zeta}}S)_{{\it\theta}}^{\prime }-\bar{J}^{\prime }(a)-{\it\nu}_{{\it\theta}}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{p}\equiv (-g_{{\it\theta}{\it\zeta}}Q+g_{{\it\zeta}{\it\zeta}}S)_{a}^{\prime }+(g_{a{\it\theta}}Q-g_{a{\it\zeta}}S)_{{\it\zeta}}^{\prime }-\bar{F}^{\prime }(a)-{\it\nu}_{{\it\zeta}}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle rD-J=0, & \displaystyle\end{eqnarray}$$
where
$D$
is defined by (2.4). The first two equations represent simply the toroidal and poloidal components of the Ampere law
$(\boldsymbol{{\rm\nabla}}\times (\boldsymbol{{\rm\nabla}}\times \boldsymbol{A}))={\it\mu}_{0}\boldsymbol{j}$
with a special form for the toroidal
$\bar{J}^{\prime }+{\it\nu}_{{\it\theta}}$
and poloidal
$\bar{F}^{\prime }+{\it\nu}_{{\it\zeta}}$
current densities required by the equilibrium equation
$\boldsymbol{{\rm\nabla}}p=(\boldsymbol{j}\times \boldsymbol{B})$
. (The bar symbol distinguishes the toroidal current
$\bar{J}$
from Jacobian
$J$
.) Two of the three 1-D functions
$\bar{p}(a),\bar{J}(a),\bar{F}(a)$
are considered as given.
The functions
$Q,S$
written in arbitrary flux coordinates
$a,{\it\theta},{\it\zeta}$
have the form
$$\begin{eqnarray}\displaystyle Q\equiv \frac{\bar{{\it\Psi}}^{\prime }+\bar{{\it\Phi}}^{\prime }{\it\eta}_{{\it\theta}}^{\prime }}{J},\quad S\equiv \frac{\bar{{\it\Phi}}^{\prime }(1+{\it\eta}_{{\it\theta}}^{\prime })}{J}. & & \displaystyle\end{eqnarray}$$
A practical choice is the use of the azimuth
${\it\varphi}$
as a toroidal angle,
${\it\zeta}\equiv {\it\varphi}$
. In this case, (4.1)–(4.3) determine three unknown functions
$r(a,{\it\theta},{\it\zeta})$
,
$z(a,{\it\theta},{\it\zeta})$
,
${\it\eta}(a,{\it\theta},{\it\zeta})$
and two profiles
$\bar{{\it\Psi}}^{\prime }(a),\bar{{\it\Phi}}^{\prime }(a)$
. In straight field-line coordinates
$a,\bar{{\it\theta}},\bar{{\it\zeta}}$
, the function
${\it\eta}$
is absent, and
${\it\varphi}(a,{\it\theta},{\it\zeta})$
replaces
${\it\eta}$
as an unknown function.
Here we consider the case of straight field-line coordinates only. Equations (4.1)–(4.3) should be complemented by the radial force balance equation which determines the oscillatory part of the current density
${\it\nu}={\it\nu}(a,\bar{{\it\theta}},\bar{{\it\zeta}})$
$$\begin{eqnarray}\displaystyle \bar{p}^{\prime }J=-\bar{{\it\Phi}}^{\prime }(\bar{F}^{\prime }+{\it\nu}_{\bar{{\it\zeta}}})+\bar{{\it\Psi}}^{\prime }(\bar{J}^{\prime }+{\it\nu}_{\bar{{\it\theta}}}). & & \displaystyle\end{eqnarray}$$
The Hamada principle for nested surface equilibrium requires the absence of the resonant harmonics
$m^{\ast },n^{\ast }$
in this expression
$$\begin{eqnarray}\displaystyle J\equiv \mathop{\sum }_{m,n}J_{mn}\text{e}^{\text{i}m\bar{{\it\theta}}-\text{i}n\bar{{\it\zeta}}},\quad J_{m^{\ast }n^{\ast }}=0. & & \displaystyle\end{eqnarray}$$
With the resonant terms absent on the right-hand side, the explicit solution to (4.5) can be written as
$$\begin{eqnarray}\displaystyle \bar{J}^{\prime }\bar{{\it\Psi}}^{\prime }-\bar{F}^{\prime }\bar{{\it\Phi}}^{\prime }=J_{00}\bar{p}^{\prime },\quad \bar{{\it\nu}}_{00}=0,\quad \bar{{\it\nu}}=\mathop{\sum }_{m,n}\frac{-\text{i}\bar{p}^{\prime }J_{mn}}{m\bar{{\it\Psi}}^{\prime }+n\bar{{\it\Phi}}^{\prime }}\text{e}^{\text{i}m\bar{{\it\theta}}-\text{i}n\bar{{\it\zeta}}}. & & \displaystyle\end{eqnarray}$$
In the case of a 1-D Jacobian
$J=J(a)$
, the function
${\it\nu}$
is eliminated.
As in the previous 2-D case, (4.1)–(4.3) in Fourier representation give a set of first-order ODE for Fourier harmonics
$r_{mn},z_{mn},{\it\varphi}_{mn}$
:
$$\begin{eqnarray}\displaystyle G_{t} & = & \displaystyle -Qr_{\bar{{\it\theta}}}^{\prime }r_{a\bar{{\it\theta}}}^{\prime \prime }+Sr_{\bar{{\it\theta}}}^{\prime }r_{a\bar{{\it\zeta}}}^{\prime \prime }-Qz_{\bar{{\it\theta}}}^{\prime }z_{a\bar{{\it\theta}}}^{\prime \prime }+Sz_{\bar{{\it\theta}}}^{\prime }z_{a\bar{{\it\zeta}}}^{\prime \prime }-Qr^{2}{\it\varphi}_{\bar{{\it\theta}}}^{\prime }{\it\varphi}_{a\bar{{\it\theta}}}^{\prime }+Sr^{2}{\it\varphi}_{\bar{{\it\theta}}}^{\prime }{\it\varphi}_{a\bar{{\it\zeta}}}^{\prime \prime }\nonumber\\ \displaystyle & & \displaystyle +\,[2(Sr{\it\varphi}_{\bar{{\it\theta}}}{\it\varphi}_{\bar{{\it\zeta}}}^{\prime }-2Qr{\it\varphi}_{\bar{{\it\theta}}}^{\prime 2})+(r_{\bar{{\it\theta}}}^{\prime }Q-Sr_{\bar{{\it\zeta}}}^{\prime })_{\bar{{\it\theta}}}^{\prime }]r_{a}^{\prime }+(z_{\bar{{\it\theta}}}^{\prime }Q-Sz_{\bar{{\it\zeta}}}^{\prime })_{\bar{{\it\theta}}}^{\prime }z_{a}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,(r^{2}{\it\varphi}_{\bar{{\it\theta}}}^{\prime }Q-Sr^{2}{\it\varphi}_{\bar{{\it\zeta}}}^{\prime })_{\bar{{\it\theta}}}^{\prime }{\it\varphi}_{a}^{\prime }-g_{\bar{{\it\theta}}\bar{{\it\theta}}}Q_{a}^{\prime }+g_{\bar{{\it\theta}}\bar{{\it\zeta}}}S_{a}^{\prime }-\bar{J}^{\prime }(a)-{\it\nu}_{\bar{{\it\theta}}}^{\prime }=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle G_{p} & = & \displaystyle -Qr_{\bar{{\it\zeta}}}^{\prime }r_{a\bar{{\it\theta}}}^{\prime \prime }+Sr_{\bar{{\it\zeta}}}^{\prime }r_{a\bar{{\it\zeta}}}^{\prime \prime }-Qz_{\bar{{\it\zeta}}}^{\prime }z_{a\bar{{\it\theta}}}^{\prime \prime }+Sz_{\bar{{\it\zeta}}}^{\prime }z_{a\bar{{\it\zeta}}}^{\prime \prime }-Qr^{2}{\it\varphi}_{\bar{{\it\zeta}}}^{\prime }{\it\varphi}_{a\bar{{\it\theta}}}^{\prime \prime }+Sr^{2}{\it\varphi}_{\bar{{\it\zeta}}}^{\prime }{\it\varphi}_{a\bar{{\it\zeta}}}^{\prime \prime }\nonumber\\ \displaystyle & & \displaystyle +\,[(Qr_{\bar{{\it\theta}}}^{\prime }-Sr_{\bar{{\it\zeta}}}^{\prime })_{\bar{{\it\zeta}}}^{\prime }+2r(S{\it\varphi}_{\bar{{\it\zeta}}}^{\prime }-Q{\it\varphi}_{\bar{{\it\theta}}}^{\prime }){\it\varphi}_{\bar{{\it\zeta}}}^{\prime }]r_{a}^{\prime }+(Qz_{\bar{{\it\theta}}}^{\prime }-Sz_{\bar{{\it\zeta}}}^{\prime })_{\bar{{\it\zeta}}}^{\prime }z_{a}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,(r^{2}Q{\it\varphi}_{\bar{{\it\theta}}}^{\prime }-r^{2}S{\it\varphi}_{\bar{{\it\zeta}}}^{\prime })_{\bar{{\it\zeta}}}^{\prime }{\it\varphi}_{a}^{\prime }-g_{{\it\theta}{\it\zeta}}Q_{a}^{\prime }+g_{{\it\zeta}{\it\zeta}}S_{a}^{\prime }-\bar{F}^{\prime }(a)-{\it\nu}_{{\it\zeta}}^{\prime }=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle rD-J=r(r_{\bar{{\it\theta}}}^{\prime }{\it\varphi}_{\bar{{\it\zeta}}}^{\prime }-r_{\bar{{\it\zeta}}}^{\prime }{\it\varphi}_{\bar{{\it\theta}}}^{\prime })z_{a}^{\prime }+r(z_{\bar{{\it\zeta}}}^{\prime }{\it\varphi}_{\bar{{\it\theta}}}^{\prime }-z_{\bar{{\it\theta}}}^{\prime }{\it\varphi}_{\bar{{\it\zeta}}}^{\prime })r_{a}^{\prime }+r(r_{\bar{{\it\zeta}}}^{\prime }z_{\bar{{\it\theta}}}^{\prime }-r_{\bar{{\it\theta}}}^{\prime }z_{\bar{{\it\zeta}}}^{\prime }){\it\varphi}_{a}^{\prime }-J=0. & \displaystyle\end{eqnarray}$$
These equations, which are linear with respect to the radial derivatives of the unknown functions, represent essentially the standard form resolved with respect to the first-order radial derivatives. Otherwise, the equations are nonlinear and the most efficient way to solve them is linearization
$$\begin{eqnarray}\displaystyle r\rightarrow r+{\it\delta}r,\quad z\rightarrow z+{\it\delta}r,\quad {\it\varphi}\rightarrow {\it\varphi}+{\it\delta}{\it\varphi},\quad Q\rightarrow Q+{\it\delta}Q,\quad S\rightarrow S+{\it\delta}S, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
which is straightforward because most of the terms contain only quadratic nonlinearities of 3-D functions. The conversion of (4.8)–(4.10) to LEE can be performed in the same way as in the 2-D case of § 3 or automatically by some software in the process of programming the equations. Because of the numerous terms involved, and with no particular additional specifics, we do not present here the resulting LEE for the 3-D case.
The boundary conditions represent a nonlinear problem
$$\begin{eqnarray}\displaystyle r=R({\it\omega},{\it\varphi}),\quad r=Z({\it\omega},{\it\varphi}) & & \displaystyle\end{eqnarray}$$
for conversion of
${\it\omega},{\it\varphi}$
to
$\bar{{\it\theta}},\bar{{\it\zeta}}$
.
The fundamental solution to the 3-D LEE can be presented as
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}r=\mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}{\it\delta}r^{kk^{\prime }}(a,\bar{{\it\theta}},\bar{{\it\zeta}})=\mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}\mathop{\sum }_{m,n}{\it\delta}r_{mn}^{kk^{\prime }}\text{e}^{-\text{i}(m\bar{{\it\theta}}-n\bar{{\it\zeta}})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}z=\mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}{\it\delta}z^{kk^{\prime }}(a,\bar{{\it\theta}},\bar{{\it\zeta}})=\mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}\mathop{\sum }_{m,n}{\it\delta}z_{mn}^{kk^{\prime }}\text{e}^{-\text{i}(m\bar{{\it\theta}}-n\bar{{\it\zeta}})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}{\it\varphi}=\mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}{\it\delta}{\it\varphi}^{kk^{\prime }}(a,\bar{{\it\theta}},\bar{{\it\zeta}})=\mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}\mathop{\sum }_{m,n}{\it\delta}{\it\varphi}_{mn}^{kk^{\prime }}\text{e}^{-\text{i}(m\bar{{\it\theta}}-n\bar{{\it\zeta}})}, & \displaystyle\end{eqnarray}$$
where
$C_{kk^{\prime }}$
are free parameters and
$\{{\it\delta}r(a,\bar{{\it\theta}},\bar{{\it\zeta}}),{\it\delta}z(a,\bar{{\it\theta}},\bar{{\it\zeta}}),{\it\delta}{\it\varphi}(a,\bar{{\it\theta}},\bar{{\it\zeta}})\}^{k,k^{\prime }}$
are independent sets of 3-D solutions.
The linearized boundary conditions
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}{\it\delta}r^{kk^{\prime }}(a,\bar{{\it\theta}},\bar{{\it\zeta}})=R({\it\omega},{\it\varphi})-r^{b}(\bar{{\it\theta}},\bar{{\it\zeta}})-r_{\bar{{\it\theta}}}^{\prime b}{\it\delta}\bar{{\it\theta}}-r_{\bar{{\it\zeta}}}^{\prime b}{\it\delta}\bar{{\it\zeta}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}{\it\delta}z^{kk^{\prime }}(a,\bar{{\it\theta}},\bar{{\it\zeta}})=Z({\it\omega},{\it\varphi})-z^{b}(\bar{{\it\theta}},\bar{{\it\zeta}})-z_{\bar{{\it\theta}}}^{\prime b}{\it\delta}\bar{{\it\theta}}-z_{\bar{{\it\zeta}}}^{\prime b}{\it\delta}\bar{{\it\zeta}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{k,k^{\prime }}C_{kk^{\prime }}{\it\delta}{\it\varphi}^{kk^{\prime }}(a,\bar{{\it\theta}},\bar{{\it\zeta}})={\it\varphi}-{\it\varphi}^{b}(\bar{{\it\theta}},\bar{{\it\zeta}})-{\it\varphi}_{\bar{{\it\theta}}}^{\prime b}{\it\delta}\bar{{\it\theta}}-{\it\varphi}_{\bar{{\it\zeta}}}^{\prime b}{\it\delta}\bar{{\it\zeta}} & \displaystyle\end{eqnarray}$$
determine the coefficients
$C_{kk^{\prime }}$
and perturbation of the angles
${\it\delta}\bar{{\it\theta}},{\it\delta}\bar{{\it\zeta}}$
, given the initial guess
$r^{b},z^{b},\bar{{\it\theta}},\bar{{\it\zeta}}$
.
This consideration shows that for both 2- and 3-D equilibrium configurations, the LEE have a transparent algorithm of calculations equivalent to a fast Newton scheme.
It should be noted that in the 3-D case, arbitrary plasma boundary shapes cannot be reproduced by fundamental solutions to the LEE. By their nature, the LEE eliminate the boundary perturbations which would generate the magnetic islands on resonant magnetic surfaces. The author is currently unable to formulate this general conclusion in an explicit form, leaving it for future studies.
5. Summary
This paper introduces a rigorous approach for the direct implementation of the Hamada principle in designing 3-D equilibrium configurations with nested magnetic surfaces. The deficiency of the Kruskal–Kulsrud principle is commented on. This paper complements Zakharov and Li (Reference Zakharov and Li2015), where the energy principle and numerical method are described for calculation of 3-D equilibria with magnetic islands using reference magnetic coordinates and Fourier transforms instead of field-line tracing. The approaches of both papers lead to fast Newton numerical schemes for solving nonlinear equilibrium equations.
The linearized equilibrium equations introduced by this paper are especially suitable for numerical implementation. Their solution is given in terms of cylindrical coordinates
$r,{\it\varphi},z$
and their first derivatives as functions of flux coordinates
$a,\bar{{\it\theta}},\bar{{\it\zeta}}$
. This information is sufficient to generate the metric tensor for the next iteration without any additional interpolation or differentiation. For both 2- and 3-D cases, such a property promises very fast and stable calculations. In the 2-D Grad–Shafranov case, the LEE codes could be suitable for the real time forecast of tokamak discharges when the fast equilibrium reconstruction and transport codes work in real time together for prediction of plasma parameters on the energy confinement time scale.
Acknowledgement
This work was partially supported by US DoE contract no. DE-AC02-09-CH11466.
