1. Introduction
The Grad–Shafranov equation is a semilinear partial differential equation (PDE) of second order, which describes the equilibrium state of toroidally symmetric plasmas in the context of ideal magnetohydrodynamics (Grad & Rubin Reference Grad and Rubin1958; Shafranov Reference Shafranov1963). External heating methods such as neutral beam injection, ion-cyclotron resonance heating, etc., induce strong pressure anisotropy in plasmas (Hole et al. Reference Hole, von Nessi, Fitzgerald, McClements and Svensson2011). These contributions are not reflected in the Grad–Shafranov equation, since ideal magnetohydrodynamics assumes an isotropic pressure throughout the plasma. The inclusion of anisotropic effects in equilibrium calculations based on the Grad–Shafranov equation is done through a modification of the plasma pressure tensor that appears in the PDE, changing the equation from being semilinear to fully quasilinear. This approach typically utilises drift-kinetic theory to derive the required plasma pressure tensor components, and has been studied by various authors in the plasma literature (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956; Dobrott & Greene Reference Dobrott and Greene1970; Freidberg Reference Freidberg1982; Iacono et al. Reference Iacono, Bondeson, Troyon and Gruber1990; Fitzgerald, Appel & Hole Reference Fitzgerald, Appel and Hole2013; Qu, Fitzgerald & Hole Reference Qu, Fitzgerald and Hole2014) in both analytic and computational contexts.
There are a plethora of codes to solve the Grad–Shafranov equation and its extensions to include anisotropic pressure in both fixed- and free-boundary cases, such as EFIT (Zwingmann, Eriksson & Stubberfield Reference Zwingmann, Eriksson and Stubberfield2001), FINESSE (Beliën et al. Reference Beliën, Botchev, Goedbloed, van der Holst and Keppens2002), EFIT-TENSOR (Fitzgerald et al. Reference Fitzgerald, Appel and Hole2013), HELENA+ATF (Qu et al. Reference Qu, Fitzgerald and Hole2014), FLOW (Guazzotto & Hameiri Reference Guazzotto and Hameiri2014), SPIDER (Ivanov et al. Reference Ivanov, Martynov, Medvedev and Poshekhonov2015) and HPE (Gorelenkov & Zakharov Reference Gorelenkov and Zakharov2018), among many others. However, literature on the well-posedness of the associated boundary value problem is minimal, with many analysis results focused on simply showing the ellipticity of the problem. Some authors show existence results comprehensively, such as Fujii & Hirai (Reference Fujii and Hirai1983), but this is done only for particular problems such as the isotropic Grad–Shafranov equation in the free-boundary case.
In this work, we focus on deriving existence and uniqueness statements for both the isotropic and anisotropic forms of the Grad–Shafranov equation. We show the existence of weak solutions to the isotropic Grad–Shafranov equation through a new fixed-point method approach, as opposed to the minimisation methods present in the existing literature, and derive conditions for uniqueness. We obtain these new results through the application of well-known functional analysis tools and results to the boundary value problem of the Grad–Shafranov equation. For the anisotropic case, we study the conditions which are required for the ellipticity of the problem, and then show the existence of solutions. These results are of importance since they are key to identifying whether the associated computational problem is well-posed, and also to identifying any required conditions on this well-posedness. Considering such conditions when approaching the problem computationally can be vital to the stability of a code, and also to obtain meaningful results.
The paper is structured as follows. Section 2 lays the physical and mathematical foundations for the problem, describing the Grad–Shafranov problem in both isotropic and anisotropic settings, and also details some useful analysis results of which we make use in later sections. Section 3 presents the main analysis of the isotropic Grad–Shafranov equation, detailing the derived existence and uniqueness results. Section 4 focuses on the anisotropic version of the Grad–Shafranov equation, examining the conditions required for its ellipticity and for the existence of solutions.
2. Preliminaries
2.1. Force balance in axisymmetry
The Grad–Shafranov equation describes a plasma in equilibrium in a toroidal domain symmetric under rotation around the vertical axis (axisymmetric). Considering a cylindrical coordinate system $(R,\varphi ,Z)$
, the toroidal angle $\varphi$
is an ignorable variable. The Grad–Shafranov equation can be derived starting with the force-balance condition of plasmas, $\boldsymbol {j}\times \boldsymbol {B}=\boldsymbol {\nabla } p$
, where $\boldsymbol {B}$
is the magnetic field, $\boldsymbol {j}=\boldsymbol {\nabla }\times \boldsymbol {B}/\mu _0$
the plasma current density, $\mu _0$
the vacuum permeability and $\boldsymbol {\nabla } p$
the pressure gradient, and exploiting the axisymmetry of the tokamak configuration. Using the divergence-free properties of both $\boldsymbol {B}$
and $\boldsymbol {j}$
in this case, we can introduce two corresponding stream functions: $\varPsi$
the poloidal flux (divided by $2{\rm \pi}$
) and $F$
the poloidal current (divided by $2{\rm \pi}$
). Applying Ampère's law and the force-balance condition, and considering the toroidal and poloidal components of each, the magnetic field can be shown to be $\boldsymbol {B} = F(\varPsi )\boldsymbol {\nabla } \varphi + \boldsymbol {\nabla } \varPsi \times \boldsymbol {\nabla } \varphi$
, where $F=F(\varPsi )$
is necessarily a function of the poloidal flux. We also obtain that the pressure $p=p(\varPsi )$
is similarly a flux function. From Ampère's law, we obtain a relationship between derivatives of $\varPsi$
and the toroidal current density $j_\varphi$
, which combined with the relationship between $j_\varphi$
and the flux functions from force balance leads to the Grad–Shafranov equation, presented below.
2.2. Grad–Shafranov equation and pressure anisotropy
By toroidal axisymmetry, this is a two-dimensional problem and we note that in all following sections, we are assuming a two-dimensional domain $\varOmega \subset \{(R,Z)\in \mathbb {R}^2\,|\,R>0\}$
. All operators used are two-dimensional unless explicitly stated otherwise, including the gradient $\boldsymbol {\nabla } = (\partial /\partial R, \partial /\partial Z)$
and Laplacian $\varDelta = \partial ^2/\partial R^2 + \partial ^2/\partial Z^2$
. The Grad–Shafranov equation is expressed as follows (Grad & Rubin Reference Grad and Rubin1958; Shafranov Reference Shafranov1963):
where $\varOmega$
is the cross-sectional region of the torus, and is generally simply considered to be an open, bounded and connected subset of the half-plane; the poloidal current $F:\mathbb {R}\to \mathbb {R}$
and the scalar pressure $p:\mathbb {R}\to \mathbb {R}$
are flux functions (functions of the poloidal flux $\varPsi$
only). This equation is semilinear, due to the nonlinear right-hand side which depends on the flux $\varPsi$
, and elliptic everywhere in the domain $\varOmega$
. This ellipticity can be seen from the fact that the coefficients of the second-order derivatives form a positive-definite matrix (when only considering the highest order of derivatives, the Grad–Shafranov operator is the same as the two-dimensional Laplacian operator; the differences appear in lower-order derivative terms). It is interesting to note that the Grad–Shafranov operator cannot be reduced to a Laplace-Beltrami type operator through any change of metric or coordinates; this is demonstrated in Appendix D.
We close the system with Dirichlet boundary conditions on $\partial \varOmega$
such that
which can be chosen as being homogeneous without loss of generality. The functions $F$
and $p$
are determined by the experimental configuration being studied, and must be given for the system to be fully defined.
While a useful model for plasma shaping and tokamak control, the Grad–Shafranov equation (2.1) is fairly limited to describe burning plasmas. Among others, one important limitation is the fact that the pressure is treated as a scalar quantity.
Heating methods such as neutral-beam injection or ion cyclotron resonance heating can result in particles with differing parallel and perpendicular velocities, and lead to pressure anisotropy – the components of the diagonalised pressure tensor are not equal (Chew et al. Reference Chew, Goldberger and Low1956; Qu et al. Reference Qu, Fitzgerald and Hole2014). The standard way to account for pressure anisotropy is to compute kinetic corrections to plasma fluid models from the guiding-centre formalism of single-particle motion (Dobrott & Greene Reference Dobrott and Greene1970). In this approach, the pressure tensor is presented as
where $\boldsymbol {B}$
is the magnetic field vector and $I$
the identity matrix. The parallel and perpendicular pressure functions, $p_\parallel (R,B,\varPsi )$
and $p_\perp (R,B,\varPsi )$
, respectively, are defined relative to the direction of the magnetic field $\boldsymbol {B}$
.
These two functions now appear in the modified system of equations for the anisotropic problem, instead of the sole flux function for pressure, $p(\psi )$
, as in the isotropic case. It is important to note that the form of these pressure functions must be given (which is typically done through an equation of state, though it is not the only possibility) in order for the system to be closed and solvable. Here, we are considering both pressure anisotropy and toroidal flow; hence the pressures are expressed as functions of three variables instead of just $\psi$
as in the isotropic and stationary case. In the CGL model (Chew et al. Reference Chew, Goldberger and Low1956), these pressures take the form $p_{\parallel } = c_1\rho ^3/B^2$
and $p_\perp = c_2\rho B$
, where $c_1,c_2 \in \mathbb {R}$
are constants and $\rho$
is the plasma density. In the formulation of the Grad–Shafranov solver HELENA+ATF (Qu et al. Reference Qu, Fitzgerald and Hole2014), which includes pressure anisotropy and toroidal flow, the pressures are related to the plasma density through two temperature functions as $p_\parallel = \rho T_\parallel (\varPsi )$
and $p_\perp = \rho T_\perp (\varPsi ,B)$
. The plasma density is fixed through the following definition:
where $H=H(\varPsi )$
denotes an enthalpy flux function, $\varGamma = \varGamma (\varPsi )$
is a rotation flux function and $\varTheta = \varTheta (\varPsi )$
is an anisotropy flux function which relates the perpendicular temperature $T_\perp$
to the parallel temperature $T_\parallel(\varPsi)$
, which is chosen as being a flux function to provide a thermal closure. An initial density $\rho _0 \in \mathbb {R}$
is also fixed. In this formulation, the functions $T_\parallel ,H,\varTheta ,\varGamma$
must be given in order for the system to be solvable.
A more general version of the Grad–Shafranov PDE can be derived to include pressure anisotropy in general, and toroidal flow, by considering the force-balance equation and taking its components in the $\boldsymbol {\nabla } \varphi$
, $\boldsymbol {\nabla }\psi$
, $\boldsymbol {\nabla } R$
and $\boldsymbol {\nabla } B$
directions (Dobrott & Greene Reference Dobrott and Greene1970; Iacono et al. Reference Iacono, Bondeson, Troyon and Gruber1990), resulting in the following form:
where the operator $\boldsymbol {\nabla }_{3D}$
denotes the three-dimensional gradient operator. The function $\sigma (p_\parallel ,p_\perp ,B):\mathbb {R}_+^3 \to \mathbb {R}$
and the toroidal component of the magnetic field strength $B_\varphi$
are defined through the following relations:
for a flux function $F(\varPsi ):\mathbb {R}\to \mathbb {R}$
. Here, the function $\sigma$
represents the strength of the anisotropy,Footnote 1 and (2.5) reduces to the standard isotropic Grad–Shafranov equation as shown in (2.1) when we have $\sigma = 0$
. Equation (2.7) relates the toroidal flux function $F(\varPsi )$
to the toroidal magnetic field strength $B_\varphi$
. The total magnetic field strength, $B$
, is defined through axisymmetry by
where $B_\varphi$
now has a form different from that of the isotropic case, as shown in (2.7). In addition, we have two constraints on the derivatives of the parallel pressure as follows:
where $\varGamma (\varPsi ): \mathbb {R} \to \mathbb {R}$
denotes the toroidal angular velocity, which is non-zero when also considering (purely) toroidal flow. It is important to note that the functions $p_\parallel ,p_\perp ,F(\varPsi ),\varGamma (\varPsi )$
must be given to ensure the solvability of the system.
The modified Grad–Shafranov equation is a quasilinear PDE, since the operator itself now depends on the flux $\psi$
, but further complexity arises compared with the isotropic case, since the anisotropy factor $\sigma$
can lead to the equation transitioning from being elliptic to hyperbolic. This is due to the matrix of second-derivative coefficients no longer necessarily being positive definite, which is the requirement to have an elliptic operator. This is explicitly shown in proposition C.5.
2.3. Definitions and useful theorems
We now give some definitions, and list useful results from functional analysis, which we will use to prove the existence of solutions to the various forms of the Grad–Shafranov equation. Firstly, we define the function spaces in which we are interested.
Definition 2.1 (Square-integrable functions, $L^2$
)
Let $\varOmega$
be a bounded open subset of $\mathbb {R}^n$
, and let $f:\varOmega \to \mathbb {R}$
be a measurable function. Then $f \in L^2(\Omega )$
, or is square-integrable, if
where we identify functions $f$
which agree almost everywhere. The operation $\lVert \cdot \rVert _2$
is the norm on this space.
Definition 2.2 (Sobolev space of weakly differentiable functions, $H^1$
)
Suppose $\varOmega$
and $f$
are as before. We say that $f \in H^1(\varOmega )$
, or is once-weakly differentiable, if $f \in L^2(\Omega )$
and there exist functions $g_1,\ldots ,g_n \in L^2(\Omega )$
such that
where $C^\infty _c(\varOmega )$
denotes the space of infinitely differentiable functions with compact support on the domain $\varOmega$
. The functions $\{g_i\}$
are then referred to as the weak derivatives of $f$
, and we use the notation $\boldsymbol {\nabla } f = (g_1, \ldots , g_n)$
.
We are especially interested in the subspace $H_0^1(\varOmega )$
, which is the closure of the space $C^\infty _c(\varOmega )$
in $H^1(\varOmega )$
.Footnote 2 Functions in this space include those which have compact support and zero traceFootnote 3 on the boundary $\partial \varOmega$
. The norm operation $\lVert \cdot \rVert _{H^1(\varOmega )}$
for functions $f \in H^1(\varOmega )$
is defined as
The following is one of the main theorems for the purpose of showing existence of solutions.
Theorem 2.3 (Schauder's fixed-point theorem) Let $\mathfrak {B}$
be a Banach space, and $M$
be a closed, convex and non-empty subset of $\mathfrak {B}$
. If $T: M \to M$
is a compact operator, then $T$
has a fixed point in $M$
.
Proof. Schauder's fixed-point theorem is proved by extension from Brouwer's fixed-point theorem – we refer the reader to the chapter of Gilbarg & Trudinger (Reference Gilbarg and Trudinger1998) on topological fixed-point theorems for a comprehensive proof.
We also make use of the following inequality in the construction of bounds.
Lemma 2.4 (Poincaré's inequality) Let $\varOmega$
be an open bounded set as before. There exists a constant $C = C(\varOmega )$
such that
for all functions $u \in H_0^1(\varOmega )$
.
Proof. The proof of Poincaré's inequality can be found in any standard functional analysis textbook; we refer the reader to chapter 9 of Brezis (Reference Brezis2010). We note that this is the specific $H^1_0$
version of a more general inequality which holds for Sobolev spaces.
3. Results for the isotropic Grad–Shafranov equation
The severe nonlinearity of the Grad–Shafranov equation, especially in the case when modified for anisotropy and flow, presents some challenges in deriving the well-posedness of the problem. We first present the analysis of the isotropic problem as shown in (3.1), which can be obtained from the modified Grad–Shafranov equation in (2.5) by setting the anisotropy factor to be $\sigma = 0$
. We also note that we are using natural units in this work and will drop physical constants such as $\mu _0$
for the analysis by setting them to have the value of unity.Footnote 4
The isotropic Grad–Shafranov equation is still nonlinear in general, since the functions $F(\varPsi )$
and $p'(\varPsi )$
typically have complex dependencies on the flux $\varPsi$
, but elliptic everywhere in the domain $\varOmega$
with no transition points. We fix homogeneous Dirichlet boundary conditions for the problem to be defined as follows:
Definition 3.1 Let $\varOmega$
be an open bounded domain such that its closure is contained in the half-plane $\mathbb {R}_+\times \mathbb {R}$
, denoted $\varOmega \subset \subset \mathbb {R}_+\times \mathbb {R}$
, and let $\partial \varOmega$
be its boundary. The isotropic problem is to find a function $\psi (R,Z) \in C^2(\varOmega )$
such that
We note the elliptic Grad–Shafranov operator effectively corresponds in $\mathbb {R}^2$
to $\varDelta ^* \varPsi = \varDelta \varPsi - \boldsymbol {\nabla } \ln R\boldsymbol {\cdot }\boldsymbol {\nabla } \varPsi$
, where the two-dimensional Laplacian $\varDelta =\partial _R^2+\partial _Z^2$
and the two-dimensional gradient $\boldsymbol {\nabla } \psi =(\partial _R \psi ,\partial _Z \psi )$
are used throughout.
In this section, we show existence of weak solutions – functions $\psi$
in the Sobolev space $H_0^1(\varOmega )$
which satisfy the weak form of the boundary value problem in (3.1), as follows:
for all test functions $\phi \in H_0^1(\varOmega )$
.
Define the bilinear form $a:H_0^1(\varOmega )\times H_0^1(\varOmega )\to \mathbb {R}$
via $a(u,v)=(\boldsymbol {\nabla } u,\boldsymbol {\nabla } v)_{L^2}+(\boldsymbol {\nabla } \ln R \boldsymbol {\cdot }\boldsymbol {\nabla } u, v)_{L^2}$
. This induces a map $\varDelta _0^*:H_0^1(\varOmega )\to H^{-1}(\varOmega )$
defined by $\varDelta _0^*(u)(v)=a(u,v),\ \forall \ v\in H_0^1(\varOmega )$
.
For convenience, we define the continuous function $f:\mathbb {R}_+\times \mathbb {R}\to \mathbb {R}$
via $f(R,u) = FF'(u) + R^2 p'(u)$
, as well as the function $\mathcal {F}_\psi :\varOmega \to \mathbb {R}$
via $\mathcal {F}_\psi (R,Z)=f(R,\psi (R,Z))$
for some fixed $\psi$
. We consider $\mathcal {F}:H^1_0(\varOmega )\to L^2(\Omega ) \subset H^{-1}(\varOmega )$
eventually to be a map between function spaces, and we show that this map is a continuous map between these spaces in lemma A.2 (Appendix A).
As standard in nonlinear analysis, we make use of the fact that these solutions are fixed points of the associated nonlinear map (Evans Reference Evans2010; Rupflin Reference Rupflin2017) given by
The linear operator $(-\varDelta ^* _0)^{-1}:H^{-1}(\varOmega )\to H_0^1(\varOmega )$
is the inverse of the Grad–Shafranov operator, and the subscript $0$
is used to remind the reader of the Dirichlet boundary conditions.
The strategy to prove existence of solutions to this problem is to show that the operator $A$
is continuous and compact in order to apply Schauder's fixed-point theorem (theorem 2.3). We do this in two parts: first, we show that the inverse operator restricted to $L^2$
functions, $(-\varDelta _0^*)^{-1}|_{L^2(\varOmega )}$
, is compact by considering a reduced problem; we then use the preceding result to obtain the compactness of $A$
, and derive the other necessary conditions in order to obtain solutions to the full problem through Schauder's fixed-point theorem.
First, we show that the inverse Grad–Shafranov operator is compact. To do this, we consider a reduced problem in the form of a modified Dirichlet–Poisson problem by fixing the right-hand side function to be some $\rho \in L^2(\varOmega )$
, for the same Dirichlet boundary conditions. The resultant boundary value problem is almost equivalent to the standard Dirichlet–Poisson problem for the Laplacian operator, albeit with the Grad–Shafranov operator instead:
By imposing a few extra conditions on the class of the right-hand side function and the boundary, we can obtain a useful characterisation of the inverse operator $(-\varDelta _0^*)^{-1}$
as follows.
Lemma 3.2 Let the boundary $\partial \varOmega$
be of class $C^2$
. Then, there exists a unique inverse operator to $\varDelta _0^*$
of the modified Dirichlet–Poisson problem which is continuous as a map $\left (-\Delta _0^*\right )^{-1}|_{L^2(\Omega )}:L^2(\Omega )\to {H_0^1(\varOmega ) \cap H^2(\varOmega )}$
.
Proof. The proof is in two parts – the first part uses the Lax–Milgram theorem (Lax & Milgram Reference Lax and Milgram1954) to obtain an existence and uniqueness statement, and the second involves the application of a regularity result to obtain higher-order differentiability of the solution.
We have that the operator $\varDelta _0^*$
is coercive with respect to the norm on $H_0^1(\varOmega )$
; see lemma B.3 for a proof. We also have that it is continuous as a map from $H_0^1(\varOmega )\to H^{-1}(\varOmega )$
, shown in lemma B.4. By the Lax–Milgram theorem, there exists a unique element $\psi \in H_0^1(\varOmega )$
such that $\psi$
satisfies the weak form of the modified Dirichlet–Poisson problem for a given linear functional $\rho \in H^{-1}(\varOmega )$
. We also obtain that the inverse operator $(-\varDelta _0^*)^{-1}$
is continuous as a map from $H^{-1}(\varOmega )$
into $H_0^1(\varOmega )$
.
By the global regularity result for linear Dirichlet problems (theorem 8.12 of Gilbarg & Trudinger (Reference Gilbarg and Trudinger1998), pp. 186–187), if $\rho \in L^2(\varOmega )$
, we have that the solution $\psi =-(\varDelta _0^*)^{-1}\rho \in H_0^1(\varOmega )$
is not only once but twice weakly differentiable, i.e. $\psi \in {H_0^1(\varOmega ) \cap H^2(\varOmega )}$
, since the boundary $\partial \varOmega$
is assumed to be of class $C^2$
. Thus, we can consider the inverse map as $\left (-\Delta _0^*\right )^{-1}|_{L^2(\Omega )}: L^2(\Omega ) \to {H_0^1(\varOmega ) \cap H^2(\varOmega )}$
, and this is a continuous map due to the Lax–Milgram result.
We will also make use of the consequence of the following lemma.
Lemma 3.3 The inverse Grad–Shafranov operator $\iota \circ \left (-\Delta _0^*\right )^{-1}|_{L^2(\Omega )}:L^2(\Omega )\to H_0^1(\varOmega )$
is a compact map, where $\iota : H_0^1(\varOmega ) \cap H^2(\varOmega ) \hookrightarrow H_0^1(\varOmega )$
denotes the Sobolev embedding.
Proof. Lemma 3.2 shows that the inverse Grad–Shafranov operator is a continuous map from $L^2(\Omega )$
into ${H_0^1(\varOmega ) \cap H^2(\varOmega )}$
. The Rellich–Kondrachov embedding theorem (see theorem 6.3 of Adams & Fournier (Reference Adams and Fournier1975), p. 168) states that the space ${H_0^1(\varOmega ) \cap H^2(\varOmega )}$
is continuously and compactly embedded in the space $H_0^1(\varOmega )$
, provided that the domain $\varOmega$
is a bounded subset of $\mathbb {R}^n$
and the boundary is locally Lipschitz continuous. This requirement is satisfied automatically since $\partial \varOmega \in C^2$
, and so the inverse Grad–Shafranov operator can be considered as a map $\iota \circ \left (-\Delta _0^*\right )^{-1}|_{L^2(\Omega )}:L^2(\Omega )\to H_0^1(\varOmega )$
which is compact.
With these results, we now show the existence of solutions to the isotropic Grad–Shafranov equation.
Theorem 3.4 Let $\varOmega \subset \subset \mathbb {R}_+\times \mathbb {R}$
be an open bounded set, and let the boundary $\partial \varOmega$
be of class $C^2$
. Let $f:\mathbb {R}_+\times \mathbb {R} \to \mathbb {R}$
be a continuous and bounded function, such that $\sup |f| = c < \infty$
, and let $\mathcal {F}_\psi \in L^2(\varOmega )$
be the function defined before. Then, there exist weak solutions $\psi \in H_0^1(\varOmega )$
to the isotropic Grad–Shafranov equation with right-hand side given by $f(R,\psi )$
.
Proof. We use theorem 2.3 to show the existence of a fixed point, which is equivalent to the existence of a solution. Define the map
which, as per (3.3a,b), is the nonlinear map associated with the Grad–Shafranov equation.
We proceed by constructing the three conditions required for Schauder's theorem. First, we show that $A$
is a continuous map between these function spaces.
Note that $\psi \mapsto \mathcal {F}_\psi$
is a continuous map from $H_0^1(\varOmega )$
into $L^2(\Omega )$
– this is a consequence of the bounded and continuous nature of the function $f$
, and the boundedness of the domain $\varOmega$
. A proof of this continuity can be found in Appendix A. The operator $\iota \circ \left (-\Delta _0^*\right )^{-1}|_{L^2(\Omega )}$
is a continuous map from $L^2(\Omega )$
into $H_0^1(\varOmega )$
, which we proved in lemma 3.2. As the composition of continuous maps is continuous, we have that the operator $A$
is continuous.
Next, we construct a set $M \subset H_0^1(\varOmega )$
which is convex, non-empty and bounded, such that $A:M \to M$
. We can construct a bound on the range of the operator $A$
by using the fact that $A\psi \in H_0^1(\varOmega )$
satisfies the variational form of the PDE:
where $a(\cdot , \cdot )$
is the bilinear form introduced prior, and the test function is chosen as $\phi = A\psi$
. Applying the Cauchy–Schwarz inequality to the right-hand side, we obtain an upper bound as follows:
By lemma B.3, we have that the Grad–Shafranov operator is coercive and so we also have a lower bound on the bilinear form
where the second inequality follows from the bound obtained in (3.7), and the constant $\alpha ^* < \infty$
is the coercivity constant. For convenience, denote $K_1 = c |\varOmega |/\alpha ^*$
, and we have that the operator $A$
is bounded in $L^2(\Omega )$
. To show that it is bounded in $H_0^1(\varOmega )$
, we also need to construct a bound for the gradient. We note that we can expand the left-hand side of (3.7) as a sum of a gradient term and another term as follows:
and since the second term of the left-hand side is strictly positive, we obtain through the use of Poincaré's inequality on the right-hand side that
If we choose $K_2 = c|\varOmega |C(\varOmega )$
, we can now construct a set $M \subset H_0^1(\varOmega )$
such that the image of any bounded subset of $M$
under $A$
is bounded in $H_0^1$
as follows:
where $K^2 = K_1^2 + K_2^2 < \infty$
.
Finally, we show that the operator $A$
is a compact operator on the set $M$
. We note that $A = (- \varDelta ^*)^{-1} \circ \mathcal {F}$
is the composition of two continuous maps, and so it is sufficient to show that $(- \varDelta ^*)^{-1}: L^2(\Omega ) \to H_0^1(\varOmega )$
is compact. The composition of a continuous linear map and a compact map is itself compact.Footnote 5 By lemma 3.3, we have that the inverse operator is compact from $L^2(\Omega )$
into $H_0^1(\varOmega )$
. As such, the map $A = (- \varDelta ^*)^{-1} \circ \mathcal {F}$
is indeed a compact map.
By these points, the conditions for Schauder's fixed-point theorem are satisfied and so we have the existence of a fixed point of $A$
in $M \subset H_0^1$
.
To show uniqueness, we demand an extra constraint on the right-hand side of the PDE in the form of Lipschitz continuity.
Lemma 3.5 Let $\mathcal {F}: H_0^1 \to L^2(\Omega )$
be bounded and Lipschitz continuous, with Lipschitz constant $L \geq 0$
, and let $\varOmega$
be a bounded, convex domain with diameter $d$
. Let the coercivity constant of the Grad–Shafranov operator be denoted $\alpha ^*$
. The isotropic Grad–Shafranov equation with right-hand side given by the map $\mathcal {F}$
has a unique solution, provided that the Lipschitz constant is smaller than the coercivity constant, $L < \alpha ^*$
.
Proof. Let $u,v \in H_0^1(\varOmega )$
be two solutions to the Grad–Shafranov equation for homogeneous Dirichlet boundary value conditions, such that they satisfy the variational form of the PDE
where $\varDelta ^* = \varDelta ^2 - \boldsymbol {\nabla } \ln R \boldsymbol {\cdot }\boldsymbol {\nabla }$
is the two-dimensional linear elliptic operator, and that $u = v = 0$
for all $x \in \partial \varOmega$
. By the linearity of the operator, we can make the following observation on the domain $\varOmega$
:
for any arbitrary function $\phi \in H_0^1$
. Choosing $\phi = u-v$
, we can manipulate this relation as follows:
where the second line follows from the coercivity of the Grad–Shafranov operator $\varDelta ^*$
(which we prove in lemma B.3), the third from the use of the Cauchy–Schwarz inequality and the last simply from the definition of Lipschitz continuity.
To show that any solution which exists is unique, we want to show from this that $\lVert u -v\rVert _{H_0^1}^2 = 0$
, since this would imply that $u = v$
for any functions $u,v \in H_0^1$
that satisfy the PDE. For this we require that $(\alpha ^*-L)\lVert u -v\rVert _{H_0^1}^2 \geq 0$
, as we already know that this quantity is bounded above by zero from (3.14). This implies that we have the condition $\alpha ^* > L$
, i.e. the Lipschitz constant of our right-hand side function is bounded above by the coercivity constant of the Grad–Shafranov operator.
As derived in lemma B.3, the coercivity constant is given by $\alpha ^* = \max \{1,1/2R_{\max }\}$
, where $R_{\max }$
is the maximal radius of the domain $\varOmega$
. Thus, we obtain a sufficient bound on the Lipschitz constant as
which ensures that $(\alpha ^*-L)\lVert u -v\rVert _{H_0^1}^2 \geq 0$
, and as this is bounded above by zero, we obtain $\lVert u -v\rVert _{H_0^1} = 0$
. Thus, any solution which exists must be unique, under the derived condition, and this completes the proof.
4. Results for the anisotropic modified Grad–Shafranov equation
In the previous section, we showed the well-posedness of the homogeneous Dirichlet boundary value problem for the isotropic Grad–Shafranov equation. For the modified Grad–Shafranov equation with anisotropy included, however, it becomes more complicated to show the existence and uniqueness of solutions since the differential operator on the left-hand side is now quasilinear and takes the form
On the left-hand side, we express the operator in the usually adopted divergence form in cylindrical coordinates, and in the two-dimensional Cartesian convention of this paper on the right-hand side. Here, the function $\sigma :\mathbb {R}^3 \to \mathbb {R}$
is given by
where the functions $p_\parallel (R,B,\psi ):\mathbb {R}^3 \to \mathbb {R}$
and $p_\perp (R,B,\psi ): \mathbb {R}^3 \to \mathbb {R}$
are the parallel and perpendicular pressure components, respectively, and $B = |\boldsymbol {B}|$
is the magnetic field strength.
This differential operator is no longer elliptic everywhere, and can transition to being hyperbolic (see proposition C.5 for a proof). Studying where this transition occurs is important, since imposing constraints to ensure the ellipticity of the equation is key to obtaining a well-posed problem in this case. We discuss sufficient conditions to restrict the equation to the elliptic regime in the following section.
4.1. Ellipticity of modified Grad–Shafranov equation
There are two primary conditions which are imposed on the modified Grad–Shafranov boundary value problem to ensure its well-posedness. These are known as the firehose and mirror stability criteria (Grad Reference Grad1967; Qu et al. Reference Qu, Fitzgerald and Hole2014), and are defined as follows.
Definition 4.1 (Firehose stability) Let $\sigma = \sigma (R,Z,\psi ,|\boldsymbol {\nabla } \psi |): \varOmega \times \mathbb {R} \times \mathbb {R}^2 \to \mathbb {R}$
represent the function defined in (4.2). The firehose stability criterion is given by the condition
Definition 4.2 (Mirror stability) Let $p_\perp (\psi ,B) = p_\perp (R,Z,\psi ,|\boldsymbol {\nabla } \psi |): \varOmega \times \mathbb {R} \times \mathbb {R}^2 \to \mathbb {R}$
denote the perpendicular pressure function and $B: \varOmega \times \mathbb {R} \times \mathbb {R}^2 \to \mathbb {R}$
the magnetic field strength. The mirror stability criterion is then defined by the following condition:
We can show that the two stability criteria together are sufficient conditions to ensure the ellipticity of the generalised Grad–Shafranov equation. We show this in the following proposition.
Proposition 4.3 The modified Grad–Shafranov operator, as in (4.1), is elliptic if the firehose and mirror stability criteria are satisfied.
Proof. A second-order differential operator is elliptic if the coefficient matrix of second order derivatives is positive definite. In this case, we obtain the following two eigenvalues (see proposition C.5 for a full derivation):
where $\tau = 1 - \sigma$
and $\dot {\tau }$
is shorthand for the derivative of $\tau$
with respect to $|\boldsymbol {\nabla }\psi |^2$
.
If the firehose stability criterion in (4.3) is satisfied, we immediately have positivity of the first eigenvalue:
To see that the two criteria are sufficient to ensure positivity of the second eigenvalue, we require the specific form of the function $\dot {\tau }$
. Note that the dependency of $\tau$
on $|\boldsymbol {\nabla }\psi |^2$
is through the function $B = \sqrt {B_\varphi ^2 + |\boldsymbol {\nabla }\psi |^2/R^2}$
, and so we obtain the expression for $\dot {\tau }$
as
Substituting this into the expression for the second eigenvalue in (4.5a,b), we have
We know that $B^2 = B_\varphi ^2 + B_p^2 \geq B_\varphi ^2$
, and so it must be that $(1 - B_\varphi ^2/B^2) \geq 0$
. If the firehose stability criterion is satisfied, we also have that $\tau > 0$
. As such, the only remaining term which would guarantee the positivity of $\lambda _2$
gives us the requirement
which is exactly the mirror stability condition from (4.4). Thus, we have shown that the firehose and mirror stability criteria are sufficient conditions to maintain the positivity of the eigenvalues $\lambda _1$
and $\lambda _2$
of the modified Grad–Shafranov operator, and this concludes the proof.
4.2. Existence of solutions
We can derive the existence of solutions to the Dirichlet boundary value problem for the anisotropic Grad–Shafranov equation by restricting it to be elliptic under the firehose and mirror stability criteria described in the previous section. It is, however, still complicated to show this because of the quasilinearity of the operator. As such, we work with the explicit form of the operator in (4.1).
First, we introduce the concept of uniform ellipticity (Gilbarg & Trudinger Reference Gilbarg and Trudinger1998).
Definition 4.4 (Uniform ellipticity) Consider a second-order quasilinear elliptic operator of the form
where $a_{ij} = a_{ji}$
such that the matrix is symmetric, $x \in \varOmega \subset \mathbb {R}^n$
where $\varOmega$
is open and connected and $u \in C^2(\varOmega )$
as standard. Let $\lambda (x,u,Du)$
and $\varLambda (x,u,Du)$
denote the minimum and maximum eigenvalues, respectively, of the coefficient matrix $[a_{ij}]_{n\times n}$
. The operator $Q$
is called uniformly elliptic if its minimum and maximum eigenvalues satisfy
for all $(x,u,Du) \in \varOmega \times \mathbb {R} \times \mathbb {R}^n$
.
Uniform ellipticity is a helpful property when showing the existence of solutions to a quasilinear PDE. It typically reduces the constraints one must impose on the problem to obtain a statement of existence, such as those on the shape of the boundary of the domain (Gilbarg & Trudinger Reference Gilbarg and Trudinger1998). We now show that under certain conditions, the anisotropic Grad–Shafranov operator is uniformly elliptic.
Lemma 4.5 Let $\varOmega \subset \mathbb {R}^2$
be an open and connected set, and let $\sigma : \varOmega \times \mathbb {R} \times \mathbb {R}^2 \to \mathbb {R}$
be a function such that $0 \leq \sigma \leq 1$
for all $(R,Z,\psi ,|\boldsymbol {\nabla }\psi |)$
, so that the firehose stability criterion is satisfied everywhere in the domain. Also suppose that the mirror stability criterion of (4.4) is satisfied everywhere. The generalised Grad–Shafranov operator,
is then uniformly elliptic.
Proof. From proposition 4.3, we know that the operator has two eigenvalues of the form
where we now use $\lambda$
and $\varLambda$
to denote the minimum and maximum eigenvalues, respectively. Their ratio is then given by
Using the expression for $\dot {\tau }$
from (4.7), we have
We note the helpful fact that
which is easily shown to be bounded when normalised by $B^2$
. Thus, we can write (4.16) as
which will be bounded when the following conditions are satisfied:
and so the generalised Grad–Shafranov operator is uniformly elliptic under these conditions.
Remark We note that in the physical systems of interest, the boundedness conditions of (4.19a,b) are always satisfied: $\tau < \infty$
as $\sigma < \infty$
by definition. The second condition is slightly harder to see. We note that the functional dependence of $p_\perp$
typically dictates $\partial p_\perp /\partial B \propto B^2$
, so that
We have that $B$
is bounded, as in $|B| < \infty$
, and so the second condition is satisfied. Thus, we can restrict the study of well-posedness of the Grad–Shafranov equation boundary value problem to the cases where the operator is uniformly elliptic.
Finally, we state a domain shape condition which will be required to show the existence of solutions (Gilbarg & Trudinger Reference Gilbarg and Trudinger1998).
Definition 4.6 (Exterior sphere condition) Consider an open bounded domain $\varOmega \subset \mathbb {R}^n$
. We say that $\varOmega$
satisfies an exterior sphere condition if, for every point $\xi \in \partial \varOmega$
, there exists a ball $B = B_R(y)$
satisfying $\bar {B} \cap \bar {\varOmega } = \xi$
, where $\bar {B}$
denotes the closed ball.Footnote 6
Remark We note that the exterior sphere condition can be related to a requirement that the domain $\varOmega$
contains no non-regular points, such as the tip of a cusp.
We can now introduce the existence theorem for the case of uniformly elliptic equations for two-dimensional domains.
Theorem 4.7 Let $\varOmega$
be a domain in $\mathbb {R}^2$
satisfying an exterior sphere condition, and let $\varphi$
be a continuous function on $\partial \varOmega$
. Then, if $Q$
is an elliptic quasilinear operator such that $Qu = au_{xx} + 2bu_{xy} + cu_{yy} + f(x,y,u,u_x,u_y)$
, a solution $u \in C^2(\varOmega ) \cap C^0(\bar {\varOmega })$
to the Dirichlet boundary value problem
exists under the following conditions:
(i) The functions $a(x,y,u,u_x,u_y), \ldots ,f(x,y,u,u_x,u_y)$
are defined for all $(x,y,u, u_x,u_y) \in \varOmega \times \mathbb {R} \times \mathbb {R}^2$, and are at least $C^0$ continuous on this domain.(ii) The operator $Q$
is uniformly elliptic in $\varOmega$ for bounded functions $u$, according to definition 4.4.(iii) The function $f$
satisfies structure conditions of the form
(4.22)\begin{gather} \frac{|f|}{\lambda} \leq \mu(|u|)(1 + |u_x| + |u_y|), \end{gather}(4.23)\begin{gather}\frac{f}{\lambda} \textrm{sign} (u) \leq \nu(1 + |u_x| + |u_y|), \end{gather}where $\mu (|u|)$ is a non-decreasing function and $\nu$ is a non-negative constant.
Proof. We refer the reader to Gilbarg & Trudinger (Reference Gilbarg and Trudinger1998) for the proof of the theorem.
By theorem 4.7, we can now show the existence of solutions to the boundary value problem of the modified Grad–Shafranov equation in two dimensions.
Theorem 4.8 Let $\varOmega \subset \subset \mathbb {R}_+\times \mathbb {R}$
be an open bounded domain satisfying the exterior sphere condition,Footnote 7 and let $F(\psi )$
, $p_\parallel (R,B,\psi )$
and $p_\perp (R,B,\psi )$
be bounded and continuous functions with bounded and continuous first and second derivatives, where $B = \sqrt {B_\varphi ^2 + |\boldsymbol {\nabla } \psi |^2/R^2} > 0$
, such that the firehose and mirror stability criteria of (4.3) and (4.4) are satisfied. The modified Grad–Shafranov equation
with Dirichlet boundary condition
for a constant $\psi _1 < \infty$
has solutions $\psi \in C^2(\varOmega ) \cap C^0(\bar {\varOmega })$
.
Proof. The proof is by application of theorem 4.7. We show that the required conditions for the theorem hold in this case.
We have that the exterior sphere condition for the domain $\varOmega$
is satisfied by construction. By lemma 4.5, we have that the differential operator is uniformly elliptic, which satisfies condition (ii) of the theorem.
To satisfy conditions (i) and (iii), we note the bounded and continuous nature of the functions $F(\psi )$
, $p_\parallel (R,B,\psi )$
and $p_\perp (R,B,\psi )$
, and their derivatives. The domain $\varOmega$
also contains no points of singularity since we do not include the origin, and $B>0$
strictly.
This satisfies the requirements of theorem 4.7, and so we have the existence of a solution $\psi \in C^2(\varOmega ) \cap C^0(\bar {\varOmega })$
for the specified Dirichlet boundary condition.
5. Conclusion
In this work, we have derived analytical results for the Grad–Shafranov equation under two different regimes. We have demonstrated the existence of solutions to the standard (or isotropic) Grad–Shafranov equation through a new fixed-point method approach from nonlinear analysis, and have derived sufficient conditions for the toroidal current functional (the right-hand side of the Grad–Shafranov equation) and the boundary for the uniqueness of solutions. We have also studied the modified Grad–Shafranov equation for pressure anisotropy and toroidal flow, and shown that the well-known firehose and mirror stability conditions are sufficient to maintain the ellipticity of the modified Grad–Shafranov equation in this regime. We detail this and obtain other conditions on the boundary and right-hand side functions in order to show the existence of solutions in this case. These results rigorously connect criteria that are well known in the plasma physics community to the underlying mathematics of the physical problem, and present a comprehensive picture of the conditions required to maintain the existence of solutions. One question that we have not addressed in this work is that of free-boundary solutions. In general, proving the existence of solutions in free-boundary cases is significantly more difficult than their fixed-boundary counterparts. The basic ingredient for such constructions, however, is the existence of solutions to the corresponding fixed-boundary problem. As such, the results of this paper still provide a valuable contribution to obtaining similar results for the free-boundary case and form the foundation for future work in this area. Finally, we also note that the conditions required for unique (fixed-boundary) solutions in the anisotropic regime are a subject of ongoing work.
Acknowledgements
The first author would like to thank Dr K. Duru for helpful discussions which prompted the beginnings of this research.
Editor Per Helander thanks the referees for their advice in evaluating this article.
Funding
This work (Z.S.Q.) was partly funded by Australian Research Council project DP170102606 and was supported by a grant from the Simons Foundation/SFARI (560651, AB).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Continuity of the inverse isotropic Grad–Shafranov operator
To prove the existence of solutions to the isotropic Grad–Shafranov equation, we introduced the nonlinear operator $A = (-\varDelta _0^*)^{-1} \circ \mathcal {F}$
, $A: H_0^1(\varOmega ) \to H_0^1(\varOmega )$
, which is the composition of two maps. A useful fact is that each of the component maps, $(-\varDelta _0^*)^{-1}: L^2(\Omega ) \to H_0^1(\varOmega )$
and $\mathcal {F}:H_0^1(\varOmega ) \to L^2(\Omega )$
, is bounded and continuous. We prove this in the following two lemmas.
Lemma A.1 Let $\varOmega \subset \subset \mathbb {R}_+\times \mathbb {R}$
be an open bounded domain, such that there exists some $R_{\min } > 0$
where $(R_{\min }, Z_{\min }) = \inf (\varOmega )$
. The map defined by the operator $(-\varDelta _0^*)^{-1} = (-\varDelta + \boldsymbol {\nabla } \ln R \boldsymbol {\cdot }\boldsymbol {\nabla })^{-1}$
is continuous from $L^2(\Omega )$
to $H_0^1(\varOmega )$
.
Proof. The operator $(-\varDelta _0^*)^{-1}$
is linear, and so it is sufficient to show that it is bounded to prove continuity. We show that there exists some positive constant $0 \leq M_1 < \infty$
such that
where $v = (-\varDelta + \ln R \boldsymbol {\cdot }\boldsymbol {\nabla })^{-1}u$
. Variationally, we have
for all test functions $\phi \in H_0^1(\varOmega )$
. Choose $\phi = v$
, so
where we have used the Cauchy–Schwarz inequality on the right-hand side to obtain an upper bound. Recognising that $\lVert v \rVert _{L^2(\Omega )} \leq \lVert v \rVert _{H_0^1(\varOmega )}$
on the right-hand side, and using the Poincaré inequality on the left-hand side, we have
where $D \geq 1$
is a constant, completing the proof.
Lemma A.2 Let $f(R,u):\mathbb {R}_+ \times \mathbb {R} \to \mathbb {R}$
be a continuous and bounded function, where $\sup |f| = a < \infty$
. Then, the map defined by $\mathcal {F}^R_u: H_0^1(\varOmega ) \to L^2(\Omega )$
is continuous, where $\mathcal {F}^R_u = f(R,u(R,Z))$
, for $u \in H_0^1(\varOmega )$
and some fixed $R \in \mathbb {R}_+$
.
Proof. First, we show that $\mathcal {F}^R_u \in L^2(\Omega )$
for fixed $R \in \mathbb {R}_+$
and $u \in H_0^1(\varOmega )$
. Consider
where we have used the boundedness of the function $f$
in the second line and the finite measure of the domain $\varOmega$
in the third. As such, we have that $\lVert \mathcal {F}_u^R\rVert _{L^2(\Omega )} < \infty$
and so $\mathcal {F}_u^R \in L^2(\Omega )$
. The continuity of the map $\mathcal {F}_u^R$
follows from the required continuity of $f$
.
Appendix B. Coercivity and continuity of the Grad–Shafranov operator
A key result required to show existence of solutions is the fact that the operator associated with the isotropic Grad–Shafranov PDE is continuous and coercive. We first show the coercivity property in the following lemma.
Lemma B.3 Let $\varOmega$
be as before, and let $R_{\min }$
and $R_{\max }$
be the infimum and supremum of $R$
in $\varOmega$
, respectively. The isotropic Grad–Shafranov operator is coercive on $H_0^1(\varOmega )$
norm.
Proof. Let $\varOmega$
denote an open bounded subset of $\mathbb {R}_+ \times \mathbb {R}$
, and let $(\cdot , \cdot )$
represent the $L^2$
inner product. The isotropic Grad–Shafranov operator can be explicitly written as follows, through the use of differential and vector calculus identities:
where the last line follows by use of integration by parts on the second term. The associated bilinear form for this operator is then given by
where we simplify the second term through repeated uses of integration by parts.
Now, we note that by the bounded nature of the domain, and the exclusion of the origin, we can write
where we use the fact that $R > 0$
. Using this bound, we obtain that
To finally obtain the coercivity condition, we must consider two cases: $1/R_{\max } > 1$
and $1/R_{\max } < 1$
. If the condition $1/R_{\max } > 1$
holds, we immediately have
and the operator is coercive with coercivity constant $\alpha ^* = 1$
. Conversely, if the condition $1/R_{\max } < 1$
holds, we obtain
i.e. the coercivity condition is satisfied with constant $\alpha ^* = 1/R_{\max }$
. Thus, there always exists some constant $\alpha ^* = \max \{1,1/2R_{\max }\} > 0$
such that the coercivity condition is satisfied, and this concludes the proof.
Lemma B.4 The isotropic Grad–Shafranov operator, $\varDelta _0^*:H_0^1(\varOmega ) \to H^{-1}(\varOmega )$
, is a continuous operator.
Proof. To show that the operator is continuous, we instead consider the induced bilinear form $a(u,v)=(\boldsymbol {\nabla } u,\boldsymbol {\nabla } v)_{L^2}+(\boldsymbol {\nabla } \ln R \boldsymbol {\cdot }\boldsymbol {\nabla } u, v)_{L^2}$
, where the map acts as $a: H_0^1(\varOmega ) \times H_0^1(\varOmega ) \to \mathbb {R}$
. Showing that this form is bounded is equivalent to showing that the associated linear operator is continuous. Consider functions $u,v \in H_0^1(\varOmega )$
; an upper bound on the bilinear form can be constructed as follows:
where the second line follows through the use of integral properties and the Cauchy–Schwarz inequality. Noting that the norm on $H_0^1(\varOmega )$
is $\lVert u \rVert _{H_0^1(\varOmega )}^2 = \lVert u \rVert _{L^2(\Omega )}^2 + \lVert \boldsymbol {\nabla } u\rVert _{L^2(\Omega )}^2$
, and that the embedding $H_0^1(\varOmega ) \hookrightarrow L^2(\Omega )$
is continuous, we can construct the bound:
which is the required result. Thus, the bilinear form $a$
is bounded and continuous, and the associated linear operator $\varDelta _0^*:H_0^1(\varOmega ) \to H^{-1}(\varOmega )$
is also continuous.
Appendix C. Anisotropic Grad–Shafranov equation is elliptic-hyperbolic
Proposition C.5 The modified Grad–Shafranov equation belongs to the class of elliptic-hyperbolic equations.
Proof. The proof is by direct computation of the coefficient matrix of second-order derivatives in the differential operator in (4.1). By the axisymmetry of the toroidal domain, the poloidal flux function takes the form $\psi = \psi (R,Z)$
. We can then note that its gradient takes the form
where we use the contracted subscript notation to indicate differentiation with respect to the subscripted variable. For convenience, we introduce the function $\tau$
such that $\tau = 1 - \sigma$
, and so our differential operator can be expanded as
We note that through the dependency of $\sigma$
on the functions $p_\parallel$
, $p_\perp$
and $B$
, the function $\tau$
is dependent on both $\psi$
and $|\boldsymbol {\nabla } \psi |^2$
and introduces the complicated quasilinearity of the operator. We can take derivatives as follows:
where we only require the term containing derivatives of $|\boldsymbol {\nabla }\psi |^2$
since the other terms will contain first-order derivatives at most. Expanding the relevant term, we have
where the dot refers to differentiation with respect to $|\boldsymbol {\nabla }\psi |^2$
. We obtain a similar expression for the derivative of $\tau$
with respect to the $Z$
coordinate:
By these two expressions, the coefficient matrix for second-order derivatives can be written as
The characteristic determinant of this matrix is then
which is not always guaranteed to be positive, in which case the matrix $[ a_{ij}]$
is no longer positive definite and the operator becomes hyperbolic.
Appendix D. Laplace–Beltrami operators on Riemannian manifolds
Let $(M,g)$
be a $D$
-dimensional Riemannian manifold. The gradient of a smooth function $f\in C^\infty (M)$
is the vector field $\text {grad}_g f:M\to TM$
:
The divergence of a vector field $X\in \varGamma (TM)$
is the function $\text {div}_g X : M\to R$
such that
where $\mu _g$
is the natural volume form associated with the metric $g$
and $\mathcal {L}_X=i_X d + d i_X$
is the Lie derivative on forms. It is immediately seen that the following Leibniz property holds:
Combining the gradient and the divergence operators, we have the Laplace–Beltrami operator on functions $f\in C^\infty (M)$
as
It is a second-order linear elliptic operator, which can be expressed locally as
where $g_{ij}=g(\partial _i,\partial _j)$
and $g^{ij}g_{jk}=\delta ^i_k$
.
Two manifolds $(M,g)$
and $(N,h)$
are isometric if there exists a diffeomorphism (smooth map with smooth inverse) $\varPhi :M\to N$
such that the metrics are related by pullback, $g=\varPhi ^* h$
, i.e. $g_p(v,w)=h_q(d\varPhi _p(v),d\varPhi _p(w))$
, where $p\in M$
, $q=\varPhi (p)\in N$
and $v,w\in T_pM$
.
Proposition D.1 The gradient, the divergence and the Laplace–Beltrami operators of $\varPhi$
-isometric Riemannian manifolds are $\varPhi$
-related, namely
for all $\bar {f},\bar {K}\in C^\infty (N)$
and $\bar {X}\in \varGamma (TN)$
.
Proof. Note that $g(X,Y)=\varPhi ^*[h(\varPhi _* X,\varPhi _*Y)]$
, where $X,Y\in \varGamma (TM)$
. We compute
Hence $\text {grad}_g\varPhi ^* \bar {f} = \varPhi ^* \text {grad}_h \bar {f}$
.
The volume form associated with the metric $g$
on $M$
is the pullback $\mu _g = \varPhi ^* \bar {\mu }_h$
of the volume form $\bar {\mu }_h$
related to $h$
on $N$
. Letting $\bar {X}\in \varGamma (TN)$
, we compute
Hence $\text {div}_g(\varPhi ^*\bar {X}) = \varPhi ^*\text {div}_h\bar {X}$
.
By combining these two results, the third claim is shown:
The final result follows as a special case with $\bar {K}=1$
.
Remark Proposition D.6 conveys the fact that a change of coordinates does not alter the nature of the Laplace–Beltrami operator. In other words, identifying a differential operator as the Laplace–Beltrami operator for some metric is a necessary condition to being able to reduce it to the (standard) Laplacian via a change of coordinates.
Consider now a generic diffusion operator (elliptic second-order linear differential operator) on the space of smooth functions on a $D$
-dimensional manifold $M$
, expressed locally in the form
where $a$
is a smooth positive function and $b^{ij}$
are smooth functions forming a symmetric positive-definite matrix at every point.
Proposition D.7 There exists a smooth function $K$
and a metric $g$
such that (D12) is
Proof. Define $K^D = a^{D-2}|b|$
and $g=({|b|}/{a^2})^{1/D} b^{-1}$
. We compute $\sqrt {g} ={|b|}/{a^2|b|}= {1}/{a}$
and
so that $Lf=a\partial _i b^{ij}\partial _j f = ({1}/{\sqrt {g}})\partial _i(K\sqrt {g}g^{ij}\partial _j f) = \text {div}_g (K \text {grad}_g f)$
.
Remark A major advantage of adopting the metric $g$
from proposition D.7 associated with the diffusion operator is that it immediately provides self-adjointness.
Proposition D.8 The function $K$
and metric $g$
of proposition D.7 are unique up to uniform scaling.
Proof. Assume two metrics $g,\tilde {g}$
and two functions $K,\tilde {K}$
represent the same diffusion operator. Then, $\text {div}_{\tilde {g}} \tilde {K} \text {grad}_{\tilde {g}}-\text {div}_g K \text {grad}_g \equiv 0$
. In local coordinates, this means
The term with second-order derivatives reads $(\tilde {K} \tilde {g}^{ij}-K g^{ij})\partial _i\partial _j f$
. Since $f$
is arbitrary, this term must vanish, implying that the metrics are conformal, $\tilde {g}=({\tilde {K}}/{K}) g =\lambda ^2 g$
. Substituting $\sqrt {\tilde {g}} = \lambda ^D \sqrt {g}$
and $\tilde {K}\tilde {g}^{ij}=Kg^{ij}$
, we have
Since $f$
is arbitrary, $\lambda = \tilde {K}/K$
must be constant, which corresponds to equal metrics and functions up to uniform scaling.
Corollary D.4 The diffusion operator $L=a\partial _i b^{ij}\partial _j$
is the Laplace–Beltrami operator associated with some metric if and only if $a^{D-2}|b|$
is constant.
Proof. The function $K$
is uniquely identified as $(a^{D-2}|b|)^{1/D}$
up to a constant factor. Equation (D13) shows that the diffusion operator matches the Laplace–Beltrami operator if and only if $dK=0$
.
D.1. Applications
To highlight the impact of these results, we apply the criterion of corollary D.9 to two important operators that appear in the context of this paper. The first one derives from expressing the three-dimensional Laplacian in cylindrical coordinates and discarding derivatives with respect to the toroidal angle by assuming axisymmetry. The operator in question is
where $g_{ij} = R \delta _{ij}$
. Since $L$
cannot be identified with the Laplace–Beltrami operator (the function $R$
is not constant), there will be no coordinate transform that morphs it into the (standard) Laplacian.
The second case is the Grad–Shafranov operator, obtained by taking the toroidal component of the Laplacian of the vector potential and discarding the derivatives with respect to the toroidal angle by assuming axisymmetry. The operator reads
where $g_{ij}=({1}/{R})\delta _{ij}$
. For the same reason as before, there will be no change of coordinates in which the Grad–Shafranov matches the (standard) Laplacian.
That being said, it is interesting to note that by feeding $f=\sqrt {R} U$
through the Grad–Shafranov operator, one obtains (Pataki et al. Reference Pataki, Cerfon, Freidberg, Greengard and O'Neil2013)
