Existence of Weakly Quasisymmetric Magnetic Fields
in Asymmetric Toroidal Domains
with Non-Tangential Quasisymmetry

Nuclear fusion is a technology with the potential to revolutionize the way energy is harvested. In the approach to nuclear fusion based on magnetic confinement, charged particles (the plasma fuel) are trapped in a doughnut-shaped (toroidal) reactor with the aid of a suitably designed magnetic field. In a classical tokamak [1], the reactor vessel is axially symmetric (see figure 1(a)). The axial symmetry is mathematically described by the independence of physical quantities, such as the magnetic field $\boldsymbol{B}$ and its modulus $B$, from the toroidal angle $\varphi$. Such symmetry is crucial to the quality of tokamak confinement, because it ensures the conservation of the angular momentum $p_{\varphi}$ of charged particles. However, the constancy of $p_{\varphi}$ is not enough to constrain particle orbits in a limited volume because, in addition to the tendency to follow magnetic field lines, particles drift across the magnetic field. This perpendicular drift eventually causes particle loss at the reactor wall, deteriorating the confinement needed to sustain fusion reactions. In a tokamak, perpendicular drifts are therefore suppressed by driving an axial electric current through the confinement region, which generates a poloidal magnetic field in addition to the external magnetic field produced by coils surrounding the confinement vessel (see figures 1(a) and 1(b)). The overall magnetic field therefore forms twisted helical field lines around the torus. Unfortunately, the control of such electric current is difficult because it is maintained by the circulation of the burning fuel itself, making steady operation of the machine a practical challenge.

In contrast to tokamaks, stellarators [2, 3] are designed to confine charged particles through a vacuum magnetic field produced by suitably crafted asymmetric coils (see figure 1(c)). In this context, symmetry is defined as invariance under continuous Euclidean isometries, i.e. transformations of three-dimensional Euclidean space that preserve the Euclidean distance between points. In practice, these transformations are combinations of translations and rotations, with three corresponding types of symmetry: translational, rotational (including axial), and helical. The magnetic field generated by the asymmetric coils of a stellarator is endowed with the field line twist required to minimize particle loss associated with perpendicular drift motion. This removes, in principle, the need to drive an electric current within the confinement region, and thus enables the reactor to operate in a condition close to a steady state (in practice currents may exist in stellarators as well, but they are sensibly smaller than those in a tokamak). Unfortunately, the loss of axial symmetry comes at a heavy price: in general, the angular momentum $p_{\varphi}$ is no longer constant, and confinement is degraded. However, a conserved momentum that spatially constrains particle orbits can be restored if the magnetic field satisfies a more general kind of symmetry, the so-called quasisymmetry [3, 4]. The essential feature of a quasisymmetric magnetic field, whose rigorous definition [5] is given in equation (1), is the invariance $\boldsymbol{u}\cdot\nabla B=0$ of the modulus $B=\left\lvert{\boldsymbol{B}}\right\rvert$ in a certain direction in space $\boldsymbol{u}$ (the quasisymmetry). For completeness, it should be noted that there exist two kinds of quasisymmetry [6, 7, 8, 9]: weak quasisymmetry (the one considered in the present paper), and strong quasisymmetry. In the former, quasisymmetry results in a conserved momentum at first order in the guiding-center expansion, while in the latter the conservation law originates from an exact symmetry of the guiding-center Hamiltonian. Furthermore, the notion of quasisymmetry can be generalized to omnigenity, a property that guarantees the suppression of perpendicular drifts on average [10].

Despite the fact that several stellarators aiming at quasisymmetry or omnigenity have been built [11, 12], that significant efforts are being devoted to stellarator optimization (see e.g. [13]), and that quasisymmetric magnetic fields have been obtained with high numerical accuracy [14], at present the existence of quasisymmetric magnetic fields lacks mathematical proof. This deficiency is rooted in the complexity of the partial differential equations governing quasisymmetry, which are among the hardest in mathematical physics. Indeed, on one hand the toroidal volume where the solution is sought is itself a variable of the problem. On the other hand, the first order nature of the equations prevents general results from being established beyond the existence of local solutions. The availability of quasisymmetric magnetic fields also strongly depends on the additional constraints that are imposed on the magnetic field. For example, if a quasisymmetric magnetic field is sought within the framework of ideal isotropic magnetohydrodynamics, the analysis of [15] suggests that such configurations do not exist (see also [16, 17, 18, 19]) due to an overdetermined system of equations where geometrical constraints outnumber the available degrees of freedom. The issue of overdetermination is less severe [20, 21, 22] if quasisymmetric mgnetic fields correspond to equilibria of ideal anisotropic magnetohydrodynamics [23, 24, 25] where scalar pressure is replaced by a pressure tensor. In this context, it has been shown [26] that local quasisymmetric magnetic fields do exist, although the local nature of the solutions is exemplified by a lack of periodicity around the torus.

The goal of the present paper is to establish the existence of weakly quasisymmetric magnetic fields in toroidal domains by constructing explicit examples. This ‘constructive’ approach has the advantage of bypassing the intrinsic difficulty of the general equations governing quasisymmetry, and hinges upon the method of Clebsch parametrization [27], which provides an effective representation of the involved variables, including the shape of the boundary enclosing the confinement region. The quasisymmetric magnetic fields reported in the present paper hold within asymmetric toroidal volumes, are smooth, have nested flux surfaces, are not invariant under continuous Euclidean isometries, and can be regarded as quilibria of ideal anisotropic magnetohydrodynamics. Nevertheless, these results come with two caveats: since the constructed solutions are optimized only to fulfill weak quasisymmetry, the resulting magnetic fields are not vacuum fields, and their quasisymmetry does not lie on toroidal flux surfaces. Whether these two properties are consistent with weak quasisymmetry therefore remains an open theoretical issue.

Let $\Omega\subset\mathbb{R}^{3}$ denote a smooth bounded domain with boundary $\partial\Omega$. In the context of stellarator design $\Omega$ represents the volume occupied by the magnetically confined plasma, while the bounding surface $\partial\Omega\simeq{\rm T}^{2}$ has the topology of a torus (a 2-dimensional manifold of genus 1). It is important to observe that, in contrast with conventional tokamak design, the vessel $\partial\Omega$ of a stellarator does not exhibit neither axial nor helical symmetry. In $\Omega$, a stationary magnetic field $\boldsymbol{B}\left({\boldsymbol{x}}\right)$ is said to be weakly quasisymmetric provided that there exist a vector field $\boldsymbol{u}\left({\boldsymbol{x}}\right)$ and a function $\zeta\left({\boldsymbol{x}}\right)$ such that the following system of partial differential equations holds,

	$\displaystyle\nabla\cdot\boldsymbol{B}$	$\displaystyle=0,~{}~{}~{}~{}\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta,~{}~{}~{}~{}\nabla\cdot\boldsymbol{u}=0,~{}~{}~{}~{}\boldsymbol{u}\cdot\nabla B^{2}=0~{}~{}~{}~{}{\rm in}~{}~{}\Omega,$		(1a)
	$\displaystyle\boldsymbol{B}\cdot\boldsymbol{n}$	$\displaystyle=0~{}~{}~{}~{}{\rm on}~{}~{}\partial\Omega,$		(1b)

where $B=\left\lvert{\boldsymbol{B}}\right\rvert$ is the modulus of $\boldsymbol{B}$, $\boldsymbol{n}$ denotes the unit outward normal to $\partial\Omega$, and $\boldsymbol{u}$ is the direction of quasisymmetry. As previously explained, system (1a) ensures the existence of a conserved momentum at first order in the guiding center ordering that is expected to improve particle confinement. Usually, the function $\zeta$ is identified with the flux function $\Psi$ so that both $\boldsymbol{B}$ and $\boldsymbol{u}$ lie on flux surfaces $\Psi={\rm constant}$ and the conserved momentum originating from the quasisymmetry is well approximated by the flux function. Although this property is highly desirable from a confinement perspective because it confines particle orbits into a bounded region, if only weak quasisymmetry (1) is sought $\zeta$ and $\Psi$ may differ (see e.g. [5]). In particular, allowing configurations with $\zeta\neq\Psi$ leaves the interesting possibility of achieving good confinement if the level sets of $\zeta$ enclose bounded regions with a topology that may depart from a torus. Mathematically, the four equations in system (1a) represent so-called Lie-symmetries of the solution, i.e. the vanishing of the Lie-derivative $\mathfrak{L}_{\boldsymbol{\xi}}T$ quantifying the infinitesimal difference between the value of a tensor field $T$ at a given point and that obtained by advecting the tensor field along the flow generated by the vector field $\boldsymbol{\xi}$. Specifically, the first equation and the third equation, which imply that both $\boldsymbol{B}$ and $\boldsymbol{u}$ are solenoidal vector fields, express conservation of volumes advected along $\boldsymbol{B}$ and $\boldsymbol{u}$ according to $\mathfrak{L}_{\boldsymbol{B}}dV=\mathfrak{L}_{\boldsymbol{u}}dV=\left({\nabla\cdot\boldsymbol{B}}\right)dV=\left({\nabla\cdot\boldsymbol{u}}\right)dV=0$, where $dV=dxdydz$ is the volume element in $\mathbb{R}^{3}$. Similarly, the second equation in (1a) expresses the invariance of the vector field $\boldsymbol{B}$ along $\boldsymbol{u}$ according to $\mathfrak{L}_{\boldsymbol{u}}\boldsymbol{B}=\boldsymbol{u}\cdot\nabla\boldsymbol{B}-\boldsymbol{B}\cdot\nabla\boldsymbol{u}=\nabla\times\left({\boldsymbol{B}\times\boldsymbol{u}}\right)=\boldsymbol{0}$, while the fourth equation expresses the invariance of the modulus $B^{2}$ along $\boldsymbol{u}$, i.e. $\mathfrak{L}_{\boldsymbol{u}}B^{2}=\boldsymbol{u}\cdot\nabla B^{2}=0$. For further details on these points see [26].

The construction of a solution of (1) requires the simultaneous optimization of $\boldsymbol{B}$, $\boldsymbol{u}$, $\zeta$ and the shape of the boundary $\partial\Omega$. Indeed, assigning the bounding surface $\partial\Omega$ from the outset will generally prevent the existence of solutions due to overdetermination (the available degrees of freedom are not sufficient to satisfy the quasisymmetry equations). A convenient way to simultaneously optimize $\boldsymbol{B}$, $\boldsymbol{u}$, $\zeta$, and $\partial\Omega$ is to use Clebsch parameters [27], which enable the enforcement of the topological requirement on $\partial\Omega$, which must be a torus, and the extraction of the remaining geometrical degrees of freedom for $\boldsymbol{B}$, $\boldsymbol{u}$, and $\zeta$. To see this, first observe that the unit outward normal $\boldsymbol{n}$ to the boundary $\partial\Omega$ can be expressed through the flux function $\Psi$, which is assumed to exist, according to $\boldsymbol{n}=\nabla\Psi/\left\lvert{\nabla\Psi}\right\rvert$. Next, parametrize $\boldsymbol{B}$ and $\boldsymbol{u}$ as

$$\boldsymbol{B}=\nabla\beta_{1}\times\nabla\beta_{2},~{}~{}~{}~{}\boldsymbol{u}=\nabla u_{1}\times\nabla u_{2},$$

(2)

where the Clebsch parameters $\beta_{1}$, $\beta_{2}$, $u_{1}$, and $u_{2}$ are (possibly multivalued) functions that must be determined from the quasisymmetry equations (1). Here, it should be noted that, due to the Lie-Darboux theorem [28], for a given smooth solenoidal vector field $\boldsymbol{v}$ one can always find single valued functions $\alpha_{1}$ and $\alpha_{2}$ defined in a sufficiently small neighborhood $U$ of a chosen point $\boldsymbol{x}\in\Omega$ such that $\boldsymbol{v}=\nabla\alpha_{1}\times\nabla\alpha_{2}$ in $U$. Using the parametrization (2), system (1) reduces to

$$\left({\nabla\beta_{1}\times\nabla\beta_{2}}\right)\times\left({\nabla u_{1}\times\nabla u_{2}}\right)=\nabla\zeta,~{}~{}~{}~{}\left\lvert{\nabla\beta_{1}\times\nabla\beta_{2}}\right\rvert^{2}=f_{B}\left({u_{1},u_{2}}\right),~{}~{}~{}~{}\Psi=\Psi\left({\beta_{1},\beta_{2}}\right).$$

(3)

In going from (1) to (3) we used the fact that the first and third equations in (1a) are identically satisfied. Furthermore, assuming $\boldsymbol{u}\neq\boldsymbol{0}$, the third equation in (1a) implies that the modulus $B^{2}$ must be a function $f_{B}\left({u_{1},u_{2}}\right)$ of $u_{1}$ and $u_{2}$. Similarly, assuming that the magnetic field $\boldsymbol{B}\neq\boldsymbol{0}$ lies on flux surfaces one has $\boldsymbol{B}\cdot\nabla\Psi=0$ in $\Omega$, which implies that $\Psi$ must be a function of $\beta_{1}$ and $\beta_{2}$. The condition $\Psi=\Psi\left({\beta_{1},\beta_{2}}\right)$ also ensures that boundary conditions (1b) are fulfilled because $\boldsymbol{n}=\nabla\Psi/\left\lvert{\nabla\Psi}\right\rvert$.

Now our task is to solve system (3) by determining $\beta_{1}$, $\beta_{2}$, $u_{1}$, $u_{2}$, $f_{B}$, $\zeta$, and $\Psi$ so that the level sets of $\Psi$ define toroidal surfaces. Direct integration of (3) is a mathematically difficult task due to the number and complexity of the geometric constraints involved. Therefore, it is convenient to start from known special solutions corresponding to axially symmetric configurations, and then perform a tailored symmetry breaking generalization. The simplest axially symmetric vacuum magnetic field is given by

$$\boldsymbol{B}_{0}=\nabla\varphi=\nabla z\times\nabla\log r.$$

(4)

The magnetic field (4) satisfies system (1) if, for example, the quasisymmetry is chosen as $\boldsymbol{u}_{0}=\boldsymbol{B}_{0}$. The corresponding flux surfaces are given by axially symmetric tori generated by level sets of the function

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

(5)

with $r_{0}$ a positive real constant representing the radial position of the toroidal axis (major radius). Comparing equation (2) with equations (4) and (5), one sees that $\beta_{1}=u_{1}=z$, $\beta_{2}=u_{2}=\log r$, $B_{0}^{2}=1/r^{2}=e^{-2u_{2}}$, and $\Psi_{0}=\frac{1}{2}\left[\left({e^{\beta_{2}}-r_{0}}\right)^{2}+\beta_{1}^{2}\right]$.

The axially symmetric torus (5) can be generalized to a larger class of toroidal surfaces [26] as

$$\Psi=\frac{1}{2}\left[\left({\mu-\mu_{0}}\right)^{2}+\mathcal{E}\left({z-h}\right)^{2}\right].$$

(6)

In this notation, $\mu$, $\mu_{0}$, $\mathcal{E}$, and $h$ are single valued functions with the following properties. For each $z$, the function $\mu$ measures the distance of a point in the $\left({x,y}\right)$ plane from the origin in $\mathbb{R}^{2}$. The simplest of such measures is the radial coordinate $r$. More generally, on each plane $z={\rm constant}$ level sets of $\mu$ may depart from circles and exhibit, for example, elliptical shape. The function $\mu_{0}$ assigns the $\mu$ value at which the toroidal axis is located. For the axially symmetric torus $\Psi_{0}$, we have $\mu_{0}=r_{0}$. The function $\mathcal{E}>0$ expresses the departure of toroidal cross sections (intersections of the torus with level sets of the toroidal angle) from circles. For example, the axially symmetric torus $\Psi_{\rm ell}=\frac{1}{2}\left[\left({r-r_{0}}\right)^{2}+2z^{2}\right]$ corresponding to $\mathcal{E}=2$ has elliptic cross section. Finally, the function $h$ can be interpreted as a measure of the vertical displacement of the toroidal axis from the $\left({x,y}\right)$ plane. Figure 2 shows different toroidal surfaces generated through (6).

The axial symmetry of the torus $\Psi_{0}$ given by (5) can be broken by introducing dependence on the toroidal angle $\varphi$ in one of the functions $\mu$, $\mu_{0}$, $\mathcal{E}$, or $h$ appearing in (6). Let us set $\mu=r$, take $\mu_{0}$ and $\mathcal{E}$ as positive constants, and consider a symmetry breaking vertical axial displacement $h=h\left({r,\varphi,z}\right)$. For the corresponding $\Psi$ to define a toroidal surface, the function $h$ must be single valued. Hence, $\varphi$ must appear in $h$ as the argument of a periodic function. The simplest ansatz for $h$ is therefore

$$h=\epsilon\sin\left[m\varphi+g\left({r,z}\right)\right].$$

(7)

Here $m\in\mathbb{Z}$ is an integer, $\epsilon$ a positive control parameter such that the standard axially symmetric magnetic field $\boldsymbol{B}_{0}$ with flux surfaces $\Psi_{0}$ can be recovered in the limit $\epsilon\rightarrow 0$, and $g$ a function of $r$ and $z$ to be determined. Now recall that from equation (3) the function $\Psi$ is related to the Clebsch potentials $\beta_{1}$ and $\beta_{2}$ generating the magnetic field $\boldsymbol{B}=\nabla\beta_{1}\times\nabla\beta_{2}$ according to $\Psi\left({\beta_{1},\beta_{2}}\right)$. Comparing with the axially symmetric case (5) we therefore deduce that the analogy holds if $\beta_{1}=z-h$ and $\beta_{2}=\log r$. Defining $\eta=m\varphi+g$, it follows that the candidate quasisymmetric magnetic field is

$$\boldsymbol{B}=\nabla\left(z-\epsilon\sin\eta\right)\times\nabla\log r=\left(1-\epsilon\cos\eta\frac{\partial g}{\partial z}\right)\nabla\varphi+\epsilon m\frac{\cos\eta}{r^{2}}\nabla z,$$

(8)

where $g$ must be determined by enforcing quasisymmetry. Next, observe that

$$B^{2}=\frac{1}{r^{2}}\left[\epsilon^{2}m^{2}\frac{\cos^{2}\eta}{r^{2}}+\left({1-\epsilon\cos\eta\frac{\partial g}{\partial z}}\right)^{2}\right].$$

(9)

An essential feature of quasisymmetry (3) is that the modulus $B^{2}$ can be written as a function of two variables only, $B^{2}=f_{B}\left({u_{1},u_{2}}\right)$. From equation (9) one sees that this result can be achieved by setting $\partial g/\partial z=q\left({r}\right)$ for some radial function $q\left({r}\right)$ so that $u_{1}=\eta$, $u_{2}=\log r$, and also

$$g\left({r,z}\right)=q\left({r}\right)z+v\left({r}\right),$$

(10)

with $v\left({r}\right)$ a radial function. The candidate direction of quasisymmetry is therefore

$$\boldsymbol{u}=\sigma\left({\eta,r}\right)\nabla\eta\times\nabla\log r=\sigma\left({\eta,r}\right)\left({q\nabla\varphi-\frac{m}{r^{2}}\nabla z}\right),$$

(11)

with $\sigma\left({\eta,r}\right)$ a function of $\eta$ and $r$ to be determined. Since by construction $B^{2}=B^{2}\left({u_{1},u_{2}}\right)$, $\Psi=\Psi\left({\beta_{1},\beta_{2}}\right)$, and both $\boldsymbol{B}$ and $\boldsymbol{u}$ as given by (8) and (11) are solenoidal, the only remaining equation in system (3) to be satisfied is the first one. In particular, we have

$$\boldsymbol{B}\times\boldsymbol{u}=\sigma\left({\nabla\varphi-\epsilon\cos\eta\nabla\eta\times\nabla\log r}\right)\times\left({\nabla\eta\times\nabla\log r}\right)=-m\frac{\sigma}{r^{3}}\nabla r.$$

(12)

Hence, upon setting $\sigma=\sigma\left({r}\right)$, system (3) is satisfied with

$$\zeta=-m\int\frac{\sigma}{r^{3}}dr.$$

(13)

Without loss of generality, we may set $\sigma=-r^{3}$ so that $\zeta=mr$ and the quasisymmetric configuration is given by

	$\displaystyle\boldsymbol{B}=$	$\displaystyle\nabla\left[z-\epsilon\sin\left({m\varphi+qz+v}\right)\right]\times\nabla\log r=\left[1-\epsilon\cos\left({m\varphi+qz+v}\right)q\right]\nabla\varphi+\epsilon m\frac{\cos\left({m\varphi+qz+v}\right)}{r^{2}}\nabla z,$		(14a)
	$\displaystyle\boldsymbol{u}=$	$\displaystyle-\frac{1}{3}\nabla\left({m\varphi+qz}\right)\times\nabla r^{3}=mr\nabla z-qr^{3}\nabla\varphi,$		(14b)
	$\displaystyle\Psi=$	$\displaystyle\frac{1}{2}\left\{\left({r-r_{0}}\right)^{2}+\mathcal{E}\left[z-\epsilon\sin\left({m\varphi+qz+v}\right)\right]^{2}\right\},$		(14c)

where $\mathcal{E}$ is a positive real constant.

For the family of solutions (14) to qualify both as quasisymmetric and without continuous Euclidean isometries, we must verify that the magnetic field (14a) is not invariant under some appropriate combination of translations and rotations. To see this, consider the case $q=1/r$ and $v=0$ corresponding to

	$\displaystyle\boldsymbol{B}=$	$\displaystyle\nabla\left[z-\epsilon\sin\left({m\varphi+\frac{z}{r}}\right)\right]\times\nabla\log r=\left[1-\epsilon\frac{\cos\left({m\varphi+\frac{z}{r}}\right)}{r}\right]\nabla\varphi+\epsilon m\frac{\cos\left({m\varphi+\frac{z}{r}}\right)}{r^{2}}\nabla z,$		(15a)
	$\displaystyle\boldsymbol{u}=$	$\displaystyle-\frac{1}{3}\nabla\left({m\varphi+\frac{z}{r}}\right)\times\nabla{r^{3}}=mr\nabla z-r^{2}\nabla\varphi,$		(15b)
	$\displaystyle\Psi=$	$\displaystyle\frac{1}{2}\left\{\left({r-r_{0}}\right)^{2}+\mathcal{E}\left[z-\epsilon\sin\left({m\varphi+\frac{z}{r}}\right)\right]^{2}\right\},$		(15c)

where $\mathcal{E}$ is a positive real constant. Notice that the magnetic field (15a) is smooth in any domain $V\subset\mathbb{R}^{3}$ not containing the vertical axis $r=0$. To exclude the existence of any continuous Euclidean isometry for (15a) it is sufficient to show that the equation

$$\mathfrak{L}_{\boldsymbol{\xi}}B^{2}=\boldsymbol{\xi}\cdot\nabla B^{2}=0,~{}~{}~{}~{}\boldsymbol{\xi}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x},$$

(16)

does not have solution for any choice of constant vector fields $\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{3}$ with $\boldsymbol{a}^{2}+\boldsymbol{b}^{2}\neq 0$. Indeed, since $\boldsymbol{\xi}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x}$ represents the generator of continous Euclidean isometries, the impossibility of satisfying (16) prevents the magnetic field $\boldsymbol{B}$ from possessing translational, axial, or helical symmetry. For further details on this point, see [26]. Next, introducing again $\eta=m\varphi+z/r$, from equation (15a) one has

$$B^{2}=\frac{1}{r^{2}}-2\epsilon\frac{\cos\eta}{r^{3}}+\epsilon^{2}\left({1+m^{2}}\right)\frac{\cos^{2}\eta}{r^{4}}.$$

(17)

It follows that

$$\begin{split}\boldsymbol{\xi}\cdot\nabla B^{2}=&\frac{2}{r^{3}}\left[-1+3\epsilon\frac{\cos\eta}{r}-2\epsilon^{2}\left({1+m^{2}}\right)\frac{\cos^{2}\eta}{r^{2}}\right]\boldsymbol{\xi}\cdot\nabla r\\ &+2\epsilon\frac{\sin\eta}{r^{3}}\left[1-\epsilon\left({1+m^{2}}\right)\frac{\cos\eta}{r}\right]\boldsymbol{\xi}\cdot\nabla\eta.\end{split}$$

(18)

Let $\left({a_{x},a_{y},a_{z}}\right)$ and $\left({b_{x},b_{y},b_{z}}\right)$ denote the Cartesian components of $\boldsymbol{a}$ and $\boldsymbol{b}$. On the surface $\eta=0$, corresponding to $z=z\left({x,y}\right)=-mr\varphi=-m\arctan\left({y/x}\right)\sqrt{x^{2}+y^{2}}$, we have $\sin\eta=0$ and $\cos\eta=1$, and therefore,

$$\boldsymbol{\xi}\cdot\nabla B^{2}=\frac{2}{r^{4}}\left[-1+\frac{3\epsilon}{r}-2\epsilon^{2}\frac{\left({1+m^{2}}\right)}{r^{2}}\right]\left[xa_{x}+ya_{y}+\left({xb_{y}-yb_{x}}\right)z\left({x,y}\right)\right].$$

(19)

This quantity vanishes provided that $a_{x}=a_{y}=b_{x}=b_{y}=0$. Consider now the surface $\eta=\pi/2$, which implies $z=z\left({x,y}\right)=r\left({\pi/2-m\varphi}\right)=\sqrt{x^{2}+y^{2}}\left({\pi/2-m\arctan\left({y/x}\right)}\right)$. In this case $\sin\eta=1$ while $\cos\eta=0$. Furthermore, since the only surviving components in $\boldsymbol{\xi}$ are those coming from $a_{z}$ and $b_{z}$, one has $\boldsymbol{\xi}\cdot\nabla r=0$, and therefore

$$\boldsymbol{\xi}\cdot\nabla B^{2}=\frac{2\epsilon}{r^{3}}\left({\frac{a_{z}}{r}+mb_{z}}\right).$$

(20)

This quantity vanishes provided that $a_{z}=b_{z}=0$. Hence, the quasisymmetric magnetic field (15a) cannot possess continuous Euclidean isometries.

Similarly, the flux function $\Psi$ defined by equation (15c) is not invariant under continuous Euclidean isometries. Indeed, the equation

$$\mathfrak{L}_{\boldsymbol{\xi}}\Psi=\boldsymbol{\xi}\cdot\nabla\Psi=0,~{}~{}~{}~{}\boldsymbol{\xi}=\boldsymbol{a}+\boldsymbol{b}\times\boldsymbol{x},$$

(21)

does not have solution for any nontrivial choice of $\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{3}$. This can be verified easily for $\left\lvert{m}\right\rvert>1$. Indeed, in this case it is sufficient to evaluate $\boldsymbol{\xi}\cdot\nabla\Psi$ over the line $r=r_{0}$, $z=0$ parametrized by $\varphi$. Here, we have

$$\begin{split}\boldsymbol{\xi}\cdot\nabla\Psi=&-\epsilon\mathcal{E}\sin\left({m\varphi}\right)\boldsymbol{\xi}\cdot\nabla\left({z-\epsilon\sin\eta}\right)\\ =&-\epsilon\mathcal{E}\sin\left({m\varphi}\right)\left[a_{z}-\frac{\epsilon a_{z}}{r_{0}}\cos\left({m\varphi}\right)+r_{0}b_{x}\sin\varphi-r_{0}b_{y}\cos\varphi-\epsilon\left({b_{x}-\frac{ma_{x}}{r_{0}}}\right)\sin\varphi\cos\left({m\varphi}\right)\right.\\ &\left.+\epsilon\left({b_{y}-\frac{ma_{y}}{r_{0}}}\right)\cos\varphi\cos\left({m\varphi}\right)-\epsilon mb_{z}\cos\left({m\varphi}\right)\right].\end{split}$$

(22)

This quantity identically vanishes provided that $a_{x}=a_{y}=a_{z}=b_{x}=b_{y}=b_{z}=0$.

Let us examine the properties of the quasisymmetric configuration (15). First, observe that level sets of (15c) define toroidal surfaces (see figure 3(a)), implying that the magnetic field (15a) has nested flux surfaces. Next, note that the function $\zeta$ such that $\boldsymbol{B}\times\boldsymbol{u}=\nabla\zeta$ is proportional to the radial coordinate, i.e. $\zeta=mr$. This function is associated with the conserved momentum $\bar{p}$ generated by the quasisymmetry. In particular, we have [5]

$$\bar{p}=-\frac{1}{\epsilon_{\rm gc}}\zeta+v_{\parallel}\frac{\boldsymbol{u}\cdot\boldsymbol{B}}{B}.$$

(23)

Here, $v_{\parallel}$ denotes the component of the velocity of a charged particle along the magnetic field $\boldsymbol{B}$ while $\epsilon_{\rm gc}\sim\rho/L$ is a small parameter associated with guiding center ordering, $\rho$ the gyroradius, and $L$ a characteristic length scale for the magnetic field. It follows that charged particles moving in the magnetic field (15a) will approximately preserve their radial position since $\bar{p}\approx-\frac{m}{\epsilon_{\rm gc}}r$. This property works in favor of good confinement, although it cannot prevent particles from drifting in the vertical direction. The situation is thus analogous to the case of an axially symmetric vacuum magnetic field $\boldsymbol{B}_{0}=\nabla\varphi$. Level sets of $\zeta=mr$ on a flux surface (15c) are shown in figure 3(b). These contours correspond to magnetic field lines because the magnetic field (15a) is such that $\boldsymbol{B}\cdot\nabla\Psi=\boldsymbol{B}\cdot\nabla r=0$, and field lines are solutions of the ordinary differential equation $\dot{\boldsymbol{x}}=\boldsymbol{B}$. In particular, observe that magnetic field lines are not twisted, and are given by the intersections of the surfaces $\Psi={\rm constant}$ and $r={\rm constant}$, implying that their projection on the $\left({x,y}\right)$ plane is a circle. Plots of the magnetic field (15a) and its modulus $B^{2}$ are given in figures 3(c) and 3(d). It is also worth noticing that the magnetic field (15a) is not a vacuum field. Indeed, it has a non-vanishing current $\boldsymbol{J}=\nabla\times\boldsymbol{B}$ given by

$$\boldsymbol{J}=\frac{\epsilon}{r^{3}}\left[-\left({1+m^{2}}\right)\sin\eta\nabla r+m\left({2r\cos\eta-z\sin\eta}\right)\nabla\varphi+\left({\cos\eta-\frac{z}{r}\sin\eta}\right)\nabla z\right].$$

(24)

Figures 3(e) and 3(f) show plots of the current field $\boldsymbol{J}$ and the corresponding modulus $J^{2}$. The Lorentz force $\boldsymbol{J}\times\boldsymbol{B}$ can be evaluated to be

$$\begin{split}\boldsymbol{J}\times\boldsymbol{B}=&\frac{\epsilon}{r^{4}}\left\{\left[\left({\cos\eta-\frac{z}{r}\sin\eta}\right)\left({\epsilon\left({1+m^{2}}\right)\frac{\cos\eta}{r}-1}\right)+\epsilon m^{2}\frac{\cos^{2}\eta}{r}\right]\nabla r\right.\\ &\left.+\epsilon m\left({1+m^{2}}\right)\sin\eta\cos\eta\nabla\varphi-\left({1+m^{2}}\right)\sin\eta\left({1-\epsilon\frac{\cos\eta}{r}}\right)\nabla z\right\}.\end{split}$$

(25)

It is not difficult to verify that the right-hand side of this equation cannot be written as the gradient of a pressure field $\nabla P$. Hence, the quasisymmetric magnetic field (15a) does not represent an equilibrium of ideal magnetohydrodynamics. Nevertheless, it can be regarded as an equilibrium of anistropic magnetohydrodynamics provided that the components of the pressure tensor are appropriately chosen (on this point, see [26]). Plots of the Lorentz force $\boldsymbol{J}\times\boldsymbol{B}$ and its modulus $\left\lvert{\boldsymbol{J}\times\boldsymbol{B}}\right\rvert^{2}$ are given in figures 3(g) and 3(h). Next, observe that the quasisymmetry $\boldsymbol{u}$ given by equation (15b) is not tangential to flux surfaces $\Psi$. Indeed,

$$\boldsymbol{u}\cdot\nabla\Psi=m\mathcal{E}\left({z-\epsilon\sin\eta}\right)r.$$

(26)

Plots of the quasisymmetry $\boldsymbol{u}$ and its modulus $u^{2}$ can be found in figures 3(i) and 3(j).

Finally, let us consider how the quasisymmetry of the configuration (15) compares with the usual understanding that the modulus of a quasisymmetric magnetic field depends on the flux function $\Psi$ and a linear combination of toroidal angle $\varphi$ and poloidal angle $\vartheta$, i.e. $B^{2}\left({\Psi,M\vartheta-N\varphi}\right)$ with $M,N$ integers. When $B^{2}=B^{2}\left({\Psi,M\vartheta-N\varphi}\right)$, on each flux surface the contours of the modulus $B^{2}$ in the $\left({\varphi,\theta}\right)$ plane form straight lines. For the quasisymmetric magnetic field (15a) we have $B^{2}=B^{2}\left({r,m\varphi+z/r}\right)$. Hence, the correspondence with the usual setting can be obtained by the identification $\Psi\rightarrow r$, $\varphi\rightarrow\varphi$, and $\vartheta\rightarrow z/r$. Figure 4 shows how the contours of the quasisymmetric magnetic field (15a) form straight lines in the $\left({m\varphi,z/r}\right)$ plane. Next, it is useful to determine how much the contours of $B^{2}$ depart from straight lines on each flux surface. To this end, observe that equation (15c) can be inverted to obtain $r\left({\Psi,z/r,\eta}\right)$ with $\eta=m\varphi+z/r$ so that the modulus (17) can be written in the form $B^{2}=B^{2}\left({r\left({\Psi,z/r,\eta}\right),\eta}\right)$. Figure 5 shows contours of $B^{2}$ on the plane $\left({m\varphi,z/r}\right)$ for a fixed value of $\Psi$ and different choices of the parameter $\epsilon$ controlling the degree of asymmetry of the solution. In particular, notice how the solution (15) approaches axial symmetry for smaller values of $\epsilon$.

In conclusion, we have demonstrated the existence of weakly quasisymmetric magnetic fields in toroidal volumes by constructing explicit examples (14) through the method of Clebsch parametrization. The obtained configurations are solutions of system (1) with the following properties. In the optimized toroidal domain $\Omega$, the magnetic field $\boldsymbol{B}$ is smooth and equipped with nested flux surfaces $\Psi$. Both $\boldsymbol{B}$ and $\Psi$ do not exhibit continuous Euclidean isometries, i.e. invariance under an appropriate combination of translations and rotations. The quasisymmetry $\boldsymbol{u}$ is not tangential to toroidal flux surfaces $\Psi$, but lies on surfaces of constant radius $r$. In particular, $\boldsymbol{B}\times\boldsymbol{u}=m\nabla r$ with $m$ an integer while $B^{2}=B^{2}\left({r,m\varphi+z/r}\right)$ in the example (15). The conserved momentum arising from the quasisymmetry is given by (23), which is approximately the radial position of a charged particle. The magnetic field $\boldsymbol{B}$ is not a vacuum field since a current $\boldsymbol{J}=\nabla\times\boldsymbol{B}\neq\boldsymbol{0}$ is present. The obtained quasisymmetric magnetic fields (14a) can be regarded as solutions of anisotropic magnetohydrodynamics if the component of the pressure tensor are appropriately chosen [26].

In addition to providing mathematical proof of existence of solutions to system (1) with the properties described above, this work offers an alternative theoretical framework for the numerical and experimental efforts devoted to modern stellarator design, and possibly paves the way to the development of semi-analytical schemes aimed at the optimization of confining magnetic fields. The next goal of the present theory would be to further improve the obtained results by ascertaining the existence of vacuum solutions $\nabla\times\boldsymbol{B}=\boldsymbol{0}$ of system (1) such that the modulus of the magnetic field can be written as a function of the flux function and a linear combination of toroidal and poloidal angles, $B^{2}=B^{2}\left({\Psi,M\vartheta-N\varphi}\right)$, and in particular to establish the existence of vacuum quasisymmetric configurations with the field line twist required to effectively trap charged particles.