Quasisymmetric Magnetic Fields in Asymmetric Toroidal Domains

In the context of ideal magnetohydrodynamics with isotropic pressure, steady magnetic confinement of a plasma is achieved by balance between Lorentz force and pressure gradient,

	$\displaystyle\left({\nabla\times\boldsymbol{B}}\right)\times\boldsymbol{B}=\nabla P,~{}~{}~{}~{}\nabla\cdot\boldsymbol{B}=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega,$		(1a)
	$\displaystyle\boldsymbol{B}\cdot\boldsymbol{n}=0~{}~{}~{}~{}{\rm on}~{}~{}\partial\Omega.$		(1b)

In the notation above, $\boldsymbol{B}$ is the magnetic field, $P$ the pressure, $\Omega\subset\mathbb{R}^{3}$ a smooth bounded domain, and $\boldsymbol{n}$ the unit outward normal to the bounding surface $\partial\Omega$. As a system of nonlinear first order partial differential equations for the Cartesian components $B_{x}$, $B_{y}$, $B_{z}$ and the pressure $P$, equation (1) is twice elliptic and twice hyperbolic [1]. Such mixed behavior makes (1) one of the hardest equations in mathematical physics. Indeed, while in the presence of a continuous Euclidean symmetry the hyperbolic part can be removed to obtain an elliptic nonlinear second order partial differential equation for the flux function (the Grad-Shafranov equation [2, 3, 4]), the existence of regular solutions of (1) without such isometries remains an unsolved problem in the theory of partial differential equations [5].

In the following, we shall always use the word symmetry to refer to continuous transformations of three dimensional Euclidean space that preserve the Euclidean distance between points, i.e. superpositions of translations and rotations. Asymmetry will then imply absence of such Euclidean isometries. Furthermore, a certain equilibrium will be deemed symmetric as long as the magnetic field is, without any requirements on other fields or boundary shape. It should be noted that a symmetric solution $\left({\boldsymbol{B},P}\right)$ of (1) does not necessarily imply a symmetric boundary $\partial\Omega$. For example, a symmetric magnetic field satisfying (1) can be embedded within an asymmetric bounded domain. A symmetric solution always implies a symmetric boundary when the boundary corresponds to an isobaric surface. On these points, see [6]. According to a conjecture due to H. Grad, well-behaved solutions of (1) with non-constant pressure in $\Omega$ and constant pressure on $\partial\Omega$ should possess a high degree of symmetry [7]. Intuitively, this is because the boundary condition for the pressure effectively forces magnetic field $\boldsymbol{B}$ and current $\nabla\times\boldsymbol{B}$ on toroidal surfaces (hairy ball theorem [8]). In the absence of a boundary symmetry, this topological obstruction cannot be trivially overcome when one tries to lie $\boldsymbol{B}$ and $\nabla\times\boldsymbol{B}$ consistently around the torus. A workaround for this problem consists in expanding the class of admissible solutions so that controlled discontinuities are allowed. In this construction, the domain $\Omega$ is partitioned into a given number of subdomains. Within each subdomain a constant pressure is assumed, and the corresponding magnetic field is a Beltrami field (an eigenstate of the curl operator, see [9, 10, 11] for the existence of such solutions), while pressure jumps occur at the boundaries among subdomains. This approach has proven effective for the modeling of three-dimensional magnetohydrodynamic equilibria aimed at stellarator design [12, 13, 14] (stellarators are candidate magnetic confinement devices for nuclear fusion applications that do not exhibit boundary symmetry).

Compared to an axially symmetric tokamak, a stellarator offers the advantage that the field line twist required to trap charged particles is not induced by currents within the plasma, but sustained through asymmetric coils, thus enabling improved steady state operation of the machine. This usually comes at the price of deteriorated confinement, and a nontrivial coil design that should be compatible with a quasisymmetric confining magnetic field [15], i.e. a magnetic field $\boldsymbol{B}$ whose strength $B=\sqrt{\boldsymbol{B}\cdot\boldsymbol{B}}$ is independent of a direction $\boldsymbol{u}$ in space, $\boldsymbol{u}\cdot\nabla B=0$ (see below for a rigorous definition of quasisymmetric vector field). Here, quasisymmetry is required because it guarantees the existence of a conserved momentum arising from the invariance along $\boldsymbol{u}$ of the guiding center Hamiltonian, which depends on the magnetic field only through the strength $B$ [16]. The conserved quantity is expected to enhance particle confinement along the magnetic field.

A quasisymmetric magnetic field is mathematically characterized by the property that there exists a solenoidal vector field $\boldsymbol{u}$ such that both $\boldsymbol{B}$ and the field strength $B$ are Lie-invariant along $\boldsymbol{u}$, i.e.

$$\mathfrak{L}_{\boldsymbol{u}}\boldsymbol{B}=\boldsymbol{0},~{}~{}~{}~{}\mathfrak{L}_{\boldsymbol{u}}B=0,~{}~{}~{}~{}\mathfrak{L}_{\boldsymbol{u}}dV=0,$$

(2)

where $\mathfrak{L}$ denotes the Lie derivative and $dV=dxdydz$ the volume element in $\mathbb{R}^{3}$. Clearly, a symmetric magnetic field is also quasisymmetric. In this case $\boldsymbol{u}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x}$ is the generator of continuous Euclidean isometries, with $\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{3}$ constant vectors and $\boldsymbol{x}$ the position vector in $\mathbb{R}^{3}$. Note that (1a) can also be expressed through the Lie derivative as $\mathfrak{L}_{\nabla\times{\boldsymbol{B}}}\boldsymbol{B}=\boldsymbol{0}$ and $\mathfrak{L}_{\boldsymbol{B}}dV=0$. The quasisymmetry condition (2) is usually written in vector notation as [17],

$$\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta,~{}~{}~{}~{}\boldsymbol{u}\cdot\nabla B=0,~{}~{}~{}~{}\nabla\cdot\boldsymbol{u}=0.$$

(3)

Here, $\zeta$ is any function. Nonetheless, in stellarator applications it is customary to identify $\zeta$ with the flux function $\Psi$ taking constant values on the bounding surface, so that both the magnetic field $\boldsymbol{B}$ and the quasisymmetry $\boldsymbol{u}$ lie on flux surfaces, and the conserved momentum is a combination of the poloidal and toroidal momenta. Unfortunately, the analysis of [18] suggests that fulfilling (1) and (3) simultaneously with a constant pressure at the bounding surface in an asymmetric configuration should not be possible (see also [19, 20, 21, 22]). This is because the Lie-invariance required for a quasisymmetric solution to exist cannot be satisfied through the available geometrical degrees of freedom. A possible strategy to overcome this difficulty is to modify the force balance equation (the first equation in (1a)) by introducing pressure anisotropy. In this regard, it has been suggested that pressure anisotropy could remove the problem of overdetermination encountered in the isotropic setting [17, 23]. In the anisotropic case, anisotropic magnetohydrodynamic equilibria are described by the boundary value problem

$$\begin{split}&\left({\nabla\times\boldsymbol{B}}\right)\times\boldsymbol{B}=\nabla\cdot\Pi,~{}~{}~{}~{}\nabla\cdot\boldsymbol{B}=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega,\\ &\boldsymbol{B}\cdot\boldsymbol{n}=0~{}~{}~{}~{}{\rm on}~{}~{}\partial\Omega,\end{split}$$

(4)

where in Cartesian coordinates $\left({x^{1},x^{2},x^{3}}\right)=\left({x,y,z}\right)$ the twice contravariant symmetric (here, symmetric does not refer to continuous Euclidean isometries but to the property that the transpose tensor equals the tensor itself) pressure tensor $\Pi$ has components

$$\Pi^{ij}=P_{\perp}\delta^{ij}+\frac{\left({P_{\parallel}-P_{\perp}}\right)}{B^{2}}B^{i}B^{j},~{}~{}~{}~{}i,j=1,2,3.$$

(5)

The functions $P_{\perp}$ and $P_{\parallel}$ are pressure fields associated with the perpendicular and parallel directions to the magnetic field $\boldsymbol{B}$. Indeed, setting $b^{i}=B^{i}/B$, $i=1,2,3$, one has

$$\nabla\cdot\Pi=\boldsymbol{b}\times\left({\nabla P_{\perp}\times\boldsymbol{b}}\right)+\left({\boldsymbol{b}\cdot\nabla P_{\parallel}}\right)\boldsymbol{b}+\left({P_{\parallel}-P_{\perp}}\right)\nabla\cdot\left({\boldsymbol{b}\boldsymbol{b}}\right).$$

(6)

Hence, the gradient of $P_{\perp}$ generates the pressure force in the direction perpendicular to $\boldsymbol{B}$, while the gradient of $P_{\parallel}$ generates the pressure force in the direction parallel to $\boldsymbol{B}$. Isotropic magnetohydrodynamic equilibria can be recovered when $P_{\perp}=P_{\parallel}=P$, with $P$ the scalar pressure. We also remark that the form of the pressure tensor (5) is rooted in kinetic theory: the pressure fields $P_{\perp}$ and $P_{\parallel}$ represent the averaged squared deviations of the charged particle velocities in the guiding center approximation across and along the magnetic field with respect to the mean perpendicular and parallel flows [7, 24, 25, 26].

Our aim in this paper is to study the boundary value problem (4) with the quasisymmetry condition (3) for the magnetic field. It will be shown that, by appropriately choosing the values of $P_{\perp}$ and $P_{\parallel}$, it is possible to construct regular quasisymmetric magnetic fields, i.e. asymmetric solutions $\boldsymbol{B}$ of system (3), (4) in a neighborhood $U\subset\Omega$ within an asymmetric toroidal volume $\Omega$ such that $\partial U\cap\partial\Omega\neq\emptyset$. We remark that for the time being we shall not be concerned with the feasibility of the obtained equilibria in the laboratory or their stability, but simply focus on their existence.

The present paper is organized as follows. In section 2, general considerations on the geometrical implications of quasisymmetry are given. In particular, equation (3) is written in terms of a set of Clebsch potentials [27]. This form will be useful to build explicit solutions. In section 3, we discuss the representation of asymmetric tori in terms of special sets of coordinates. In section 4, we derive regular symmetric equilibria that solve system (4) within asymmetric toroidal domains. In section 5, we construct regular quasisymmetric solutions of system (4) in a neighborhood $U\subset\Omega$ of an asymmetric toroidal volume $\Omega$ that satisfy tangential boundary conditions on a portion of the boundary $\partial U\cap\partial\Omega\neq\emptyset$. Concluding remarks are given in section 5.

Consider Cartesian coordinates $\left({x^{1},x^{2},x^{3}}\right)=\left({x,y,z}\right)$ in $\Omega$. It is convenient to write the pressure fields $P_{\perp}$ and $P_{\parallel}$ as follows:

$$P_{\perp}=\mathcal{P}-\frac{1}{2}\sigma B^{2},~{}~{}~{}~{}P_{\parallel}=\mathcal{P}+\frac{1}{2}\sigma B^{2}.$$

(7)

Here, the functions $\mathcal{P}$ and $\sigma$ can be interpreted as a reference pressure field and a deviation (anisotropy) term respectively. The Cartesian components of the pressure tensor become

$$\Pi^{ij}=\left({\mathcal{P}-\frac{1}{2}\sigma B^{2}}\right)\delta^{ij}+\sigma B^{i}B^{j},~{}~{}~{}~{}i,j=1,2,3.$$

(8)

In the following, we denote with $\partial_{i}$, $i=1,2,3$, the tangent vector in the $x^{i}$ direction. Then, using $\nabla\cdot\boldsymbol{B}=0$, one can verify the identity

$$\nabla\cdot\Pi=\left({\partial_{i}\Pi^{ij}}\right)\partial_{j}=\nabla\left({\mathcal{P}-\frac{1}{2}\sigma B^{2}}\right)+\boldsymbol{B}\cdot\nabla\left({\sigma\boldsymbol{B}}\right).$$

(9)

Hence, force balance now takes the form

$$\left({1-\sigma}\right)\left({\nabla\times\boldsymbol{B}}\right)\times\boldsymbol{B}=\nabla\mathcal{P}-\frac{1}{2}B^{2}\nabla\sigma+\left({\boldsymbol{B}\cdot\nabla\sigma}\right)\boldsymbol{B}.$$

(10)

In the presence of symmetry, it is known that system (10) can be reduced to a well-posed Grad-Shafranov type equation under appropriate ellipticity conditions, among which the requirement $1-\sigma>0$ [7, 26, 28]. Isotropic equilibria correspond to the case $\mathcal{P}=P$ and $\sigma=0$. More generally, solutions with constant $\sigma\neq 1$ are qualitatively equivalent to isotropic configurations with $P=\mathcal{P}/(1-\sigma)$. Next, let $P_{0}\in\mathbb{R}$ be a constant. The force balance equation (10) can be trivially satisfied by setting $\sigma=1$ and $\mathcal{P}=P_{0}/2$ so that

$$P_{\perp}=\frac{1}{2}\left({P_{0}-B^{2}}\right),~{}~{}~{}~{}P_{\parallel}=\frac{1}{2}\left({P_{0}+B^{2}}\right),$$

(11)

With this choice, it follows that finding a solution of (4) under the quasisymmetry condition (3) now amounts to solving the boundary value problem

$$\begin{split}&\nabla\cdot\boldsymbol{B}=0,~{}~{}~{}~{}\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta,~{}~{}~{}~{}\boldsymbol{u}\cdot\nabla B=0,~{}~{}~{}~{}\nabla\cdot\boldsymbol{u}=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega,\\ &\boldsymbol{B}\cdot\boldsymbol{n}=0~{}~{}~{}~{}{\rm on}~{}~{}\partial\Omega.\end{split}$$

(12)

Notice that these are just the equations satisfied by a quasisymmetric vector field tangential to $\partial\Omega$. Furthermore, observe that in this construction the pressure fields $P_{\perp}$ and $P_{\parallel}$ are not constrained at the boundary. Instead, they are functions of the magnetic energy density. In particular, the perpendicular pressure field $P_{\perp}$ exactly balances the magnetic pressure since $P_{\perp}+B^{2}/2=P_{0}/2$, the mechanical pressure along the magnetic field $P_{\parallel}$ satisfies $P_{\parallel}+B^{2}/2=P_{0}/2+B^{2}$, while the total energy density associated with pressure and magnetic fields is $P_{\perp}+P_{\parallel}+B^{2}/2=P_{0}+B^{2}/2$. If one further demands that, for example, $P_{\perp}=0$ on $\partial\Omega$, this results in an additional boundary condition for the magnetic field strength, $\boldsymbol{B}^{2}=P_{0}$ on $\partial\Omega$. We will briefly return to this point later on, although the present study will be mainly focused on the case in which the pressure fields are not constrained at the bounding surface. We also remark that, if the quasisymmetry condition (3) is not imposed, any solenoidal magnetic field satisfying tangential boundary conditions represents a solution of anisotropic magnetohydrodynamics (4) as long as the pressure fields are chosen as in (11).

It is useful to briefly discuss the physical consistency of equation (11). First, observe that both $P_{\perp}$ and $P_{\parallel}$ can be made non-negative by taking a sufficiently large constant $P_{0}$. Now suppose that the curvature $\left\lvert{\boldsymbol{b}\cdot\nabla\boldsymbol{b}}\right\rvert$ is a small quantity, and rewrite the force balance equation $\left({\nabla\times\boldsymbol{B}}\right)\times\boldsymbol{B}=\nabla\cdot\Pi$ in the following form:

$$-\frac{1}{2}\boldsymbol{b}\times\left({\nabla B^{2}\times\boldsymbol{b}}\right)+B^{2}\boldsymbol{b}\cdot\nabla\boldsymbol{b}=\boldsymbol{b}\times\left({\nabla P_{\perp}\times\boldsymbol{b}}\right)+\left({P_{\parallel}-P_{\perp}}\right)\boldsymbol{b}\cdot\nabla\boldsymbol{b}+\left[\boldsymbol{b}\cdot\nabla P_{\parallel}+\left({P_{\parallel}-P_{\perp}}\right)\left({\nabla\cdot\boldsymbol{b}}\right)\right]\boldsymbol{b}.$$

(13)

Collecting small terms involving the curvature $\left\lvert{\boldsymbol{b}\cdot\nabla\boldsymbol{b}}\right\rvert$, one thus arrives at the system

	$\displaystyle B^{2}=P_{\parallel}-P_{\perp},$		(14a)
	$\displaystyle\boldsymbol{b}\times\left[\nabla\left({P_{\perp}+\frac{B^{2}}{2}}\right)\times\boldsymbol{b}\right]=\boldsymbol{0},$		(14b)
	$\displaystyle\boldsymbol{b}\cdot\nabla P_{\parallel}+\left({P_{\parallel}-P_{\perp}}\right)\left({\nabla\cdot\boldsymbol{b}}\right)=0.$		(14c)

By eliminating $P_{\parallel}$, these equations reduce to

	$\displaystyle\boldsymbol{b}\times\left[\nabla\left({P_{\perp}+\frac{B^{2}}{2}}\right)\times\boldsymbol{b}\right]=\boldsymbol{0},$		(15a)
	$\displaystyle\boldsymbol{b}\cdot\nabla\left({P_{\perp}+\frac{B^{2}}{2}}\right)=0,$		(15b)

which imply equation (11) as desired. Hence, equation (11) corresponds to a magnetic configuration in which the magnetic field curvature $\left\lvert{\boldsymbol{b}\cdot\nabla\boldsymbol{b}}\right\rvert$ represents a higher-order term in the force balance equation.

Now consider system (12). We must distinguish two cases.

1. 1.

   The vector fields $\boldsymbol{B}$ and $\boldsymbol{u}$ are linearly dependent. Then, $\nabla\zeta=\boldsymbol{0}$, and the magnetic field $\boldsymbol{B}$ is self-quasisymmetric. Indeed, system (12) reduces to the magnetic differential equation

   $$\begin{split}&\nabla\cdot\boldsymbol{B}=0,~{}~{}~{}~{}\boldsymbol{B}\cdot\nabla B=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega,\\ &\boldsymbol{B}\cdot\boldsymbol{n}=0~{}~{}~{}~{}{\rm on}~{}~{}\partial\Omega.\end{split}$$

   (16)

   Note that self-quasisymmetry implies that field strength does not change along field lines. Hence, guiding center motion is not accelerated in the direction of the magnetic field, and the conserved momentum is the momentum parallel to $\boldsymbol{B}$. Due to the regularity of the boundary $\partial\Omega$, the unit outward normal $\boldsymbol{n}$ can be expressed as

   $$\boldsymbol{n}=\frac{\nabla\Psi}{\left\lvert{\nabla\Psi}\right\rvert},$$

   (17)

   where $\Psi$ is a single-valued function corresponding to the usual flux function. For the magnetic field $\boldsymbol{B}$ to possess nested flux surfaces, we demand that

   $$\boldsymbol{B}\cdot\nabla\Psi=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

   (18)

   or, denoting with $\Theta$ a possibly multivalued (angle) variable,

   $$\boldsymbol{B}=\nabla\Psi\times\nabla\Theta~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

   (19)

   Recall that, due to the solenoidal nature of $\boldsymbol{B}$, the variable $\Theta$ can always be chosen as a single-valued function in a sufficiently small region $U\subset\Omega$ (Lie-Darboux theorem [29]). In the following, we shall refer to functions associated with the representation of vector fields such as $\Psi$ and $\Theta$ as Clebsch potentials. It should be noted that Clebsch representations of the form (19) always correspond to helicity-free configurations since the vector potential $\boldsymbol{A}$ associated with the magnetic field $\boldsymbol{B}=\nabla\times\boldsymbol{A}$ can be chosen as an integrable vector field, $\boldsymbol{A}=\Psi\nabla\Theta$. Using (19), system (16) reduces to a single nonlinear first-order partial differential equation for the variable $\Theta$,

   $$\left\lvert{\nabla\Psi\times\nabla\Theta}\right\rvert^{2}=E\left({\Psi,\Theta}\right)~{}~{}~{}~{}{\rm in}~{}~{}\Omega,$$

   (20)

   where $E=E\left({\Psi,\Theta}\right)$ is any non-negative single-valued function of the variables $\Psi$ and $\Theta$. Assume a non-vanishing magnetic field, $E\neq 0$. Then, by defining the function

   $$\Theta^{\prime}=\int\frac{d\Theta}{\sqrt{E\left({\Psi,\Theta}\right)}},$$

   (21)

   equation (20) can be written as

   $$\left\lvert{\nabla\Psi\times\nabla\Theta^{\prime}}\right\rvert^{2}=1~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

   (22)

   Notice that equation (22) resembles the eikonal equation. Indeed, it is equivalent to

   $$\left\lvert{\nabla\Theta^{\prime}}\right\rvert^{2}=\frac{1}{\left\lvert{\nabla\Psi}\right\rvert^{2}}+\frac{\left({\nabla\Psi\cdot\nabla\Theta^{\prime}}\right)^{2}}{\left\lvert{\nabla\Psi}\right\rvert^{2}}~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

   (23)

   In general, the solution of a first-order partial differential equation such as (23) involves the integration of the associated characteristic system of ordinary differential equations. Even if such solution is obtained, it usually has a local nature. Therefore, the existence of a global self-quasisymmetric magnetic field in $\Omega$ is contingent upon the possibility of extending and/or patching local solutions consistently within $\Omega$. We also remark that, if one does not enforce boundary conditions, constructing self-quasisymmetric fields becomes a simpler task. For example, denoting with $\left({r,\phi,z}\right)$ cylindrical coordinates, the vector field

   $$\boldsymbol{B}=f\left({r,\frac{z}{r}-\phi}\right)\nabla r\times\nabla\left({\frac{z}{r}-\phi}\right),$$

   (24)

   with $f$ an arbitrary function of $r$ and $z/r-\phi$, satisfies $\nabla\cdot\boldsymbol{B}=0$, $\boldsymbol{B}\cdot\nabla B=0$, and it does not possess continuous Euclidean symmetries in general because the equation $\mathfrak{L}_{\boldsymbol{u}}\boldsymbol{B}=\boldsymbol{0}$ with $\boldsymbol{u}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x}$ does not have solutions such that $\boldsymbol{u}\neq\boldsymbol{0}$. In particular, the quantity $z/r-\phi$ does not correspond to an Euclidean helical symmetry due to the radial dependence. In addition, the vector field (24) is tangential to cylindrical surfaces (see figure 1). Finally, note that in order to validate the asymmetry of a solution in most cases it is easier to verify the violation of the condition $\mathfrak{L}_{\boldsymbol{u}}B^{2}=\boldsymbol{u}\cdot\nabla B^{2}=0$. Indeed, choosing a curvilinear coordinate system $\left({y^{1},y^{2},y^{3}}\right)$ such that $\boldsymbol{u}=\partial/\partial y^{3}$ and denoting with $B^{k}$ and $g_{k\ell}$, $k,\ell=1,2,3$, the components of magnetic field and the Euclidean metric tensor with respect to the new coordinates, it follows that $\boldsymbol{u}\cdot\nabla B^{2}=2g_{k\ell}B^{k}\partial B^{\ell}/\partial y^{3}$ since by hypothesis $\boldsymbol{u}$ is an Euclidean isometry such that $\partial g_{k\ell}/\partial y^{3}=0$, $k,\ell=1,2,3$. On the other hand, $\mathfrak{L}_{\boldsymbol{u}}\boldsymbol{B}=\left({\partial B^{i}/\partial y^{3}}\right)\left({\partial/\partial y^{i}}\right)$, which implies that $\mathfrak{L}_{\boldsymbol{u}}\boldsymbol{B}$ cannot vanish as long as $\mathfrak{L}_{\boldsymbol{u}}B^{2}$ is different from zero.

   

2. 2.

   The vector fields $\boldsymbol{B}$ and $\boldsymbol{u}$ are linearly independent, implying that $\nabla\zeta\neq\boldsymbol{0}$. From the second equation in (12), one sees that both $\boldsymbol{B}$ and $\boldsymbol{u}$ must be orthogonal to $\nabla\zeta$. However, notice that on the boundary $\nabla\zeta$ is not necessarily aligned with the unit outward normal $\boldsymbol{n}$. Hence, in general, $\zeta$ should not be identified with the flux function $\Psi$ taking constant values on $\partial\Omega$. Denoting with $\vartheta$ and $\rho$ two possibly multivalued (angle) variables, and recalling that $\boldsymbol{B}$ and $\boldsymbol{u}$ are solenoidal vector fields, the space of solutions of system (12) can be restricted to

   $$\boldsymbol{B}=\nabla\zeta\times\nabla\vartheta,~{}~{}~{}~{}\boldsymbol{u}=\nabla\zeta\times\nabla\rho.$$

   (25)

   Since $\boldsymbol{B}\cdot\nabla\Psi=0$ in $\Omega$, from (25) one has

   $$\Psi=\Psi\left({\zeta,\vartheta}\right).$$

   (26)

   For a given flux function $\Psi=\Psi\left({\boldsymbol{x}}\right)$, we assume that equation (26) can be inverted to obtain

   $$\zeta=\zeta\left({\vartheta,\boldsymbol{x}}\right).$$

   (27)

   Then, upon substituting equations (25) and (27) into equation (12), system (12) reduces to the following coupled nonlinear first order partial differential equations for the Clebsch potentials $\vartheta$ and $\rho$

   $$\nabla\zeta\times\nabla\vartheta\cdot\nabla\rho=1,~{}~{}~{}~{}\left\lvert{\nabla\zeta\times\nabla\vartheta}\right\rvert^{2}=F\left({\zeta,\rho}\right)~{}~{}~{}~{}{\rm in}~{}~{}\Omega,$$

   (28)

   where $F$ is any non-negative single-valued function of the Clebsch potentials $\zeta\left({\vartheta,\boldsymbol{x}}\right)$ and $\rho$. Notice that, since equation (28) is a first-order system, the main hurdle toward finding an integral is again represented by the possibility of extending a local solution to the whole $\Omega$. If Clebsch potentials $\vartheta$ and $\rho$ can be determined so that (28) is satisfied, the corresponding vector fields (25) define a quasisymmetric solution of anisotropic magnetohydrodynamics. Finally, observe that system (28) can be further reduced to a single nonlinear second-order partial differential equation in regions $V\subset\Omega$ where $\nabla\zeta\times\nabla B\neq\boldsymbol{0}$. Indeed, the condition $\boldsymbol{u}\cdot\nabla B=0$ then implies that $\rho=\rho\left({\zeta,B^{2}}\right)$ with $B^{2}=\left\lvert{\nabla\zeta\times\nabla\theta}\right\rvert^{2}$. Therefore, the remaining equation $\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta$ can be written as

   $$\nabla\zeta\times\nabla\theta\cdot\nabla\left\lvert{\nabla\zeta\times\nabla\theta}\right\rvert^{2}=\frac{1}{\frac{\partial\rho}{\partial B^{2}}}~{}~{}~{}~{}{\rm in}~{}~{}V.$$

   (29)

   For a given $\rho=\rho\left({\zeta,B^{2}}\right)$, a quasisymmetric solution in $V$ can thus be obtained by solving the equation above for the unknown variable $\theta$. We remark that equation (29) is degenerate because it is invariant under the transformation $\theta\rightarrow\theta+\xi$, with $\xi=\xi\left({\zeta}\right)$ an arbitrary function of $\zeta$. Unfortunately, equation (29) cannot hold in the whole $\Omega$ (that is we never have $V=\Omega$). To see this, consider again system (12) and suppose that $\nabla\zeta\times\nabla B\neq\boldsymbol{0}$ in $\Omega$. Since both $\boldsymbol{B}$ and $\boldsymbol{u}$ lie on surfaces of constant $\zeta$, the conditions $\nabla\cdot\boldsymbol{u}=0$ and $\boldsymbol{u}\cdot\nabla B=0$ imply that the quasisymmetry must satisfy

   $$\boldsymbol{u}=\sigma\left({\zeta,B}\right)\nabla\zeta\times\nabla B,$$

   (30)

   where $\sigma=\sigma\left({\zeta,B}\right)$ is some function of $\zeta$ and $B$. Using the property $\boldsymbol{B}\cdot\nabla\zeta=0$, the equation $\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta$ therefore reduces to

   $$\sigma\boldsymbol{B}\cdot\nabla B=1.$$

   (31)

   Due to the regularity of the magnetic field and its derivatives, we must have $\sigma\neq 0$ for (31) to hold. Using $\nabla\cdot\boldsymbol{B}=0$ and $\boldsymbol{B}\cdot\boldsymbol{n}=0$ on $\partial\Omega$, the equation above implies that

   $$0=\int_{\Omega}\nabla\cdot\left({B\boldsymbol{B}}\right)dV=\int_{\Omega}\frac{dV}{\sigma}.$$

   (32)

   For the quasisymmetry $\boldsymbol{u}$ to be a continuous function, the fraction $1/\sigma$ must be continuous. Hence, equation (32) can be satisfied only if there exists at least one point $\boldsymbol{x}^{\ast}\in\Omega$ such that $1/\sigma\left({\boldsymbol{x}^{\ast}}\right)=0$, implying $\left\lvert{\boldsymbol{u}\left({\boldsymbol{x}^{\ast}}\right)}\right\rvert=\infty$. This contradicts the regularity of $\boldsymbol{u}$. Therefore, a global regular quasisymmetric configuration cannot exist in the case $\nabla\zeta\times\nabla B\neq\boldsymbol{0}$. In other words, if a global regular quasisymmetric configuration exists, there must be regions/points where $\nabla\zeta\times\nabla B=\boldsymbol{0}$. Evidently, the property $\rho=\rho\left({\zeta,B^{2}}\right)$ breaks down at those points where $\nabla\zeta\times\nabla B=\boldsymbol{0}$, and equation (29) does not hold there.

In this section, we are concerned with the representation of asymmetric tori through special sets of coordinates. Such representation will be used in the following sections to derive symmetric and quasisymmetric anisotropic magnetohydrodynamic equilibria in asymmetric tori. Let $\left({r,\phi,z}\right)$ denote cylindrical coordinates. First, consider an axially symmetric torus. Such surface is defined as a level set of the function

$$\Psi_{\rm ax}=\frac{1}{2}\left[\left({r-r_{0}}\right)^{2}+z^{2}\right],$$

(33)

with $r_{0}$ a positive real constant representing the major radius of the torus (the distance from the origin to the center of the torus). The axial symmetry is described by the vector field $\boldsymbol{u}=\partial_{\phi}=r^{2}\nabla\phi$. Indeed, $\mathfrak{L}_{\partial_{\phi}}\Psi_{\rm ax}=0$. The axially symmetric torus (33) can be thought as a special case of a more general family of toroidal surfaces in $\mathbb{R}^{3}$, which correspond to level sets of the function

$$\Psi_{\rm T}=\frac{1}{2}\left[\left({\mu-\mu_{0}}\right)^{2}+\left({z-h}\right)^{2}\right].$$

(34)

Here, $\mu=\mu\left({x,y}\right)$ is a single-valued function in the $\left({x,y}\right)$ plane measuring the distance from the origin in $\mathbb{R}^{2}$ and $\mu_{0}$ a positive real constant corresponding to the major radius of the torus (more generally $\mu_{0}$ can be a single-valued function, although this case is not considered here). The more the level sets of $\mu$ differ from circles, the more the torus $\Psi_{\rm T}$ departs from axial symmetry in the $\left({x,y}\right)$ plane. Similarly, the single-valued function $h=h\left({\boldsymbol{x}}\right)$ expresses the deviation of the toroidal axis from the $\left({x,y}\right)$ plane. Figure 2 shows an axially symmetric torus (33) and three asymmetric tori (34) obtained from different choices of $\mu$ and $h$. The asymmetry of the tori (b), (c), and (d) in figure 1 can be verified by observing that these surfaces are not Lie-invariant with respect to the generator of continuous Euclidean isometries, that is the equation $\mathfrak{L}_{\boldsymbol{u}}\Psi_{\rm T}=0$ with $\boldsymbol{u}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x}$ does not have solution for any nontrivial choice of the constant vectors $\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{3}$.

If quasisymmetry is not imposed, solutions of (4) in asymmetric tori (34) can be easily obtained by enforcing (11) and setting

$$\boldsymbol{B}=\nabla\Psi_{\rm T}\times\nabla\Phi,$$

(35)

with $\Phi$ an arbitrary function whose gradient field is linearly independent of $\nabla\Psi_{\rm T}$. In this case, a sufficient condition for the perpendicular pressure $P_{\perp}=\left({P_{0}-B^{2}}\right)/2$ to be constant on $\partial\Omega$ (which corresponds to a level set of $\Psi_{\rm T}$ by construction) is that $\Phi$ is chosen by requiring that $\nabla\Psi_{\rm T}\times\nabla P_{\perp}=\boldsymbol{0}$, or, substituting the expression for $P_{\perp}$,

$$\left\lvert{\nabla\Psi_{\rm T}\times\nabla\Phi}\right\rvert^{2}=E\left({\Psi_{\rm{}_{T}}}\right)~{}~{}~{}~{}{\rm in}~{}~{}\Omega,$$

(36)

where $E=E\left({\Psi_{\rm T}}\right)$ is a non-negative function of $\Psi_{\rm T}$. Of course, if $P_{\perp}$ is constant on $\partial\Omega$, so is $P_{\parallel}$ due to (11). Regarding the solvability of (36), considerations analogous to those pertaining to equation (20) apply. Furthermore, by comparison with equation (20), it is clear that a magnetic field fulfilling (36) is also self-quasisymmetric since $\boldsymbol{B}\cdot\nabla B=0$ (recall equation (16)). Solutions of (36) can be obtained easily in axially symmetric tori (33). Indeed, first observe that the vector fields

$$\boldsymbol{B}=\lambda\nabla\phi,$$

(37)

with $\lambda=\lambda\left({r,z}\right)$ are self-quasisymmetric because they satisfy (16) when $\Omega$ is an axially symmetric torus (33). Then, equation (36) holds provided that $\lambda^{2}=r^{2}E\left({\Psi_{\rm ax}}\right)$. Of course, the self-quasisymmetry of (37) corresponds to the usual rotational symmetry because $\mathfrak{L}_{\partial_{\phi}}\boldsymbol{B}=\boldsymbol{0}$. However, the construction leading to (37) can be generalized to a certain class of asymmetric tori and, in particular, to derive symmetric solutions in asymmetric toroidal domains. To see this, notice that the orthogonal coordinates $\left({\log{r},\phi,z}\right)$ represent a special case of a larger family of orthogonal harmonic coordinates $\left({x^{1},x^{2},x^{3}}\right)=\left({\mu,\nu,z}\right)$ satisfying the following properties in $\Omega$:

$$\nabla x^{i}\cdot\nabla x^{j}=\delta^{ij},~{}~{}~{}~{}\Delta x^{i}=0,~{}~{}~{}~{}\left\lvert{\nabla x^{1}}\right\rvert=\left\lvert{\nabla x^{2}}\right\rvert,~{}~{}~{}~{}\left\lvert{\nabla x^{3}}\right\rvert=1,~{}~{}~{}~{}i,j=1,2,3.$$

(38)

The coordinates $\mu$ and $\nu$ can be constructed as follows. Let $\Gamma_{1}$ and $\Gamma_{2}$ denote two smooth non-intersecting Jordan (simple and closed) curves in $\mathbb{R}^{2}$. Let $\Sigma_{1}$ and $\Sigma_{2}$ be the areas enclosed by $\Gamma_{1}$ and $\Gamma_{2}$ respectively. Assume that $\Sigma_{2}\subset\Sigma_{1}$ and define $\Gamma=\Gamma_{1}\cup\Gamma_{2}$ and $\Sigma=\Sigma_{1}-\mkern 1.5mu\overline{\mkern-1.5mu\Sigma\mkern-1.5mu}\mkern 1.5mu_{2}$, where the bar denotes the closure of a set. Then, the coordinate $\mu=\mu\left({x,y}\right)$ is obtained as solution of the following Dirichlet boundary value problem for the Laplace equation,

$$\Delta\mu=0~{}~{}~{}~{}{\rm in}~{}~{}\Sigma,~{}~{}~{}~{}\mu=c_{1}~{}~{}~{}~{}{\rm on}~{}~{}\Gamma_{1},~{}~{}~{}~{}\mu=c_{2}~{}~{}~{}~{}{\rm on}~{}~{}\Gamma_{2},~{}~{}~{}~{}c_{1},c_{2}\in\mathbb{R},~{}~{}~{}~{}c_{1}\neq c_{2}.$$

(39)

Observe that $\mu$ can be thought of as a measure of the distance from the origin in $\mathbb{R}^{2}$. Next, the coordinate $\nu=\nu\left({x,y}\right)$ is chosen as the harmonic conjugate of $\mu$. In particular, denoting with $\tilde{\boldsymbol{n}}$ the unit outward normal to the boundary $\Gamma$, the coordinate $\nu$ is a multivalued (angle) variable whose gradient $\nabla\nu$ is a harmonic vector field in $\Sigma$, i.e.

$$\Delta\nu=0~{}~{}~{}~{}{\rm in}~{}~{}\Sigma,~{}~{}~{}~{}\nabla\nu\cdot\tilde{\boldsymbol{n}}=0~{}~{}~{}~{}{\rm on}~{}~{}\Gamma.$$

(40)

Finally, equation (38) is satisfied in virtue of the Cauchy-Riemann equations.

As an example, suppose that $\Gamma_{1}$ and $\Gamma_{2}$ are the boundaries of the horizontal sections of two confocal elliptic cylinders centered at the origin in $\mathbb{R}^{3}$. Then, the orthogonal harmonic coordinates $\left({\mu,\nu,z}\right)$ correspond to the so-called elliptic cylindrical coordinates, which are related to the Cartesian coordinates $\left({x,y,z}\right)$ by the transformation

$$x=a\cosh\mu\cos\nu,~{}~{}~{}~{}y=a\sinh\mu\sin\nu,$$

(41)

with $a$ a positive real constant such that the foci of the elliptic sections of the cylinders are located at $\boldsymbol{x}=\left({a,0,z}\right)$ and $\boldsymbol{x}=\left({-a,0,z}\right)$. Introducing the quantity

$$\delta=\frac{1}{a^{2}\left\lvert{\nabla\mu}\right\rvert^{2}}=\frac{1}{a^{2}\left\lvert{\nabla\nu}\right\rvert^{2}}=\sin^{2}\nu+\sinh^{2}\mu=\sqrt{\left({1-\frac{x^{2}+y^{2}}{a^{2}}}\right)^{2}+\frac{4y^{2}}{a^{2}}},$$

(42)

one can derive the following expressions for the quadrant $x,y\geq 0$,

	$\displaystyle\mu=\operatorname{arcsinh}\left[\frac{1}{\sqrt{2}}\sqrt{-1+\frac{x^{2}+y^{2}}{a^{2}}+\delta}\right],$		(43a)
	$\displaystyle\nu=\arcsin\left[\frac{1}{\sqrt{2}}\sqrt{1-\frac{x^{2}+y^{2}}{a^{2}}+\delta}\right],$		(43b)

and also

	$\displaystyle\nabla\mu=\frac{\sinh\mu\cos\nu\nabla x+\cosh\mu\sin\nu\nabla y}{a\left({\sin^{2}\nu+\sinh^{2}\mu}\right)}=\frac{\frac{\sinh\mu}{\cosh\mu}x\nabla x+\frac{\cosh\mu}{\sinh\mu}y\nabla y}{a^{2}\left({\sin^{2}\nu+\sinh^{2}\mu}\right)},$		(44a)
	$\displaystyle\nabla\nu=\frac{\sinh\mu\cos\nu\nabla y-\cosh\mu\sin\nu\nabla x}{a\left({\sin^{2}\nu+\sinh^{2}\mu}\right)}=\frac{\frac{\sinh\mu}{\cosh\mu}x\nabla y-\frac{\cosh\mu}{\sinh\mu}y\nabla x}{a^{2}\left({\sin^{2}\nu+\sinh^{2}\mu}\right)}.$		(44b)

In the other quadrants, the coordinate $\nu$ can be defined so that $\nu\in[0,2\pi)$ when moving around a $\mu$ contour such as $\Gamma_{1}$ or $\Gamma_{2}$. In particular, $\pi-\nu$ for $x<0$ and $y\geq 0$, $2\pi-\nu$ for $x\geq 0$ and $y<0$, and $\pi+\nu$ for $x,y<0$ with $\nu$ given by (43b). In figure 3, level sets of cylindrical coordinates $\log{r}$, $\phi$ and elliptic cylindrical coordinates $\mu$, $\nu$ as defined in (43) are shown together with the respective gradient fields.

By substituting (43a) into the expression (34), one obtains elliptic toroidal surfaces. In the following sections, we will see how to embed symmetric and quasisymmetric magnetic fields within such domains.

Let $\Omega$ denote the volume enclosed by the elliptic torus

$$\Psi_{\rm el}=\frac{1}{2}\left[\left({\mu-\mu_{0}}\right)^{2}+z^{2}\right],$$

(45)

with $\mu$ given by (43a). Again, the surface $\Psi_{\rm el}$ does not possess continuous Euclidean symmetries because the equation $\mathfrak{L}_{\boldsymbol{u}}\Psi_{\rm el}=0$ with $\boldsymbol{u}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x}$ does not have nontrivial solutions $\boldsymbol{u}\neq\boldsymbol{0}$. In particular, note that horizontal cuts of $\Psi_{\rm el}$ are ellipsoidal. As discussed in the previous section, in the absence of quasisymmetry and boundary conditions on the pressure fields, anisotropic magnetohydrodynamic equilibria within $\Omega$ can be obtained by enforcing (11) and by setting $\boldsymbol{B}=\nabla\Psi_{\rm el}\times\nabla\Phi$ with $\Phi$ an arbitrary function whose gradient field is linearly independent of $\nabla\Psi_{\rm el}$. Observing that $\nabla\Psi_{\rm el}\cdot\nabla\nu=0$, a family of solutions with a more familiar form is

$$\boldsymbol{B}=\nabla\Psi_{\rm el}\times\nabla\nu+\lambda\nabla\nu,$$

(46)

with $\lambda=\lambda\left({\mu,z}\right)$. Next, we seek for a translationally symmetric solution such that the direction of symmetry is $\boldsymbol{u}=\nabla z$. For simplicity, assume that $\boldsymbol{B}=\lambda\nabla\nu$. Evidently, $\nabla\cdot\boldsymbol{B}=\nabla\cdot\boldsymbol{u}=0$. Hence, equation (12) is satisfied provided that

$$\lambda\nabla\nu\times\nabla z=\nabla\zeta,~{}~{}~{}~{}\frac{\partial}{\partial z}\lambda^{2}\left\lvert{\nabla\nu}\right\rvert^{2}=0.$$

(47)

Since by construction $\nabla\nu\times\nabla z=\nabla\mu$ and $\partial\left\lvert{\nabla\nu}\right\rvert/\partial z=0$, these conditions are identically satisfied for any $\lambda=\lambda\left({\mu}\right)$. Then, $\zeta=\int\lambda d\mu$. This example shows that a symmetric solution can be fitted within an asymmetric domain. In figure 4 a plot of a symmetric magnetic field constructed as described above is given.

Let us now consider the existence of self-quasisymmetric magnetic fields (16) within $\Omega$ such that the perpendicular pressure $P_{\perp}$ is constant on the boundary. Such solution can be obtained by solving (36) with $\Psi_{\rm T}=\Psi_{\rm el}$ for the Clebsch potential $\Phi$. Assuming that $B\neq 0$, it is convenient to introduce the quantity $\Phi^{\prime}=\Phi/\sqrt{E\left({\Psi_{\rm el}}\right)}$ so that the equation to be solved becomes

$$\left\lvert{\nabla\Psi_{\rm el}\times\nabla\Phi^{\prime}}\right\rvert^{2}=1~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

(48)

Substituting the expression for $\Psi_{\rm el}$, one arrives at the following nonlinear first-order partial differential equation

$$\left({\frac{\partial\Phi^{\prime}}{\partial\nu}}\right)^{2}\left[z^{2}+\frac{\left({\mu-\mu_{0}}\right)^{2}}{a^{2}\delta}\right]+\left[z\frac{\partial\Phi^{\prime}}{\partial\mu}-\left({\mu-\mu_{0}}\right)\frac{\partial\Phi^{\prime}}{\partial z}\right]^{2}=a^{2}\delta~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

(49)

To further simplify the equation, it is convenient to introduce the change of variables $\left({\mu,z}\right)\rightarrow\left({\rho,\theta}\right)$ defined by

$$\mu-\mu_{0}=\rho\cos\theta,~{}~{}~{}~{}z=\rho\sin\theta.$$

(50)

Note that $\rho^{2}=2\Psi_{\rm el}$. Then, equation (49) reads as

$$\left({\frac{\partial\Phi^{\prime}}{\partial\nu}}\right)^{2}\rho^{2}\left(\sin^{2}\theta+\frac{\cos^{2}\theta}{a^{2}\delta}\right)+\left({\frac{\partial\Phi^{\prime}}{\partial\theta}}\right)^{2}=a^{2}\delta~{}~{}~{}~{}{\rm in}~{}~{}\Omega,$$

(51)

where $\delta=\sin^{2}\nu+\sinh^{2}\left(\mu_{0}+\rho\cos\theta\right)$. For given $\rho$, (51) is a nonlinear first-order partial differential equation with two independent variables quadratic in the derivatives on a flux surface $\Psi_{\rm el}$. Indeed, it has the form

$$F\left({p,q,\nu,\theta}\right)=fp^{2}+q^{2}-g=0,$$

(52)

where

$$p=\frac{\partial\Phi^{\prime}}{\partial\nu},~{}~{}~{}~{}q=\frac{\partial\Phi^{\prime}}{\partial\theta},~{}~{}~{}~{}f=\rho^{2}\left(\sin^{2}\theta+\frac{\cos^{2}\theta}{a^{2}\delta}\right),~{}~{}~{}~{}g=a^{2}\delta.$$

(53)

A local solution of (51) can be found by considering the Cauchy problem for the characteristic system of ordinary differential equations associated with (51) (see e.g. [30]). To this end, denote with $\xi$ a parameter and assign initial conditions as below:

$$\nu=\nu_{0}\left({\xi}\right),~{}~{}~{}~{}\theta=\theta_{0}\left({\xi}\right),~{}~{}~{}~{}\Phi^{\prime}=\Phi^{\prime}_{0}\left({\xi}\right),~{}~{}~{}~{}p=p_{0}\left({\xi}\right),~{}~{}~{}~{}q=q_{0}\left({\xi}\right),$$

(54)

where the functions $\nu_{0}$, $\theta_{0}$, $\Phi^{\prime}_{0}$, $p_{0}$, and $q_{0}$ are determined so that the following non-degeneracy, compatibility, and transversality conditions are fulfilled:

	$\displaystyle\left({\frac{\partial F}{\partial p}}\right)^{2}+\left({\frac{\partial F}{\partial q}}\right)^{2}\neq 0,$		(55a)
	$\displaystyle\left({\frac{d\nu_{0}}{d\xi}}\right)^{2}+\left({\frac{d\theta_{0}}{d\xi}}\right)^{2}\neq 0,$		(55b)
	$\displaystyle F\left({p_{0},q_{0},\nu_{0},\theta_{0}}\right)=0,$		(55c)
	$\displaystyle\frac{d\Phi^{\prime}_{0}}{d\xi}=p_{0}\frac{d\nu_{0}}{d\xi}+q_{0}\frac{d\theta_{0}}{d\xi},$		(55d)
	$\displaystyle\frac{\partial F}{\partial p}\frac{d\theta_{0}}{d\xi}-\frac{\partial F}{\partial q}\frac{d\nu_{0}}{d\xi}\neq 0.$		(55e)

Equation (55a) is identically satisfied as long as $f$, $p$, and/or $q$ are different from zero. Indeed,

$$\left({\frac{\partial F}{\partial p}}\right)^{2}+\left({\frac{\partial F}{\partial q}}\right)^{2}=4f^{2}p^{2}+4q^{2}\neq 0.$$

(56)

Equation (55c) can be solved for $q_{0}$. We have

$$q_{0}^{2}=g_{0}-f_{0}p^{2}_{0}.$$

(57)

Here, $g_{0}=g\left({\nu_{0},\theta_{0}}\right)$ and $f_{0}=f\left({\nu_{0},\theta_{0}}\right)$. The other equations (55b), (55d), and (55e) can be satisfied, for example, by demanding that $\nu_{0}=\xi$, $\theta_{0}=c_{\theta}\in\mathbb{R}$, and $p_{0}=0$. In this case (55b), (55d), and (55e) become

	$\displaystyle\left({\frac{d\nu_{0}}{d\xi}}\right)^{2}+\left({\frac{d\theta_{0}}{d\xi}}\right)^{2}=1\neq 0,$		(58a)
	$\displaystyle\frac{d\Phi^{\prime}_{0}}{d\xi}=p_{0}\frac{d\nu_{0}}{d\xi}+q_{0}\frac{d\theta_{0}}{d\xi}=0,$		(58b)
	$\displaystyle\frac{\partial F}{\partial p}\frac{d\theta_{0}}{d\xi}-\frac{\partial F}{\partial q}\frac{d\nu_{0}}{d\xi}=-2q_{0}\neq 0.$		(58c)

Hence, a set of consistent initial conditions is

$$\nu_{0}=\xi,~{}~{}~{}~{}\theta_{0}=c_{\theta},~{}~{}~{}~{}\Phi^{\prime}_{0}=c_{\Phi^{\prime}},~{}~{}~{}~{}p_{0}=0,~{}~{}~{}~{}q_{0}=\sqrt{g_{0}}.$$

(59)

with $c_{\Phi^{\prime}}\in\mathbb{R}$. Then, solving the characteristic Lagrange-Charpit system

$$\frac{d\nu}{\frac{\partial F}{\partial p}}=\frac{d\theta}{\frac{\partial F}{\partial q}}=\frac{d\Phi^{\prime}}{p\frac{\partial F}{\partial p}+q\frac{\partial F}{\partial q}}=-\frac{dp}{\frac{\partial F}{\partial\nu}}=-\frac{dq}{\frac{\partial F}{\partial\theta}}=d\tau,$$

(60)

one obtains a parametric solution $\nu=\nu\left({\xi,\tau}\right)$, $\theta\left({\xi,\tau}\right)$, and $\Phi^{\prime}\left({\xi,\tau}\right)$ in a neighborhood of the line defined by the initial conditions (59). This result implies that, for a given function $E\left({\Psi_{\rm el}}\right)>0$, a self-quasisymmetric solution $\boldsymbol{B}=\sqrt{E\left({\Psi_{\rm el}}\right)}\nabla\Psi_{\rm el}\times\nabla\Phi^{\prime}$ of anisotropic magnetohydrodynamics with constant perpendicular pressure on the boundary can be constructed in a neighborhood $U\subset\Omega$ such that the boundary condition $\boldsymbol{B}\cdot\boldsymbol{n}=0$ is satisfied on $\partial U\cap\partial\Omega\neq\emptyset$. As shown by the example above, in the context of anisotropic magnetohydrodynamics the existence of quasisymmetric solutions ultimately translates into the global consistency of locally obtained solutions, i.e. whether the solution defined in $U$ can be prolonged to cover the whole $\Omega$. We therefore conjecture that rigorous results concerning the existence of quasisymmetric vector fields cannot be obtained unless the issue of global extension of solutions of first-order nonlinear partial differential equations is carefully addressed. Finally, we observe that the results discussed in the present section apply to asymmetric tori generated trough general orthogonal harmonic coordinates, i.e. they are not restricted to elliptic geometry.

The purpose of the present section is to derive quasisymmetric magnetic fields. We will see that these vector fields can be interpreted as local solutions of anisotropic magnetohydrodynamics (12) in an asymmetric torus with perpendicular and parallel pressure fields given by (11). In particular, these solutions hold in a neighborhood $U\subset\Omega$ with $\partial U\cap\partial\Omega\neq\emptyset$, they do not have continuous Euclidean symmetries, no boundary conditions are imposed on the pressure fields, and $\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta\neq\boldsymbol{0}$ (the fields are not self-quasisymmetric). To achieve this goal, we solve equation (28) in a neighborhood $U$ of a toroidal domain $\Omega$ whose boundary $\partial\Omega$ corresponds to a level set of a function $\Psi_{\rm T}$.

First, we prescribe the toroidal domain $\Omega$. Instead of breaking axial symmetry by modifying the distance function $r$ with a new variable $\mu$, it is convenient to deform the torus by introducing a non-zero $h$ in ${\Psi_{\rm T}}$ of (34). To simplify calculations, we postulate that $h=h\left({r,\phi}\right)$. We thus have

$$\Psi_{\rm T}=\frac{1}{2}\left[\left({r-r_{0}}\right)^{2}+\left({z-h}\right)^{2}\right].$$

(61)

For the magnetic field to possess nested flux surfaces, we further demand that

$$\boldsymbol{B}=f\left({r,z-h}\right)\nabla r\times\nabla\left({z-h}\right).$$

(62)

Here, $f$ is a function of $r$ and $z-h$ to be determined together with $h$ by imposing the quasisymmetry condition. Notice that $\boldsymbol{B}\cdot\nabla\Psi_{\rm T}=0$ implying $\boldsymbol{B}\cdot\boldsymbol{n}=0$ on $\partial\Omega$. The magnetic field also satisfies $\nabla\cdot\boldsymbol{B}=0$. Hence, the remaining equations in (12) are

$$\boldsymbol{B}\times\boldsymbol{u}=\frac{dg}{d\zeta}\nabla\zeta,~{}~{}~{}~{}\boldsymbol{u}\cdot\nabla B=0,~{}~{}~{}~{}\nabla\cdot\boldsymbol{u}=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega.$$

(63)

In the equation above we have replaced $\zeta$ with a function of $\zeta$, $g=g\left({\zeta}\right)$, for later convenience. Assuming the Clebsch representation $\boldsymbol{u}=\nabla\zeta\times\nabla\rho$ with $\zeta=z-h$, from (63) we thus arrive at (28), which now reads as a system of equations that must be solved for the Clebsch potential $\rho$ and the displacement $h$,

	$\displaystyle\nabla r\times\nabla\left({z-h}\right)\cdot\nabla\rho=1,$		(64a)
	$\displaystyle\nabla\left({z-h}\right)\times\nabla\rho\cdot\nabla\left[\frac{1}{r^{2}}\left({\frac{\partial h}{\partial\phi}}\right)^{2}\right]=0.$		(64b)

Observe that in deriving (64a) and (64b) we set $dg/d\zeta=f$ and eliminated the dependence of $f$ on $r$. Next, define the quantity

$$\eta=\frac{1}{r}\frac{\partial h}{\partial\phi}.$$

(65)

Then, assuming that $\nabla\left({z-h}\right)$ and $\nabla\left({r^{-1}\partial h/\partial\phi}\right)^{2}$ are linearly independent, equation (64b) implies that

$$\rho=\rho\left({z-h,\eta}\right).$$

(66)

Substituting (66) into (64a) we arrive at

$$\frac{\partial\rho}{\partial\eta}\nabla r\times\nabla\left({z-h}\right)\cdot\nabla\eta=1,$$

(67)

which, after some manipulations, becomes

$$-\frac{\partial\rho}{\partial\eta}\frac{1}{r^{2}}\frac{\partial^{2}h}{\partial\phi^{2}}=1,$$

(68)

Since $h$ does not depend on $z$, we put $\rho=\rho\left({\eta}\right)$. Then, for given $\rho\left({\eta}\right)$, equation (68) must be solved for $h$. Now recall that $\partial h/\partial\phi=r\eta$. Therefore, equation (67) is equivalent to

$$\frac{\partial\rho}{\partial\phi}=-r.$$

(69)

Assuming that the function $\rho=\rho\left({\eta}\right)$ can be inverted with inverse $\rho^{-1}$ and integrating with respect to $\phi$ gives

$$\frac{\partial h}{\partial\phi}=r\rho^{-1}\left({\alpha-r\phi}\right),$$

(70)

where $\alpha=\alpha\left({r}\right)$ is an arbitrary integration factor. Once $\rho^{-1}$ is assigned, a further integration with respect to $\phi$ gives the desired solution. However, due to the radial dependence of $\alpha-r\phi$, the quantity $\partial h/\partial\phi$ is not periodic in $\phi$. Therefore, the components of the corresponding quasisymmetric magnetic field (62) and the flux function (61) become multivalued functions, and the obtained quasisymmetric solution is well-defined only locally. In particular, the solution is defined in some region $U\subset\Omega$, and the irregularities of the bounding surface $\partial\Omega$ arising from the multivalued nature of $\Psi_{\rm T}$ can be removed by repairing (smoothing) the flux function $\Psi_{\rm T}$ outside $U$. To see this, set $\rho\left({\eta}\right)=k^{-1}\eta$, with $k>0$ a real constant, and choose $f=\sin\left({z-h}\right)$. Then, we have

$$\frac{\partial h}{\partial\phi}=kr\left({\alpha-r\phi}\right),~{}~{}~{}~{}h=kr\phi\left({\alpha-\frac{1}{2}r\phi}\right)+\beta,$$

(71)

with $\beta=\beta\left({r}\right)$ an arbitrary integration factor, and also

	$\displaystyle\boldsymbol{B}=-\sin\left[z-\beta-kr\phi\left({\alpha-\frac{1}{2}r\phi}\right)\right]\left[r\nabla\phi+k\left({\alpha-r\phi}\right)\nabla z\right],$		(72a)
	$\displaystyle\boldsymbol{u}=r\left({\frac{\partial\alpha}{\partial r}-\phi}\right)\nabla\phi+\nabla r+\frac{\partial}{\partial r}\left(k\frac{\alpha^{2}}{2}+\beta\right)\nabla z.$		(72b)

For the component of the magnetic field to be single-valued, one must determine $\alpha$ and $\beta$ so that

$$\boldsymbol{B}\rvert_{\phi=0}=\boldsymbol{B}\rvert_{\phi=2\pi},$$

(73)

which implies

	$\displaystyle z-\beta=z-\beta-2\pi kr\left({\alpha-\pi r}\right)+2\pi m,~{}~{}~{}~{}m\in\mathbb{Z}$		(74a)
	$\displaystyle k\alpha=k\left({\alpha-2\pi r}\right).$		(74b)

However, excluding the case $k=0$ which corresponds to axial symmetry, these conditions cannot be satisfied. Hence, the local nature of the solution.