# Shape optimization of harmonic helicity in toroidal domains

Rémi Robin *¹ and Robin Roussel †²

¹Laboratoire de Physique de l'École Normale Supérieure, Mines Paris, Inria, CNRS, ENS-PSL, Sorbonne Université, PSL Research University, Paris, France

²Laboratoire Jacques-Louis Lions, Sorbonne Université, Inria, Paris

October 7, 2024

## Abstract

In this paper, we introduce a new shape functional defined for toroidal domains that we call harmonic helicity, and study its shape optimization. Given a toroidal domain, we consider its associated harmonic field. The latter is the magnetic field obtained uniquely up to normalization when imposing zero normal trace and zero electrical current inside the domain. We then study the helicity of this field, which is a quantity of interest in magneto-hydrodynamics corresponding to the $L^2$ product of the field with its image by the Biot-Savart operator. To do so, we begin by discussing the appropriate functional framework and an equivalent PDE characterization. We then focus on shape optimization, and we identify the shape gradient of the harmonic helicity. Finally, we study and implement an efficient numerical scheme to compute harmonic helicity and its shape gradient using finite elements exterior calculus.

## Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 Prerequisites</td><td>3</td></tr><tr><td>2.1 Functional spaces</td><td>3</td></tr><tr><td>2.2 Harmonic fields</td><td>4</td></tr><tr><td>2.3 Vector potentials and Bevir-Gray formula</td><td>6</td></tr><tr><td>3 Harmonic helicity and its shape derivative</td><td>8</td></tr><tr><td>3.1 Pullback of the De Rham complex on a fixed domain</td><td>9</td></tr><tr><td>3.2 Differentiability of the PDEs</td><td>13</td></tr><tr><td>3.3 Proof of Theorem 1</td><td>15</td></tr><tr><td>4 Approximation by finite element exterior calculus</td><td>15</td></tr><tr><td>4.1 Classical finite elements exterior calculus families</td><td>16</td></tr><tr><td>4.2 Discretization of the De Rham complex</td><td>17</td></tr><tr><td>4.3 Numerical convergence of the harmonic helicity</td><td>18</td></tr><tr><td>5 Numerical implementation and results</td><td>22</td></tr><tr><td>5.1 Specificity of simulations for stellarators</td><td>22</td></tr><tr><td>5.2 Implementation</td><td>22</td></tr><tr><td>5.3 Numerical tests</td><td>22</td></tr><tr><td>5.4 Two optimization programs</td><td>24</td></tr><tr><td>6 Conclusion and perspectives</td><td>25</td></tr><tr><td>A Translation from differential forms and Hodge decomposition</td><td>27</td></tr></table>

<footnote>*remi.robin@inria.fr</footnote>

<footnote>†robin.roussel@sorbonne-universite.fr</footnote>
# 1 Introduction

For a given vector field $F$ on a three-dimensional domain $\Omega$, we define its helicity (also known as Biot–Savart helicity) by the formula

$$H(F) = \frac{1}{4\pi} \int_{\Omega \times \Omega} F(y) \cdot \left( F(x) \times \frac{y-x}{|y-x|^3} \right) dxdy. \quad (1)$$

This quantity plays an important role in plasma physics, fluid dynamics and magnetohydrodynamics (see e.g. [Arn66; AK21]). In the context of electromagnetism, helicity of a magnetic field (called magnetic helicity) can be seen as a scalar quantifying the linkage and twist of the magnetic field [Mof69; Arn14]. Equation (1) can be interpreted as a volumic version of the writhe of a curve.

When working with magnetic fields, that is divergence free vector fields, tangent to the boundary of $\Omega$, a natural connection with vector potentials appears. First note that if one introduces the Biot–Savart operator of $F$

$$BS(F)(y) = \frac{1}{4\pi} \int_{\Omega} \frac{F(x) \times (y-x)}{|y-x|^3} dx, \quad (2)$$

one obtains that the helicity of $F$ is the $L^2$ inner product of $F$ with $BS(F)$. For simply connected domains, the latter remains true if we replace $BS(F)$ by any vector potential of $F$ (we recall that $\mathrm{curl}BS(F) = F$). For three-dimensional domains that are not simply connected, the connection between helicity and vector potentials is slightly more involved and have been established by Bevir and Gray in [BG80] for toroidal domains<fnref n="1" /> and [MV19] for more general ones. We recall the definition and main properties of the magnetic helicity in Section 2.3.

A classical mathematical problem related to helicity studied in e.g. [Can+00; Can+99; Val19; Mon23] is the maximization of the helicity on $H_0(\mathrm{div}^0, \Omega)^2$ with fixed $L^2$ norm. The critical vector field of this optimization problem are in fact eigenfields of the curl operator which can be seen thanks to a modified Biot–Savart operator denoted here $BS'$ [Can+00].

Elements of shape optimization were also discussed in [Can+00] and [Can+99] to characterize the domain with the highest eigenvalue of $BS'$ for a given volume. Some properties of the maximizing fields on such a domain were given, but it is also unclear whether the optimal shape exists, as some computations suggest that it would have to be a singular sphere, with North and South Pole collapsed to a single point. Recent works on the existence of an optimal shape and its characterization can be found in [EP23; Ger23b; Ger23a].

In this paper, we are interested in a slightly different problem. Given a toroidal domain $\Omega$, we consider the set of harmonic fields, that is the set of vector fields that are divergence free, curl free and tangent to the boundary. By classical results of Hodge theory, this set is a one-dimensional vector space. We are then interested in the helicity of a normalized harmonic field of $\Omega$, where the normalization is related to total flux of currents through the central hole of the torus. Thus, for any regular enough toroidal domain, we define a scalar quantity that we call the **harmonic helicity** of the domain $\Omega$.

Designing a numerical scheme to compute this shape functional is not obvious. A close problem is the spectral approximation of the curl operator in multiply connected domains; this has been tackled in [LRV15; Alo+18] using finite elements methods. For efficiency considerations, it is important to avoid the computation of the double integral in Eq. (1) and use another vector potential than the Biot–Savart. Using classical results on vector potentials characterizations [Amr+98] and tools from finite elements exterior calculus [AFW10], we provide efficient numerical approximation schemes and implementation for the harmonic helicity.

Physical motivations for considering harmonic helicity arise from the design of stellarators, advanced nuclear fusion devices that rely on the confinement of intensely hot plasma through a sophisticated magnetic field. A significant challenge arises due to the inherent impossibility of creating a non-zero magnetic field with constant magnitude on an axisymmetric toroidal domain. Additionally, variations in the magnetic field amplitude, typically inversely proportional to the major radius, result in a vertical drift, which can be mitigated through the implementation of a twisted magnetic field. Optimization of the shape of the coils is a very active field [Pau+18; PRS22]. A measure of the twisted nature of the magnetic field is expressed by the magnetic helicity. We refer to [IPW20] for a very nice introduction to the topic.

In contrast to Tokamaks, which are axisymmetric devices inducing a current inside the plasma to generate a twisting magnetic field, stellarators aim for stability without requiring a current within the plasma. Consequently, the magnetic field employed to stabilize the plasma in a stellarator can be reasonably approximated

<footnote n="1">¹ That is bounded open subsets of $\mathbb{R}^3$ homeomorphic to a full torus.</footnote>

<footnote n="2">² we refer to Section 2.1 for the definition of this Hilbert space.</footnote>
as a harmonic field within the domain representing the plasma. Hence, we believe that the optimization of the shape of the plasma to increase the harmonic helicity could give rise to interesting new forms of plasma.

This paper is organized as follows:

* In Section 2, we recall classical notions used throughout the paper. We begin by the definitions of some functional spaces. Then we properly define the harmonic fields, give two equivalent constructions and define their circulations. Next, we give two PDE formulations to characterize vector potentials. Finally, we make the connection between vector potentials and the helicity through the Bevir-Gray formula.

* Section 3 contains the main contribution of the paper. We begin by introducing the precise definition of the harmonic helicity of a toroidal domain. The rest of the section is dedicated to the computation of the shape derivative of the harmonic helicity. This is done using shape differentiation of several PDE problems and a subtle use of Piola transforms. To the best of the authors' knowledge, the employed methodology is original and may hold applicability in addressing new shape differentiation problems involving other Hilbert complexes.

* In Section 4, we recall the framework of finite element exterior calculus. We then use classical results on approximations of Hodge Laplacian problems. The adaptations of these tools to our problem is not straightforward and provides a method to compute both the helicity and the shape gradient.

* In Section 5, we provide numerical results of the proposed numerical methods on specific shapes motivated by the study of stellarator plasmas. Then we present two numerical experiments to improve the harmonic helicity of a standard plasma shape.

* In Appendix A, we recall classical results of Hodge theory used throughout the paper. In particular the Hodge decomposition and its connection with the De Rham cohomology.

# 2 Prerequisites

In this section, $\Omega$ denotes a Lipschitz toroidal domain of $\mathbb{R}^3$, that is $\partial\Omega$ is locally the graph of a Lipschitz function and $\bar{\Omega}$ is homeomorphic to $D^2 \times S^1$ with $D^2$ the closed unit disk and $S^1$ the unit circle.

## 2.1 Functional spaces

We recall the definitions of the following classical functional spaces:

$$ 
\begin{aligned}
H(\mathrm{curl}, \Omega) &= \{V \in L^2(\Omega)^3 \mid \mathrm{curl} V \in L^2(\Omega)^3\}, \\
H(\mathrm{div}, \Omega) &= \{V \in L^2(\Omega)^3 \mid \mathrm{div} V \in L^2(\Omega)\}.
\end{aligned} 
\quad (3) 
$$

On $H(\mathrm{curl}, \Omega)$ and $H(\mathrm{div}, \Omega)$, the tangential and normal traces $V \times n : \partial\Omega \to \mathbb{R}^3$ and $V \cdot n : \partial\Omega \to \mathbb{R}$ are defined respectively by

$$ \int_{\partial\Omega} (V \times n) \cdot \varphi = \int_{\Omega} V \cdot \mathrm{curl} \varphi - \int_{\Omega} \mathrm{curl} V \cdot \varphi, \quad (4) $$

for every $\varphi$ in $H^1(\Omega)^3$, and

$$ \int_{\partial\Omega} (V \cdot n)\varphi = \int_{\Omega} V \cdot \nabla\varphi + \int_{\Omega} \mathrm{div} V\varphi, \quad (5) $$

for every $\varphi$ in $H^1(\Omega)$. Since the traces of $H^1(\Omega)^3$ and $H^1(\Omega)$ are $H^{1/2}(\partial\Omega)^3$ and $H^{1/2}(\partial\Omega)$ respectively, $V \times n$ can be defined in $H^{-1/2}(\partial\Omega)^{33}$, and $V \cdot n$ in $H^{-1/2}(\partial\Omega)$. Then, we can define

$$ 
\begin{aligned}
H_0(\mathrm{curl}, \Omega) &= \{V \in L^2(\Omega)^3 \mid \mathrm{curl} V \in L^2(\Omega)^3, V \times n = 0\}, \\
H_0(\mathrm{div}, \Omega) &= \{V \in L^2(\Omega)^3 \mid \mathrm{div} V \in L^2(\Omega), V \cdot n = 0\}.
\end{aligned} 
\quad (6) 
$$

<footnote n="3">³In fact, as was shown in [BCS02], the space of tangential traces of $H(\mathrm{curl}, \Omega)$ is $H^{-1/2}(\mathrm{div}_{\Gamma}, \partial\Omega)$, so that $H^{-1/2}(\partial\Omega)^3$ is a proper subspace.</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of the curves $\gamma$ and $\gamma'$, and the surfaces $\Sigma$ and $\Sigma'$</figcaption></figure>

We also denote by $L_0^2(\Omega)$ the set of functions in $L^2(\Omega)$ which have zero average in $\Omega$. Introducing the following spaces where the differential operator vanishes will also prove to be useful

$$ 
\begin{aligned}
H(\mathrm{curl}^0, \Omega) &= \{V \in H(\mathrm{curl}, \Omega) | \mathrm{curl} V = 0\}, \\
H(\mathrm{div}^0, \Omega) &= \{V \in H(\mathrm{div}, \Omega) | \mathrm{div} V = 0\}, \\
H_0(\mathrm{curl}^0, \Omega) &= H(\mathrm{curl}^0, \Omega) \cap H_0(\mathrm{curl}, \Omega), \\
H_0(\mathrm{div}^0, \Omega) &= H(\mathrm{div}^0, \Omega) \cap H_0(\mathrm{div}, \Omega).
\end{aligned} 
\quad (7) $$

## 2.2 Harmonic fields

The set of harmonic fields in $\Omega$ is defined in the following way

$$ 
\begin{aligned}
\mathcal{K}(\Omega) &= \{V \in L^2(\Omega)^3 | \mathrm{div} V = 0, \mathrm{curl} V = 0, V \cdot n = 0\} \\
&= H(\mathrm{curl}^0, \Omega) \cap H_0(\mathrm{div}^0, \Omega),
\end{aligned} 
\quad (8) $$

Through the identifications with differential forms given in Appendix A, it is classical that this space is isomorphic to the first De Rham cohomology space of $\Omega$, with each harmonic field giving a natural representent of the corresponding cohomology class. Another characterization based on Hodge decomposition is the following

$$ \mathcal{K}(\Omega) = H(\mathrm{curl}^0, \Omega) \cap \nabla H^1(\Omega)^{\perp}. $$

We refer to Proposition 15 of Appendix A.

We introduce a poloidal cut $\Sigma$ of $\Omega$, that is, a Lipschitz surface included in $\Omega$, with boundary $\gamma$ contained in $\partial\Omega$ which generates the first homology group of $\Omega^c = \mathbb{R}^3 \setminus \Omega$. Similarly, we introduce $\Sigma'$ a Lipschitz surface in $\Omega^c$ with boundary $\gamma'$ that generates the first homology group of $\bar{\Omega}$. Fig. 1 illustrates these objects. We also define $t$ (resp. $t'$) as unit tangent vector fields on $\gamma$ (resp. $\gamma'$), and $n_{\Sigma}$ (resp. $n_{\Sigma'}$) as a unit normal vector field on $\Sigma$ (resp. $\Sigma'$). These vector fields define orientations on $\Sigma$, $\Sigma'$, $\gamma$ and $\gamma'$ which are compatible with each other. We refer to [Alo+18, Section 1] and [Amr+98, Section 3.a] for the constructions and more precise definitions of these objects.

Formally, we normalize the harmonic field $B$ by fixing its circulation along the toroidal loop $\gamma'$ to be equal to $2\pi$. Hence, as shown in [Amr+98, Proposition 3.14], we can define the normalized harmonic field $B$ on $\Omega$ by $B = \nabla u$ with $u \in H^1(\Omega \setminus \Sigma)$ defined as follows. $u$ is the unique solution of the variational problem

$$ 
\left\{ 
\begin{array}{l}
\int_{\Omega \setminus \Sigma} \nabla u \cdot \nabla v = 0, \\
[\![u]\!]_{\Sigma} = 2\pi,
\end{array} 
\right. 
\quad (9) $$

for all $v \in H^1(\Omega)$, where $[\![u]\!]_{\Sigma}$ denote the jump of $u$ across $\Sigma$, that is the difference of traces of $u$ in the direction of $\gamma'$ [Amr+98, Notation 3.9.i]. Note that $\Omega \setminus \Sigma$ is a pseudo-Lipschitz domain, we refer to [Amr+98, Definition 3.1] for further details.
Once we have the normalized harmonic field of $\Omega$, we can define rigorously the circulation of vector fields in $H(\text{curl}, \Omega)$ with curl tangent to the boundary. To do this, we define similarly the harmonic field of $\Omega^{\text{ext}} = \mathcal{B} \cap \bar{\Omega}^c$, denoted $B^{\text{ext}}$, where $\mathcal{B}$ is an open ball containing $\bar{\Omega}$. The normalized harmonic field of $\Omega^{\text{ext}}$ is the gradient of $u^{\text{ext}} \in H^1(\Omega^{\text{ext}}\setminus\Sigma')$, which verifies

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega \setminus \Sigma} \nabla u^{\text{ext}} \cdot \nabla v = 0, \\ \displaystyle \left\| u \right\|_{\Sigma'} = 2\pi, \end{array} \right. \qquad (10) $$

for all $v$ in $H^1(\Omega^{\text{ext}})$. Then, as in [Alo+18, Section 2], we are able to define the circulations of a vector field $V$ in $H(\text{curl}, \Omega)$ with curl tangent to the boundary as

$$ \int_{\gamma'} V \cdot t = \frac{1}{2\pi} \int_{\Omega^{\text{ext}}} \text{curl} \tilde{V} \cdot B^{\text{ext}}, \qquad (11) $$

$$ \int_{\gamma} V \cdot t' = \frac{1}{2\pi} \int_{\Omega} \text{curl} V \cdot B, \qquad (12) $$

where $\tilde{V}$ is a continuous extension of $V$ from $H(\text{curl}, \Omega)$ to $H(\text{curl}, \mathcal{B})$. Similarly, we can also write

$$ \int_{\gamma'} V \cdot t = -\frac{1}{2\pi} \int_{\partial\Omega} (V \times n) \cdot B^{\text{ext}}, \qquad (13) $$

$$ \int_{\gamma} V \cdot t' = -\frac{1}{2\pi} \int_{\partial\Omega} (V \times n) \cdot B. \qquad (14) $$

**Remark 1.** Although this characterization of the normalized harmonic field makes use of $\Sigma$, choosing another cutting surface may change $B$ only up to a sign. Indeed, $B$ is in $\mathcal{K}(\Omega)$ which is one-dimensional. Therefore, changing $\Sigma$ may change $B$ only by multiplying it by a constant. Furthermore, since $B$ integrates to $2\pi$ along $\gamma'$, we see that this constant is equal to 1 if we preserve the orientation of $\Sigma$ and $\gamma'$, and equal to -1 otherwise.

Note that the definition of the normalized harmonic field in Eq. (9) is equivalent to the following mixed formulation that will prove to be useful, both for the computation of the shape gradient, and for the numerical scheme.

**Proposition 1.** here exists a unique solution to the following problem. Find $(B_{\text{div}}, u_{\text{div}}) \in H_0(\text{div}, \Omega) \times L_0^2(\Omega)$ such that, for all $(\tau, v) \in H_0(\text{div}, \Omega) \times L^2(\Omega)$ we have

$$ \begin{array}{l} \displaystyle \int_{\Omega} (\text{div } B_{\text{div}}) v = 0, \\ \displaystyle \int_{\Omega} B_{\text{div}} \cdot \tau + \int_{\Omega} u_{\text{div}} (\text{div } \tau) = 2\pi \int_{\Sigma} \tau \cdot n_{\Sigma}. \end{array} \qquad (15) $$

Furthermore, we have $B_{\text{div}} = B$ defined in Eq. (9).

**Proof.** For the well-posedness, we first notice that this problem is equivalent to the following one. Find $(B_{\text{div}}, u_{\text{div}}, p) \in H_0(\text{div}, \Omega) \times L^2(\Omega) \times \mathbb{R}$ such that

$$ \begin{array}{l} \displaystyle \int_{\Omega} (\text{div } B_{\text{div}}) v + p \int_{\Omega} v = 0, \\ \displaystyle \int_{\Omega} B_{\text{div}} \cdot \tau + \int_{\Omega} u_{\text{div}} (\text{div } \tau) = 2\pi \int_{\Sigma} \tau \cdot n_{\Sigma}, \\ \displaystyle q \int_{\Omega} u_{\text{div}} = 0, \end{array} $$

for all $(\tau, v, q) \in H_0(\text{div}, \Omega) \times L^2(\Omega) \times \mathbb{R}$. Indeed, we find $p = 0$ by choosing $v = 1$ in the first equation. This problem is then equivalent to the mixed Hodge Laplacian studied in [AFW06, Section 7] with essential boundary conditions and $k = 3$. The corresponding bilinear form is then known to satisfy inf-sup conditions [AFW06, Section 7.5, Remark], and we only need to verify that the right-hand side is continuous with respect to $\tau$ to prove well-posedness. This comes from the continuity of the normal trace from $H(\text{div}, \Omega)$ to $H^{-1/2}(\Sigma)$. As a consequence, we get

$$ \begin{aligned} \left| \int_{\Sigma} \tau \cdot n_{\Sigma} \right| &\le \|1\|_{H^{1/2}(\Sigma)} \|\tau \cdot n_{\Sigma}\|_{H^{-1/2}(\Sigma)} \\ &\le C \|\tau\|_{H(\text{div}, \Omega)}. \end{aligned} $$
To prove that $B_{\text{div}}$ and $B$ coincide, we show that $B_{\text{div}}$ is a harmonic field, and that they are normalized in the same way. Equality then follows from the fact that the set of harmonic fields is one-dimensional. To find that $B_{\text{div}}$ is a harmonic field, we prove it is orthogonal to $\nabla H^1(\Omega)$ and curl $H_0(\text{curl}, \Omega)$.

The first fact follows from the first equation of (15). Indeed, we find that $\text{div} B_{\text{div}} = 0$, so that $B_{\text{div}} \in H_0(\text{div}^0, \Omega) = \nabla H^1(\Omega)^\perp$.

Now, take $\tau = \text{curl} \sigma$ in the second equation of (15), with $\sigma \in H_0(\text{curl}, \Omega)$. As a consequence, we have $\text{div} \tau = 0$, and

$$ \int_{\Sigma} \tau \cdot n_{\Sigma} = \int_{\gamma} \sigma \cdot t = 0. $$

This gives us $\int_{\Omega} B_{\text{div}} \cdot \tau = 0$ for all $\tau \in \text{curl} H_0(\text{curl}, \Omega)$. Since $B_{\text{div}}$ is orthogonal to both $\nabla H^1(\Omega)$ and $\text{curl} H_0(\text{curl}, \Omega)$, we get that it is a harmonic field by Hodge decomposition (see Proposition 15 of Appendix A).

Finally, we set $B_{\text{div}} = \lambda B$, and want to prove that $\lambda = 1$. First, by plugging $\tau = B_{\text{div}}$ in Eq. (15), we get

$$ \|B_{\text{div}}\|^2 = 2\pi \int_{\Sigma} B_{\text{div}} \cdot n_{\Sigma}. $$

Now, using the jump condition on $u$ as defined in Eq. (9), and an integration by parts, we get

$$ \begin{aligned} \|B\|^2 &= \int_{\Omega \setminus \Sigma} B \cdot \nabla u \\ &= 2\pi \int_{\Sigma} B \cdot n_{\Sigma}. \end{aligned} $$

Now, on the one side we have

$$ 2\pi \int_{\Sigma} B_{\text{div}} \cdot n_{\Sigma} = 2\pi \lambda \int_{\Sigma} B \cdot n_{\Sigma}, $$

and on the other side

$$ \begin{aligned} 2\pi \int_{\Sigma} B_{\text{div}} \cdot n_{\Sigma} &= \|B_{\text{div}}\|^2 \\ &= \lambda^2 \|B\|^2 \\ &= 2\pi \lambda^2 \int_{\Sigma} B \cdot n_{\Sigma}. \end{aligned} $$

Since $B$ and $B_{\text{div}}$ are both nonzero, we come to the conclusion that $\lambda = 1$. $\square$

## 2.3 Vector potentials and Bevir–Gray formula

As was mentioned earlier, the numerical computation of the Biot–Savart operator can be very costly. As a consequence, we chose to compute the helicity of the normalized harmonic field by substituting $\text{BS}(B)$ by an appropriate vector potential of $B$. Indeed, since $\text{curl} \text{BS}(B) = B$, we know that any vector potential $A$ of $B$ can only differ from $\text{BS}(B)$ by the sum of a gradient and a harmonic field (see Eqs. (38) and (43) in Appendix A). Since vector fields of $H_0(\text{div}^0, \Omega)$ are orthogonal to gradient vector fields, we know that $\int_{\Omega} B \cdot A$ can differ from $H(B)$ only through the harmonic part of the difference between $A$ and $\text{BS}(B)$. However, this difference can be accounted for by modifying the formula of the magnetic helicity, giving a quantity which is invariant under a change of the vector potential. This is given by the well known Bevir–Gray formula [BG80] in toroidal domains, which was later generalized to a large class of non-simply connected domains in [MV19]. In our case, for a vector field $V$ in $H_0(\text{div}^0, \Omega)$ and $A$ any of its vector potentials, this invariant quantity is given by

$$ H(V) = \int_{\Omega} V \cdot A - \int_{\gamma} A \cdot t \int_{\gamma'} A \cdot t'. \quad (16) $$

Note that $\text{BS}(V)$ has zero circulation along $\gamma'$[CDG01, Section III.A], so that this invariant quantity does correspond to the usual Biot–Savart helicity.

From this formulation, a natural problem is to find good vector potentials, which are simple to study, both theoretically and numerically. Getting back to the helicity of the normalized harmonic field, we will study two natural vector potential candidates. First, a vector potential given by a classical vector Laplacian problem, which is orthogonal to $B$, therefore canceling the first term of Eq. (16). Second, a vector potential given in
[Val19], which is of zero circulation along the toroidal loop $\gamma'$, therefore canceling the second term. As we will see in Section 5, using these two vector potentials also allows us to stay in a finite elements exterior calculus setting throughout the numerical computation of the harmonic helicity, and its shape gradient.

The first vector potential, which we denote by $A^1$, is given by the following proposition.

**Proposition 2.** Let $B$ be the normalized harmonic field of $\Omega$ as defined in Section 2.2. There exists a unique $(A^1, u) \in H(\text{curl}, \Omega) \times H(\text{div}, \Omega)$ such that, for all $(\tau, v) \in H(\text{curl}, \Omega) \times H(\text{div}, \Omega)$

$$ \int_{\Omega} A^1 \cdot \tau = \int_{\Omega} \text{curl} \tau \cdot u, \qquad (17) $$

$$ \int_{\Omega} \text{curl} A^1 \cdot v + \int_{\Omega} (\text{div} u)(\text{div} v) = \int_{\Omega} B \cdot v. \qquad (18) $$

Furthermore, $A^1$ verifies

1. $A^1$ is in $H_0(\text{div}^0, \Omega)$,

2. $\text{curl } A^1 = B$,

3. $\int_{\Omega} A^1 \cdot B = 0$.

*Proof.* Since the second De Rham cohomology space of $\Omega$ is trivial, the space of harmonic 2-forms of $\Omega$ vanishes, and this variational problem is equivalent to [AFW06, Equation (7.1)] for the case $k = 2$. Existence and uniqueness is then given by [AFW06, Theorem 7.2].

To prove Item 1, we simply take $\tau$ in $\nabla H^1(\Omega)$ in Eq. (17). This gives us $\int_{\Omega} A^1 \cdot \tau = 0$, so that $A^1$ is orthogonal to $\nabla H^1(\Omega)$, and therefore is in $H_0(\text{div}^0, \Omega)$ by Eq. (40) in Appendix A. To find Item 2, we proceed with the same splitting analysis used in [AFW06]. We define using the Hodge decomposition Eq. (39), $u = u_{\nabla} + u_{\text{curl}}$ with $u_{\text{curl}} \in \text{curl } H(\text{curl}, \Omega)$ and $u_{\nabla} \in \nabla H_0^1(\Omega)$. We now want to prove that $u_{\nabla} = 0$. Choosing $v \in H(\text{div}^0, \Omega)^{\perp} \cap H(\text{div}, \Omega)$ in Eq. (18), we get

$$ \int_{\Omega} (\text{div } u_{\nabla}) (\text{div } v) = \int_{\Omega} B \cdot v = 0. $$

Using the Poincaré inequality from Proposition 16 of Appendix A, we get $u_{\nabla} = 0$, so that $\text{div } u = 0$ as claimed. Finally, since $\text{curl } B = 0$, we simply get by choosing $\tau = B$ in Eq. (17)

$$ \int_{\Omega} A^1 \cdot B = 0, $$

proving Item 3.

For the second vector potential $A^2$, we first need to introduce the following spaces.

$$ \mathcal{X}(\Omega) = \{V \in H(\text{curl}, \Omega) \mid \text{curl } V \cdot n = 0\}, $$

$$ \mathcal{Z}(\Omega) = \{V \in \mathcal{X}(\Omega) \mid \int_{\gamma'} V \cdot t' = 0\}. $$

We then introduce the second vector potential $A^2$, given by the following result proven in [Val19]

**Proposition 3.** The following problem has a unique solution. Find $(A^2, u) \in \mathcal{Z}(\Omega) \times \nabla H^1(\Omega)$ such that for all $(\tau, v) \in \mathcal{Z}(\Omega) \times \nabla H^1(\Omega)$

$$ \int_{\Omega} \text{curl } A^2 \cdot \text{curl } \tau + \int_{\Omega} u \cdot \tau = \int_{\Omega} B \cdot \text{curl } \tau, \qquad (19) \\ \int_{\Omega} A^2 \cdot v = 0. $$

Furthermore, we have $\text{curl } A^2 = B$ and $A^2 \in H_0(\text{div}^0, \Omega)$.

Of course, the two choices of potential vectors are related to each other. Since $A^1$ and $A^2$ are in $H_0(\text{div}^0, \Omega)$ and have the same curl, we get that $A^1 - A^2$ is in $\mathcal{K}(\Omega)$. The fact that $A^2$ has zero circulation along $\gamma'$, and that $B$ has a circulation of $2\pi$ allows us to find the relation

$$ A^2 = A^1 - \frac{1}{2\pi} \left( \int_{\gamma'} A^1 \cdot t' \right) B. \qquad (20) $$
### 3 Harmonic helicity and its shape derivative

As we have seen in the previous section, with each Lipschitz toroidal domain $\Omega$, we are able to associate a normalized harmonic field $B(\Omega)$ and a vector potential of $B(\Omega)$, denoted $A^2(\Omega)$, with zero circulation along $\gamma'$. As a consequence of the Bevir-Gray formula Eq. (16), the magnetic helicity of $B(\Omega)$ is then given by $H(B(\Omega)) = \int_{\Omega} B(\Omega) \cdot A^2(\Omega)$. In turn, this allows us to define the helicity of the normalized magnetic field of $\Omega$, which we refer to as the harmonic helicity of $\Omega$ for simplicity. As was noted in Remark 1, $B(\Omega)$ is actually defined as a function of $\Omega$ only up to a sign. However, since the helicity is a quadratic form, this sign indetermination is not relevant, and the harmonic helicity is well-defined as a shape functional.

**Definition 1.** Let $\Omega$ be a Lipschitz toroidal domain. Then, the harmonic helicity of $\Omega$ is defined as

$$ \mathcal{H}(\Omega) = \int_{\Omega} B(\Omega) \cdot A^2(\Omega), $$

where $B(\Omega)$ and $A^2(\Omega)$ are given by the solutions to Eq. (9) and Eq. (19).

The aim of this section is to study how the harmonic helicity varies as a function of the domain $\Omega$. More precisely, we prove that the harmonic helicity is shape Fréchet differentiable under Lipschitz deformation, and we give a formula for its shape derivative.

Before studying the shape differentiability of $\mathcal{H}$, we state the following properties of harmonic helicity. The proof will be given in Section 3.1 as we will need to define ways to pullback vector fields beforehand.

**Proposition 4.** The following scaling and symmetry properties hold:

* Given $\lambda > 0$, we have $\mathcal{H}(\lambda\Omega) = \lambda\mathcal{H}(\Omega)$.

* Given a planar reflection $R \in O_3$, we have $\mathcal{H}(R\Omega) = -\mathcal{H}(\Omega)$.

**Corollary 1.** If $\Omega$ admits a planar symmetry, we have $\mathcal{H}(\Omega) = 0$. In particular, this is the case if $\Omega$ is axisymmetric.

As was mentioned in the introduction, another notion for the helicity of a domain which was studied in [Can+00; Can+99; Val19; Mon23] is the maximal value of the helicity on the unit $L^2$ sphere for fields in $H_0(\text{div}^0, \Omega)$, that is

$$ \tilde{\mathcal{H}}(\Omega) = \sup \{ \text{H}(V) \mid V \in H_0(\text{div}^0, \Omega), \|V\|_{L^2} = 1 \}. $$

Since the helicity is quadratic, it is therefore clear that we have the following inequality

$$ \mathcal{H}(\Omega) \le \|B(\Omega)\|_{L^2}^2 \tilde{\mathcal{H}}(\Omega). $$

To remove the $L^2$ norm in this inequality, one may normalize the harmonic field using the $L^2$ norm instead of the circulation. Although this makes the comparison between the two quantities more direct, this normalization is less convenient to work with in our context. We note however that one may easily obtain the shape derivative of $\|B(\Omega)\|_{L^2}$ using the methods we will introduce in this section, so that this change of normalization may be done without issues for computing the shape derivative of harmonic helicity. We also note that, although $\tilde{\mathcal{H}}$ is known to be bounded when fixing the volume (see [Can+00, Theorem E]), it is not immediately clear whether this is the case for $\mathcal{H}$. The optimization problems we will consider numerically will be described more precisely in Section 5.

In order to write the shape derivative of the harmonic helicity as a surface integral, we need to assume that $\Omega$ is s-regular for some $s > 1/2$ [AFW06, Section 7.7], that is we have the continuous embedding

$$ H(\text{curl}, \Omega) \cap H_0(\text{div}, \Omega) \hookrightarrow H^s(\Omega). \qquad (21) $$

Indeed, we need $B(\Omega)$ and $A^2(\Omega)$ to have traces in $L^2(\partial\Omega)^3$. For example $\Omega$ being Lipschitz-polyhedral or $C^{1,1}$ is sufficient [Amr+98, Prop. 3.7 or Th. 2.9].

**Theorem 1.** Let $\Omega$ be a s-regular toroidal domain, with $s > 1/2$, and $\theta$ a vector field in $W^{1,\infty}(\mathbb{R}^3)$ with $\|D\theta\|_{L^\infty} < 1$. Denoting $\Omega_\theta = (I + \theta)(\Omega)$, we have

$$ \mathcal{H}(\Omega_\theta) = \mathcal{H}(\Omega) + \mathcal{H}'(\Omega; \theta) + o(\|\theta\|_{W^{1,\infty}}), \qquad (22) $$

where

$$ \mathcal{H}'(\Omega; \theta) = 2 \int_{\partial\Omega} (B(\Omega) \cdot A^2(\Omega)) (\theta \cdot n). \qquad (23) $$
In order to prove this theorem, we begin by introducing some transformations in Section 3.1. These transformations, which correspond to usual pullbacks in the language of differential forms, allow us to transform functions and vector fields on the deformed domain $\Omega_\theta$ onto the fixed domain $\Omega$. Since these transformations have good commutation properties with the differential operators, we are then able to use them to differentiate the variational formulations of $B(\Omega_\theta)$ and $A^2(\Omega_\theta)$ in Section 3.2. Once this is done, we can prove Theorem 1 in Section 3.3 by pulling back the integral defining $\mathcal{H}$ onto a fixed domain, and using the differentiated vector fields from the previous section.

Throughout the rest of this section, $\theta$ denotes a vector field in $W^{1,\infty}(\mathbb{R}^3)$ with $\|D\theta\|_{L^\infty} < 1$, $\Omega_\theta$ the deformed domain $(I + \theta)(\Omega)$, $D\theta$ denotes the Jacobian matrix of $\theta$ in Cartesian coordinates, and $J_\theta = \det(I + D\theta)$ the determinant of the Jacobian of the transformation. More generally, given a function (resp. vector field) $u$, we denote by $Du$ its Jacobian matrix in Cartesian coordinates, which is a row vector (resp. square matrix) at each point where it is defined. To have notations which are less cumbersome, we will also denote $B = B(\Omega)$, $B_\theta = B(\Omega_\theta)$ the solutions to Eq. (15) in $\Omega$ and $\Omega_\theta$ respectively, and $A = A^2(\Omega)$, $A_\theta = A^2(\Omega_\theta)$ the solutions to Eq. (19) in $\Omega$ and $\Omega_\theta$ respectively. We also denote by $n_\theta$ the unit outward pointing vector field on $\partial\Omega_\theta$, and define $\gamma_\theta$, $\gamma'_\theta$, $\Sigma_\theta$ and $\Sigma'_\theta$ in the deformed domain $\Omega_\theta$ as the image by $(I + \theta)$ of the corresponding geometrical constructions given in Section 2.2, with tangent and normal vector fields $t_\theta$, $t'_\theta$, $n_{\Sigma_\theta}$, $n_{\Sigma'_\theta}$ respectively.

## 3.1 Pullback of the De Rham complex on a fixed domain

In this section, we define maps which allow us to transform elements of our functional spaces on a deformed domain $\Omega_\theta$ back to $\Omega$. These transformations correspond to pullbacks in the language of differential forms, and the corresponding formulas are therefore similar to the Piola mappings used in finite elements exterior calculus. This non-naive treatment is needed as composition by $I + \theta$ does not map $H(\text{curl}, \Omega_\theta)$ (resp. $H(\text{div}, \Omega_\theta)$) to $H(\text{curl}, \Omega)$ (resp. $H(\text{div}, \Omega)$).

**Definition 2.** Let $u_0$, $u_1$, $u_2$ and $u_3$ be in $H^1(\Omega_\theta)$, $H(\text{curl}, \Omega_\theta)$, $H(\text{div}, \Omega_\theta)$ and $L^p(\Omega_\theta)$ respectively, with $p$ being 1 or 2. We define

$$ 
\begin{aligned}
\Phi_{\theta}^{0} u_{0} &= u_{0} \circ (I + \theta), \\
\Phi_{\theta}^{1} u_{1} &= (I + D\theta^{T}) u_{1} \circ (I + \theta), \\
\Phi_{\theta}^{2} u_{2} &= J_{\theta}(I + D\theta)^{-1} u_{2} \circ (I + \theta), \\
\Phi_{\theta}^{3} u_{3} &= J_{\theta} u_{3} \circ (I + \theta).
\end{aligned}
$$

Before giving the key properties of these transformations, we start by stating the following useful algebraic identities. These can be found by direct computations.

**Lemma 1.** We have the following identities.

$$ 
\begin{aligned}
(\Phi_{\theta}^{0} u) (\Phi_{\theta}^{3} v) &= \Phi_{\theta}^{3} (uv) & \forall u \in H^{1}(\Omega_{\theta}), v \in L^{2}(\Omega_{\theta}), \\
(\Phi_{\theta}^{1} u) \cdot (\Phi_{\theta}^{2} v) &= \Phi_{\theta}^{3} (u \cdot v) & \forall u \in H(\text{curl}, \Omega_{\theta}), v \in H(\text{div}, \Omega_{\theta}), \\
(\alpha(\theta)\Phi_{\theta}^{1} u) \cdot (\Phi_{\theta}^{1} v) &= \Phi_{\theta}^{3} (u \cdot v) & \forall u \in H(\text{curl}, \Omega_{\theta}), v \in H(\text{curl}, \Omega_{\theta}), \\
(\alpha(\theta)^{-1}\Phi_{\theta}^{2} u) \cdot (\Phi_{\theta}^{2} v) &= \Phi_{\theta}^{3} (u \cdot v) & \forall u \in H(\text{div}, \Omega_{\theta}), v \in H(\text{div}, \Omega_{\theta}),
\end{aligned}
$$

where $\alpha(\theta) = J_\theta(I + D\theta)^{-1}(I + D\theta^T)^{-1}$.

**Proposition 5.** The diagram

$$ 
\begin{aligned}
H^1(\Omega_\theta) &\xrightarrow{\nabla} H(\text{curl}, \Omega_\theta) \xrightarrow{\text{curl}} H(\text{div}, \Omega_\theta) \xrightarrow{\text{div}} L^2(\Omega_\theta) \\
\Phi_\theta^0 \Big|_{\Omega_\theta} &\xrightarrow{\nabla} \Phi_\theta^1 \Big|_{\Omega_\theta} \xrightarrow{\text{curl}} \Phi_\theta^2 \Big|_{\Omega_\theta} \xrightarrow{\text{div}} \Phi_\theta^3 \Big|_{\Omega_\theta} \\
H^1(\Omega) &\xrightarrow{\nabla} H(\text{curl}, \Omega) \xrightarrow{\text{curl}} H(\text{div}, \Omega) \xrightarrow{\text{div}} L^2(\Omega)
\end{aligned}
$$

is commutative. Furthermore, the same diagram can be made with the corresponding traceless spaces.

*Proof.* For the proof of the first statement, we choose smooth functions or vector fields for simplicity (both for $u$ and $\theta$). The general cases can then be obtained by density arguments.

We begin by proving $\nabla \circ \Phi_\theta^0 = \Phi_\theta^1 \circ \nabla$. We choose $u$ in $C^\infty(\bar{\Omega}_\theta)$, and compute using the chain rule

$$ D(\Phi_{\theta}^{0} u) = Du \circ (I + \theta)(I + D\theta), $$
Since, in Cartesian coordinates, the gradient of a function is the transpose of its Jacobian matrix, we get

$$ 
\begin{aligned}
\nabla \Phi_{\theta}^{0} u &= (I + D\theta^{T})\nabla u \circ (I + \theta) \\
&= \Phi_{\theta}^{1} \nabla u.
\end{aligned}
$$

Now, we prove that $\mathrm{curl} \circ \Phi_{\theta}^{1} = \Phi_{\theta}^{2} \circ \mathrm{curl}$ by taking $u$ in $C^{\infty}(\bar{\Omega}_{\theta})^{3}$. We have

$$ 
\begin{aligned}
D(\Phi_{\theta}^{1} u) &= D((I + D\theta^{T})u \circ (I + \theta)) \\
&= D(D\theta^{T})u \circ (I + \theta) + (I + D\theta^{T})Du \circ (I + \theta)(I + D\theta),
\end{aligned}
$$

where $D(D\theta^{T})u$ is a symmetric matrix given by

$$ (D(D\theta^T)u)_{i,j} = \sum_{k=1}^{3} \frac{\partial^2 \theta^k}{\partial x^i \partial x^j} u_k, $$

Therefore, we find

$$ D(\Phi_{\theta}^{1} u) - (D(\Phi_{\theta}^{1} u))^{T} = (I + D\theta^{T})(Du - Du^{T}) \circ (I + \theta)(I + D\theta). $$

Now, defining $\mathrm{cr} : \mathfrak{so}(3) \to \mathbb{R}^{3}$, where $\mathfrak{so}(3)$ is the space of skew-symmetric matrices, by

$$ \mathrm{cr} \begin{pmatrix} 0 & -v_3 & v_2 \\ v_3 & 0 & -v_1 \\ -v_2 & v_1 & 0 \end{pmatrix} = \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix}, $$

we can prove that, for any invertible matrix $B$ and skew-symmetric matrix $A$, $\mathrm{cr}(B^{T}AB) = \det(B)B^{-1}\mathrm{cr}(A)$. Since $\mathrm{curl} u = \mathrm{cr}(Du - Du^{T})$, we find

$$ 
\begin{aligned}
\mathrm{curl} \Phi_{\theta}^{1} u &= J_{\theta}(I + D\theta)^{-1}\mathrm{curl} u \circ (I + \theta) \\
&= \Phi_{\theta}^{2}\mathrm{curl} u.
\end{aligned}
$$

For the last commutativity relation $\mathrm{div} \circ \Phi_{\theta}^{2} = \Phi_{\theta}^{3} \circ \mathrm{div}$, we proceed by duality. Take $u$ in $C^{\infty}(\bar{\Omega}_{\theta})^{3}$ and $v$ in $C_{0}^{\infty}(\Omega_{\theta})$. First, we notice that the $\Phi_{\theta}^{k}$ are isomorphisms, as $(\Phi_{\theta}^{k})^{-1} = \Phi_{(I+\theta)^{-1}-I}^{k}$. Therefore,

$$ 
\begin{aligned}
\int_{\Omega} (\mathrm{div} \Phi_{\theta}^{2} u) v &= - \int_{\Omega} (\Phi_{\theta}^{2} u) \cdot \nabla v \\
&= - \int_{\Omega} (\Phi_{\theta}^{2} u) \cdot [\Phi_{\theta}^{1} (\nabla (\Phi_{\theta}^{0})^{-1} v)],
\end{aligned}
$$

where we used $\Phi_{\theta}^{1} \nabla v = \nabla \Phi_{\theta}^{0} v$. Using Lemma 1, we then find

$$ 
\begin{aligned}
\int_{\Omega} (\mathrm{div} \Phi_{\theta}^{2} u) v &= - \int_{\Omega} \Phi_{\theta}^{3} (u \cdot \nabla (\Phi_{\theta}^{0})^{-1} v) \\
&= - \int_{\Omega_{\theta}} u \cdot \nabla (\Phi_{\theta}^{0})^{-1} v \\
&= \int_{\Omega_{\theta}} \mathrm{div} u \cdot (\Phi_{\theta}^{0})^{-1} v \\
&= \int_{\Omega} (\Phi_{\theta}^{3} \mathrm{div} u) v.
\end{aligned}
$$

Now, we prove that traceless functions or vector fields are preserved by the maps $\Phi_{\theta}^{k}$. For $k = 0$, we simply take $u \in C_{0}^{\infty}(\Omega)$, and verify that $u \circ (I + \theta)$ is in $H_{0}^{1}(\Omega)$. We then find the desired result by density.

The proofs for $k = 1$ and $k = 2$ being similar, we only write the first one. Take $u$ in $H_{0}(\mathrm{curl}, \Omega_{\theta})$ and $v$ in $H(\mathrm{curl}, \Omega)$. Then, we have

$$ 
\begin{aligned}
\int_{\partial\Omega} ((\Phi_{\theta}^{1} u) \times n) \cdot v &= \int_{\Omega} \mathrm{curl} (\Phi_{\theta}^{1} u) \cdot v - \int_{\Omega} (\Phi_{\theta}^{1} u) \cdot \mathrm{curl} v \\
&= \int_{\Omega} \Phi_{\theta}^{2} (\mathrm{curl} u) \cdot \Phi_{\theta}^{1} ((\Phi_{\theta}^{1})^{-1} v) - \int_{\Omega} (\Phi_{\theta}^{1} u) \cdot \Phi_{\theta}^{2} (\mathrm{curl} (\Phi_{\theta}^{1})^{-1} v) \\
&= \int_{\Omega_{\theta}} \mathrm{curl} u \cdot ((\Phi_{\theta}^{1})^{-1} v) - \int_{\Omega_{\theta}} u \cdot (\mathrm{curl} (\Phi_{\theta}^{1})^{-1} v) \\
&= 0,
\end{aligned}
$$
so that $\Phi_{\theta}^{1}u$ is in $H_{0}(\text{curl}, \Omega)$ as claimed.

As the $\Phi_{\theta}^{k}$ transformations correspond to pullbacks of k-forms, it is not surprising that $\Phi_{\theta}^{1}$ preserves circulations, and $\Phi_{\theta}^{2}$ preserves fluxes. This result is given by the following lemma.

**Lemma 2.** Let $u$ be in $H(\text{curl}, \Omega_{\theta})$ with $\text{curl } u \cdot n = 0$, and $v$ be in $H_{0}(\text{div}, \Omega_{\theta})$. We then have the following identities

$$ \int_{\gamma'_{\theta}} u \cdot t'_{\theta} = \int_{\gamma'} \Phi_{\theta}^{1} u \cdot t', \\ \int_{\Sigma_{\theta}} v \cdot n_{\Sigma_{\theta}} = \int_{\Sigma} \Phi_{\theta}^{2} v \cdot n_{\Sigma}. $$

*Proof.* For the first equality, we suppose that $u$ is smooth so that we can study the circulation in the usual sense. The general case then follows from a density argument, and the continuity of the circulation as defined in [Alo+18, Section 2] with respect to the $H(\text{curl}, \Omega)$ norm. We identify $\gamma'$ with a Lipschitz embedding $\gamma' : S^{1} \to \partial\Omega$, so that the derivative $\dot{\gamma}'(s)$ is defined for a.e. $s \in S^{1}$. $\gamma'_{\theta}$ is then identified as the embedding $(I + \theta) \circ \gamma'$, for which the derivative is given by

$$ \dot{\gamma}'_{\theta} = ((I + D\theta) \circ \gamma') \dot{\gamma}'. $$

The circulation of $u$ along $\gamma'_{\theta}$ is then given by

$$ \int_{\gamma'_{\theta}} u \cdot t'_{\theta} = \int_{S^{1}} [u \circ \gamma'_{\theta}] (s) \cdot \dot{\gamma}'_{\theta}(s) ds. $$

Using the formula for $\gamma'_{\theta}$ and its derivative as an embedding, we get

$$ \begin{aligned} \int_{\gamma'_{\theta}} u \cdot t'_{\theta} &= \int_{S^{1}} [u \circ (I + \theta) \circ \gamma'] (s) \cdot [((I + D\theta) \circ \gamma') \dot{\gamma}'] (s) ds \\ &= \int_{S^{1}} [(I + D\theta^{T}) u \circ (I + \theta)] \circ \gamma'(s) \cdot \dot{\gamma}'(s) ds \\ &= \int_{\gamma'} \Phi_{\theta}^{1} u \cdot t'. \end{aligned} $$

For the second equality, we choose $v \in H_{0}(\text{div}, \Omega_{\theta})$ and use the following equality given in [Alo+18, Lemma 1] for $\psi \in H^{1}(\Omega_{\theta} \setminus \Sigma_{\theta})$

$$ \int_{\Sigma_{\theta}} v \cdot n_{\Sigma_{\theta}} [[\psi]]_{\Sigma_{\theta}} = \int_{\Omega_{\theta} \setminus \Sigma_{\theta}} v \cdot \nabla \psi + \int_{\Omega_{\theta} \setminus \Sigma_{\theta}} (\text{div } v) \psi. \quad (24) $$

Taking $u_{\theta}$ to be the solution of Eq. (9) in $\Omega_{\theta}$, we get

$$ \begin{aligned} 2\pi \int_{\Sigma_{\theta}} v \cdot n_{\Sigma_{\theta}} &= \int_{\Omega_{\theta} \setminus \Sigma_{\theta}} v \cdot \nabla u_{\theta} + \int_{\Omega_{\theta} \setminus \Sigma_{\theta}} (\text{div } v) u_{\theta} \\ &= \int_{\Omega \setminus \Sigma} (\Phi_{\theta}^{2} v) \cdot \nabla (\Phi_{\theta}^{0} u_{\theta}) + \int_{\Omega \setminus \Sigma} (\text{div } \Phi_{\theta}^{2} v) (\Phi_{\theta}^{0} u_{\theta}) \\ &= \int_{\Sigma} (\Phi_{\theta}^{2} v) \cdot n_{\Sigma} [[\Phi_{\theta}^{0} u_{\theta}]]_{\Sigma}. \end{aligned} $$

From the definition of $\Phi_{\theta}^{0}$, we then see that $[[\Phi_{\theta}^{0} u_{\theta}]]_{\Sigma} = [[u_{\theta}]]_{\Sigma_{\theta}} \circ (I + \theta) = 2\pi$, which proves the desired equality. $\square$

We know from Lemma 1 that when we pullback products in $H(\text{curl}, \Omega)$ or $H(\text{div}, \Omega)$, there is a non-homogeneous term $\alpha(\theta)$ which appears. We will have to consider such products in the next sections to differentiate the harmonic helicity. The next lemma shows how $\theta \mapsto \alpha(\theta)$ behaves under differentiation when integrated against products of vector fields.

**Lemma 3.** The mapping $\theta \mapsto \alpha(\theta)$ from $W^{1,\infty}(\mathbb{R}^{3})^{3}$ to $L^{\infty}(\mathbb{R}^{3}; \mathcal{M}_{3}(\mathbb{R}))$ defined in Lemma 1 is smooth, and if $u, v$ are in $H(\text{curl}, \Omega) \cap H_{0}(\text{div}, \Omega)$, its differential at 0 verifies

$$ \int_{\Omega} (\alpha'(0; \theta)u) \cdot v = \int_{\Omega} \text{div} u\theta \cdot v + \int_{\Omega} u \cdot \theta \text{div} v + \int_{\Omega} u \times \theta \cdot \text{curl} v - \int_{\Omega} \text{curl} u \times \theta \cdot v + \int_{\partial\Omega} (u \cdot v)(\theta \cdot n). $$
Proof. We recall that $\alpha(\theta) = J_0(I + D\theta)^{-1} (I + D\theta^T)^{-1}$. To find that $\alpha$ is smooth, we simply notice that the maps defined on the unit ball of $M_3(\mathbb{R})$ given by $A \mapsto (I + A)^{-1}$, $A \mapsto (I + A^T)^{-1}$ and $A \mapsto \det(I + A)$ are smooth. We also have that $\theta \mapsto D\theta$ is linear and bounded from $W^{1,\infty}(\mathbb{R}^3)^3$ to $L^\infty(\mathbb{R}^3; M_3(\mathbb{R}))$. We conclude by composition that $\alpha$ is smooth.

Using that the differential of the determinant and the inverse at the identity are, respectively, the trace and minus the identity, we find

$$ \alpha'(0; \theta) = \operatorname{div} \theta I - D\theta - D\theta^T. $$

To prove the final point, we choose $u$ and $v$ in $H^1(\Omega)^3$ such that $u \cdot n = v \cdot n = 0$. Then, $u \cdot \theta$ and $u \times \theta$ are in $H^1$, and we have

$$ \begin{aligned} \nabla(u \cdot \theta) &= Du^T \theta + D\theta^T u, \\ \operatorname{curl}(u \times \theta) &= Du\theta - D\theta u - \operatorname{div} u\theta + \operatorname{div} \theta u. \end{aligned} $$

Combining these formulas, we find

$$ \begin{aligned} \nabla(u \cdot \theta) - \operatorname{curl}(u \times \theta) &= -Du\theta + D\theta u + \operatorname{div} u\theta - \operatorname{div} \theta u + Du^T \theta + D\theta^T u \\ &= -\alpha'(0; \theta)u - (Du - Du^T)\theta + \operatorname{div} u\theta \\ &= -\alpha'(0; \theta) - \operatorname{curl} u \times \theta + \operatorname{div} u\theta, \end{aligned} $$

so that

$$ \alpha'(0; \theta)u = \operatorname{div} u\theta + \operatorname{curl}(u \times \theta) - \operatorname{curl} u \times \theta \cdot v - \nabla(u \cdot \theta). $$

Now, using integration by parts

$$ \begin{aligned} \int_{\Omega} \alpha'(0; \theta) u \cdot v &= \int_{\Omega} \operatorname{div} u\theta \cdot v + \operatorname{curl}(u \times \theta) \cdot v - \operatorname{curl} u \times \theta \cdot v - \nabla(u \cdot \theta) \cdot v \\ &= \int_{\Omega} \operatorname{div} u\theta \cdot v + \int_{\Omega} u \times \theta \cdot \operatorname{curl} v - \int_{\partial\Omega} (u \times \theta) \times n \cdot v - \int_{\Omega} \operatorname{curl} u \times \theta \cdot v + \int_{\Omega} u \cdot \theta \operatorname{div} v \\ &= \int_{\Omega} \operatorname{div} u\theta \cdot v + \int_{\Omega} u \times \theta \cdot \operatorname{curl} v + \int_{\partial\Omega} (u \cdot v)(\theta \cdot n) - \int_{\Omega} \operatorname{curl} u \times \theta \cdot v + \int_{\Omega} u \cdot \theta \operatorname{div} v, \end{aligned} $$

where we used the formula $(a \times b) \times c = (c \cdot a)b - (c \cdot b)a$ and the fact that $u \in H_0(\operatorname{div}, \Omega)$. Finally, by density, taking any $u$ and $v$ in $H(\operatorname{curl}, \Omega) \cap H_0(\operatorname{div}, \Omega)$, we get the desired formula. Note that we need s-regularity Eq. (21) with $s > 1/2$, to ensure that $u, v \in L^2(\partial\Omega)^3$. $\square$

Before studying the shape differentiability of harmonic helicity in the next sections, we prove Proposition 4.

Proof. We begin by proving the first point. Let $u(\Omega)$ and $u(\lambda\Omega)$ be the solution to Eq. (9) in $\Omega$ and $\lambda\Omega$ respectively. Then, it is clear that $u(\lambda\Omega) = u(\cdot/\lambda)$. Therefore, we have $B(\lambda\Omega) = \nabla u_\lambda = \frac{1}{\lambda}B(\Omega)(\cdot/\lambda)$. Similarly, we have $A^2(\lambda\Omega) = \frac{1}{\lambda}A(\Omega)(\cdot/\lambda)$. Writing

$$ \begin{aligned} \mathcal{H}(\lambda\Omega) &= \int_{\lambda\Omega} B(\lambda\Omega) \cdot A^2(\lambda\Omega) \\ &= \lambda^3 \int_{\Omega} \left( \frac{1}{\lambda} B(\Omega) \right) \cdot \left( \frac{1}{\lambda} A^2(\Omega) \right) \\ &= \lambda \mathcal{H}(\Omega) \end{aligned} $$

For the second point, let $R$ be a planar reflection. We introduce as in Definition 2

$$ \begin{aligned} \Phi_R^1 u &= R^T u \circ R, \\ \Phi_R^2 u &= (\operatorname{det} R) R^{-1} u \circ R, \end{aligned} $$

for $u$ a vector field of $R\Omega$. Since $R$ is a planar reflection, we have $\operatorname{det} R = -1$ and $R^{-1} = R^T$ so that $\Phi_R^2 u = -\Phi_R^1 u$. Let $u(\Omega)$ and $u(R\Omega)$ be the solutions to Eq. (9) in $\Omega$ and $R\Omega$ respectively. Then, since $R$ is an isometry, we have

$$ u(\Omega) = u(R\Omega) \circ R, $$
so that using the definition of the normalized harmonic field and Proposition 5, we have $B(\Omega) = \Phi_R^1 B(R\Omega)$. We now prove that $A^2(\Omega) = \Phi_R^2 A^2(R\Omega)$. First, since $\Phi_R^2 = -\Phi_R^1$, it is clear from Proposition 5 that $\Phi_R^2 A^2(R\Omega)$ is in $H_0(\operatorname{div}^0, \Omega)$. Using once again Proposition 5, we also have

$$ 
\begin{aligned}
\operatorname{curl} \Phi_R^2 A^2(R\Omega) &= -\operatorname{curl} \Phi_R^1 A^2(R\Omega) \\
&= -\Phi_R^2 B(R\Omega) \\
&= \Phi_R^1 B(R\Omega) \\
&= B(\Omega).
\end{aligned}
 $$

Moreover, using Lemma 2, we know that $\Phi_R^2 A^2(R\Omega) = -\Phi_R^1 A^2(R\Omega)$ is in $\mathcal{Z}(\Omega)$. We then deduce from [Val19, Proposition 2] that $\Phi_R^2 A^2(R\Omega) = A^2(\Omega)$. Finally, since $R$ is an isometry, we have

$$ 
\begin{aligned}
\mathcal{H}(R\Omega) &= \int_{R\Omega} B(R\Omega) \cdot A^2(R\Omega) \\
&= \int_{\Omega} (R^T B(R\Omega) \circ R) \cdot (R^T A^2(R\Omega) \circ R), \\
&= -\int_{\Omega} \Phi_R^1 B(R\Omega) \cdot \Phi_R^2 A^2(R\Omega) \\
&= -\mathcal{H}(\Omega).
\end{aligned}
 $$

□

## 3.2 Differentiability of the PDEs

In this section we prove the differentiability of $\theta \mapsto \Phi_\theta^2 B_\theta$ and $\theta \mapsto \Phi_\theta^1 A_\theta$. As is classical in shape variation of PDEs (see for example [HP18, Chapter 5]), this is done using an implicit function argument on a pulled back version of the variational formulations. Although it is quite common in the literature to recover Eulerian derivatives from this, that is, the differential of $\theta \mapsto B_\theta$ and $\theta \mapsto A_\theta$ directly, we chose to skip this step here. Indeed, we will see that this step is not necessary to recover the shape derivative of the harmonic helicity. Furthermore, these Eulerian derivatives often suffer from a loss of regularity compared to the material ones, and satisfy affine variational formulations which are less practical to deal with in our case.

**Proposition 6.** Let $(B_\theta, u_\theta)$ be the solution to the following variational problem. Find $(B_\theta, u_\theta) \in H_0(\operatorname{div}, \Omega_\theta) \times L_0^2(\Omega_\theta)$ such that, for all $(\tau, v) \in H_0(\operatorname{div}, \Omega_\theta) \times L_0^2(\Omega_\theta)$ we have

$$ \int_{\Omega_\theta} B_\theta \cdot \tau + \int_{\Omega_\theta} u_\theta \operatorname{div} \tau = 2\pi \int_{\Sigma_\theta} \tau \cdot n_{\Sigma_\theta}, \qquad (25) $$

$$ \int_{\Omega_\theta} (\operatorname{div} B_\theta) v = 0. \qquad (26) $$

Then, $\theta \mapsto (\Phi_\theta^2 B_\theta, \Phi_\theta^0 u_\theta)$ is $C^1$ in a neighborhood of 0, and its differential at zero $(B', u')$ verifies

$$ 
\begin{aligned}
\int_{\Omega} B' \cdot \tau + \int_{\Omega} u' \operatorname{div} \tau &= \int_{\Omega} (\alpha'(0; \theta) B) \cdot \tau, \\
\int_{\Omega} (\operatorname{div} B') v &= 0,
\end{aligned}
 \qquad (27) $$

for all $(\tau, v)$ in $H_0(\operatorname{div}, \Omega) \times L_0^2(\Omega)$.

**Proof.** We know from Proposition 1 that $(B_\theta, u_\theta)$ is defined uniquely. Pulling back the integrals of Eq. (25) onto $\Omega$, and using Lemmas 1 to 3, we get

$$ 
\begin{aligned}
\int_{\Omega} \Phi_{\theta}^{3}(B_{\theta} \cdot \tau) + \int_{\Omega} \Phi_{\theta}^{3}(u_{\theta} \operatorname{div} \tau) &= 2\pi \int_{\Sigma} (\Phi_{\theta}^{2}\tau) \cdot n_{\Sigma}, \\
\int_{\Omega} (\alpha(\theta)^{-1}\Phi_{\theta}^{2}B_{\theta}) \cdot (\Phi_{\theta}^{0}\tau) + \int_{\Omega} (\Phi_{\theta}^{0}u_{\theta}) (\operatorname{div} \Phi_{\theta}^{2}\tau) &= 2\pi \int_{\Sigma} (\Phi_{\theta}^{2}\tau) \cdot n_{\Sigma}.
\end{aligned}
 $$

Similarly, we get from Eq. (26)

$$ \int_{\Omega} (\operatorname{div} \Phi_{\theta}^{2} B_{\theta}) (\Phi_{\theta}^{0} v) = 0. $$
Of course, since the $\Phi_{\theta}^k$ define isomorphisms, we can take test functions $(\tau, v)$ in $H_0(\text{div}, \Omega) \times L_0^2(\Omega)$. We now define

$$F : W^{1,\infty}(\mathbb{R}^3)^3 \times (H_0(\text{div}, \Omega) \times L_0^2(\Omega)) \to (H_0(\text{div}, \Omega) \times L_0^2(\Omega))'$$

by

$$F(\theta; \sigma, u)(\tau, v) = \int_{\Omega} (\alpha(\theta)^{-1}\sigma) \cdot \tau + \int_{\Omega} u \operatorname{div} \tau + \int_{\Omega} (\operatorname{div} \sigma) v - 2\pi \int_{\Sigma} \tau \cdot n_{\Sigma},$$

so that $(B_{\theta}, u_{\theta})$ solving Eqs. (25) and (26) is equivalent to $F(\theta; \Phi_{\theta}^2 B_{\theta}, \Phi_{\theta}^0 u_{\theta}) = 0$.

Now, we know from Lemma 3 that $\theta \mapsto \alpha(\theta)^{-1}$ is $C^1$. Since $(\sigma, u) \mapsto F(\theta; \sigma, u)$ is linear and bounded, we deduce that $F$ is $C^1$. Furthermore, denoting by $D_{\sigma,u}F$ the differential of $F$ with respect to the $(\sigma, u)$ variables,

$$D_{\sigma,u}F(0; B_0, u_0)(\sigma', u')(\tau, v) = \int_{\Omega} \sigma' \cdot \tau + \int_{\Omega} u' \operatorname{div} \tau + \int_{\Omega} (\operatorname{div} \sigma') v.$$

We therefore find that $D_{\sigma,u}F(0; B_0, u_0)$ is an isomorphism from $H_0(\text{div}, \Omega) \times L_0^2(\Omega)$ to $(H_0(\text{div}, \Omega) \times L_0^2(\Omega))'$ by the same inf-sup inequalities used to prove the well-posedness of Proposition 1. Using the implicit function theorem, we deduce that for small enough $\theta$, there is a unique $C^1$ mapping $\theta \mapsto (\sigma(\theta), u(\theta))$ such that $F(\theta; \sigma(\theta), u(\theta)) = 0$. From uniqueness, we find that $(\sigma(\theta), u(\theta)) = (\Phi_{\theta}^2 B_{\theta}, \Phi_{\theta}^0 u_{\theta})$.

Now, to get a variational formulation, we simply differentiate $F(\theta; \Phi_{\theta}^2 B_{\theta}, \Phi_{\theta}^0 u_{\theta}) = 0$ at $\theta = 0$. Denoting $\Phi_{\theta}^2 B_{\theta} = B + B' + o(||\theta||_{C^{1,1}})$ and $\Phi_{\theta}^0 u_{\theta} = u + u' + o(||\theta||_{C^{1,1}})$, this gives us

$$\int_{\Omega} B' \cdot \tau + \int_{\Omega} u' \operatorname{div} \tau + \int_{\Omega} (\operatorname{div} B') v = \int_{\Omega} (\alpha'(0; \theta) B_0) \cdot \tau,$$

which concludes the proof.

□

**Proposition 7.** Let $(A_{\theta}, u_{\theta})$ be the solution to the following variational problem.

Find $(A_{\theta}, u_{\theta})$ in $\mathcal{Z}(\Omega_{\theta}) \times \nabla H^1(\Omega_{\theta})$ such that, for all $(\tau, v)$ in $\mathcal{Z}(\Omega_{\theta}) \times \nabla H^1(\Omega_{\theta})$ we have

$$\int_{\Omega_{\theta}} \operatorname{curl} A_{\theta} \cdot \operatorname{curl} \tau + \int_{\Omega_{\theta}} u_{\theta} \cdot \tau = \int_{\Omega_{\theta}} B_{\theta} \cdot \operatorname{curl} \tau, \qquad (28)$$

$$\int_{\Omega_{\theta}} A_{\theta} \cdot v = 0. \qquad (29)$$

Then, $\theta \mapsto (\Phi_{\theta}^1 A_{\theta}, \Phi_{\theta}^2 u_{\theta})$ is $C^1$ in a neighborhood of 0.

*Proof.* We proceed with an implicit function theorem argument, similar to the one used for Proposition 6. First, we note that from Proposition 5 and Lemma 2, we get $\Phi_{\theta}^1 \mathcal{Z}(\Omega_{\theta}) = \mathcal{Z}(\Omega)$, so that the functional spaces of Eqs. (28) and (29) are preserved by the pullbacks.

Pulling back Eqs. (28) and (29) onto $\Omega$ and using Lemma 1, we get

$$\int_{\Omega} (\alpha(\theta)^{-1} \operatorname{curl} \Phi_{\theta}^{1} A_{\theta}) \cdot (\operatorname{curl} \Phi_{\theta}^{1} \tau) + \int_{\Omega} (\alpha(\theta) \Phi_{\theta}^{1} u_{\theta}) \cdot (\Phi_{\theta}^{1} \tau) = \int_{\Omega} (\alpha(\theta)^{-1} \Phi_{\theta}^{2} B_{\theta}) \cdot (\operatorname{curl} \Phi_{\theta}^{1} \tau),$$

$$\int_{\Omega} (\alpha(\theta) \Phi_{\theta}^{1} A_{\theta}) \cdot (\Phi_{\theta}^{1} v) = 0.$$

We therefore define

$$G : W^{1,\infty}(\mathbb{R}^3)^3 \times (\mathcal{Z}(\Omega) \times \nabla H^1(\Omega)) \to (\mathcal{Z}(\Omega) \times \nabla H^1(\Omega))'$$

by

$$G(\theta; \sigma, u)(\tau, v) = \int_{\Omega} (\alpha(\theta)^{-1} \operatorname{curl} \sigma) \cdot \operatorname{curl} \tau + \int_{\Omega} (\alpha(\theta)u) \cdot \tau + \int_{\Omega} (\alpha(\theta)\sigma) \cdot v - \int_{\Omega} (\alpha(\theta)^{-1} \Phi_{\theta}^2 B_{\theta}) \cdot \operatorname{curl} \tau,$$

so that $(A_{\theta}, u_{\theta})$ solves Eqs. (28) and (29) if and only if $G(\theta; \Phi_{\theta}^1 A_{\theta}, \Phi_{\theta}^1 u_{\theta}) = 0$.

By Lemma 3, we know that $\theta \mapsto \alpha(\theta)$ and $\theta \mapsto \alpha(\theta)^{-1}$ are $C^1$. Furthermore, we know from Proposition 6 that $\theta \mapsto \Phi_{\theta}^2 B_{\theta}$ is $C^1$. Therefore, by linearity and continuity of $G$ with respect to $(\sigma, u)$, we know that $G$ is $C^1$. We have

$$D_{\sigma,u}G(0; A_0, v_0)(\sigma', u')(\tau, v) = \int_{\Omega} \operatorname{curl} \sigma' \cdot \operatorname{curl} \tau + \int_{\Omega} u' \cdot \tau + \int_{\Omega} \sigma' \cdot v,$$

so that $D_{\sigma,u}G(0; A_0, v_0)$ is an isomorphism by the inf-sup conditions proven in [Val19, Section IV]. This proves, by the implicit function theorem, that for $\theta$ small enough there is a unique mapping $\theta \mapsto (\sigma(\theta), u(\theta))$ such that $G(\theta; \sigma(\theta), u(\theta)) = 0$. By uniqueness, we get $(\sigma(\theta), u(\theta)) = (\Phi_{\theta}^1 A_{\theta}, \Phi_{\theta}^1 u_{\theta})$. □
## 3.3 Proof of Theorem 1

Now that we have introduced ways to pullback functions and vector fields onto the fixed domain $\Omega$, and that we have derived differentiability of $\Phi_{\theta}^2 B_{\theta}$ and $\Phi_{\theta}^1 A_{\theta}$, we can conclude the proof of Theorem 1. To do so, we simply pullback the integral of $B_{\theta}$ against $A_{\theta}$ onto $\Omega$ using Lemma 1, and use the differentiability results from the last section.

*Proof.* We have

$$ 
\begin{align*}
\mathcal{H}(\Omega_{\theta}) &= \int_{\Omega_{\theta}} B_{\theta} \cdot A_{\theta} \\
&= \int_{\Omega} \Phi_{\theta}^{3}(B_{\theta} \cdot A_{\theta}) \\
&= \int_{\Omega} (\Phi_{\theta}^{2} B_{\theta}) \cdot (\Phi_{\theta}^{1} A_{\theta}).
\end{align*}
$$

From Propositions 6 and 7, we know that $\Phi_{\theta}^2 B_{\theta}$ and $\Phi_{\theta}^1 A_{\theta}$ are differentiable at zero, so that $\mathcal{H}$ is differentiable at $\Omega$ and, denoting their differentials in the direction $\theta$ as $B'$ and $A'$, we have

$$ \mathcal{H}'(\Omega; \theta) = \int_{\Omega} B' \cdot A + \int_{\Omega} B \cdot A'. $$

Since $B = \mathrm{curl} A$, and $A, A'$ are both in $\mathcal{Z}(\Omega)$, we know from [Alo+18, Lemma 7] that

$$ \int_{\Omega} B \cdot A' = \int_{\Omega} \mathrm{curl} A \cdot A' = \int_{\Omega} \mathrm{curl} A' \cdot A. $$

Furthermore, since $\mathrm{curl} A_{\theta} = B_{\theta}$, we get $\mathrm{curl} \Phi_{\theta}^1 A_{\theta} = \Phi_{\theta}^2 B_{\theta}$. By differentiating, we get $\mathrm{curl} A' = B'$, so that

$$ \mathcal{H}'(\Omega; \theta) = 2 \int_{\Omega} B' \cdot A. $$

Using the first equation in Eq. (27), we get

$$ \int_{\Omega} B' \cdot A = \int_{\Omega} (\alpha'(0; \theta) B) \cdot A. $$

Finally, since $B$ and $A$ are in $H(\mathrm{curl}, \Omega) \cap H_0(\mathrm{div}, \Omega)$, we can use Lemma 3 to get

$$ 
\begin{align*}
\int_{\Omega} (\alpha'(0; \theta) B') \cdot A &= \int_{\Omega} \mathrm{div} B \cdot A + \int_{\Omega} B \times \theta \cdot \mathrm{curl} A - \int_{\Omega} \mathrm{curl} B \times \theta \cdot A + \int_{\Omega} B \cdot \theta \mathrm{div} A + \int_{\partial\Omega} (B \cdot A)(\theta \cdot n) \\
&= \int_{\partial\Omega} (B \cdot A)(\theta \cdot n).
\end{align*}
$$

# 4 Approximation by finite element exterior calculus

In this section we assume that $\Omega$ is a toroidal domain with polyhedral boundary and that $(\mathcal{T}_h)_{h>0}$ is a quasi-uniform family of tetrahedron meshes of $\bar{\Omega}$, with $h$ the largest diameter of the cells. From [Amr+98, Proposition 3.7], we know that we can choose $s > 1/2$ so that $\Omega$ is $s$-regular. We denote by $\Delta^0(\mathcal{T}_h)$, $\Delta^1(\mathcal{T}_h)$, $\Delta^2(\mathcal{T}_h)$ and $\Delta^3(\mathcal{T}_h)$ the sets of points, edges, faces and cells of the mesh $\mathcal{T}_h$, respectively. For a $k$-simplex $S$ of the mesh $\mathcal{T}_h$, we define $\Delta^i(S)$ as the set of $i$-simplices which have non-empty intersection with $\bar{S}$ for $i < k$. Moreover, we assume that the cutting surface $\Sigma$ and toroidal curve $\gamma'$ introduced in Section 2.2 are compatible with the mesh, that is, are given by unions of elements of $\Delta^2(\mathcal{T}_h)$ and $\Delta^1(\mathcal{T}_h)$ respectively. Throughout this section, inequality constants denoted by $C$ are independent of $h$, but may depend on the domain $\Omega$.

Our goal is to provide a finite elements scheme to compute both the harmonic helicity of the domain $\Omega$ and its shape gradient. In Section 4.1, we recall the main properties of classical families of elements coming from finite elements exterior calculus. Then, we recall in Section 4.2 some notions related to the discrete version of the De Rham complex, and the related discrete harmonic fields. Finally, in Section 4.3 we prove the convergence of the numerical harmonic helicity.
## 4.1 Classical finite elements exterior calculus families

Here, we introduce some classical families in finite elements exterior calculus. The main idea is to define a discretized version of each functional space in the De Rham complex, and stable projections onto these spaces which commute with the differential operators. All these notions were first introduced in [RT77; Ned80], and later generalized for differential forms in [AFW06].

For $r$ an integer and $T$ a tetrahedral domain in $\mathbb{R}^3$, we define $\mathcal{P}_r(T)$ as the set of polynomials on $T$ with degree at most $r$, and $\tilde{\mathcal{P}}_r(T)$ the set of homogeneous polynomials on $T$ of degree $r$. Then, $\mathcal{P}_r^{-,\text{curl}}(T)$ and $\mathcal{P}_r^{-,\text{div}}$ are defined as

$$ \begin{aligned} \mathcal{P}_r^{-,\text{curl}}(T) &= \mathcal{P}_{r-1}(T)^3 \oplus \{p \in \tilde{\mathcal{P}}_r(T)^3 \mid p \cdot \mathbf{x} = 0\}, \\ \mathcal{P}_r^{-,\text{div}}(T) &= \mathcal{P}_{r-1}(T)^3 \oplus \{\mathbf{x}p \mid p \in \tilde{\mathcal{P}}_{r-1}(T)\}. \end{aligned} $$

For $r$ positive, we then define the discretizations

$$ \begin{aligned} V_h^0(\Omega) &= \{u \in H^1(\Omega) \mid u|_K \in \mathcal{P}_r(T) \ \forall T \in \mathcal{T}_h\}, \\ V_h^1(\Omega) &= \{u \in H(\text{curl}, \Omega) \mid u|_K \in \mathcal{P}_r^{-,\text{curl}}(T) \ \forall T \in \mathcal{T}_h\}, \\ V_h^2(\Omega) &= \{u \in H(\text{div}, \Omega) \mid u|_K \in \mathcal{P}_r^{-,\text{div}}(T) \ \forall T \in \mathcal{T}_h\}, \\ V_h^3(\Omega) &= \{u \in L^2(\Omega) \mid u|_K \in \mathcal{P}_{r-1}(T) \ \forall T \in \mathcal{T}_h\}. \end{aligned} $$

We then obtain the following sequence

$$ V_h^0(\Omega) \xrightarrow{\nabla} V_h^1(\Omega) \xrightarrow{\text{curl}} V_h^2(\Omega) \xrightarrow{\text{div}} V_h^3(\Omega). $$

These finite element spaces correspond respectively to the Lagrange $P^r$ elements, the Nedelec first kind elements of order $r$, the Raviart Thomas elements of order $r$, and the discontinuous $P^{r-1}$ elements.

One can check easily, using integration by parts, that $u$ is in $V_h^1(\Omega)$ if and only if the tangential trace of $u$ is continuous along all shared faces of $\Omega$. More precisely, if $T_1$ and $T_2$ are in $\Delta^3(\mathcal{T}_h)$ and $S_1, S_2$ are in $\Delta^2(T_1)$ and $\Delta^2(T_2)$ respectively with $T_1 \cap T_2 = S_1 = S_2$, we have

$$ u \times n_{S_1} + u \times n_{S_2} = 0. $$

Similarly, we have that $u$ is in $V_h^2(\Omega)$ if and only the normal trace of $u$ is continuous along all shared faces of $\Omega$.

We also introduce the discrete affine space $V_h^{\text{aff}}$ to solve Eq. (9) numerically. It is defined by

$$ V_h^{\text{aff}}(\Omega) = \{v \in H^1(\Omega \setminus \Sigma) \mid v|_T \in \mathcal{P}_r(T) \ \forall T \in \mathcal{T}_h, \text{ and } \llbracket v \rrbracket_{\Sigma} = 2\pi\} $$

where $\llbracket v \rrbracket_{\Sigma}$ denotes the jump of $v$ across $\Sigma$. The corresponding linear space is $V_h^0(\Omega)$.

We also use the smoothed quasi interpolators built in [AFW06, Section 5.4.] for differential forms, denoted by $\Pi_h^k$. It is known that these quasi interpolators are stable, and that they make the following diagram commute

$$ \begin{aligned} H^1(\Omega) &\xrightarrow{\nabla} H(\text{curl}, \Omega) \xrightarrow{\text{curl}} H(\text{div}, \Omega) \xrightarrow{\text{div}} L^2(\Omega) \\ \Pi_h^0 &\xrightarrow{\nabla} \Pi_h^1 &\xrightarrow{\text{curl}} \Pi_h^2 &\xrightarrow{\text{div}} \Pi_h^3. \\ V_h^0(\Omega) &\xrightarrow{\nabla} V_h^1(\Omega) \xrightarrow{\text{curl}} V_h^2(\Omega) \xrightarrow{\text{div}} V_h^3(\Omega) \end{aligned} $$

We then have the following approximation estimates on the quasi interpolators (see for example [EG18, Theorem 2.2. and 2.3.]).

**Proposition 8.** We have, for all $k$ between 0 and 3 and $u \in H^s$ regular,

$$ \|u - \Pi_h^k u\|_{L^2} \leq Ch^s \|u\|_{H^s}. $$
## 4.2 Discretization of the De Rham complex

In this section we state some results about the discretization of the De Rham complex. We borrow here most of our notations and lemmas from [AFW06]. We begin by defining discrete equivalents of closed and exact fields, which allow us to define the discrete harmonic fields. This then allows us to derive Hodge decompositions in the discrete setting, and uniform Poincaré inequalities. We then state some lemmas which will be useful for the coming convergence results.

First, we use some notations from the differential forms setting to unify some definitions and results. We denote

$$H\Lambda^0(\Omega) = H^1(\Omega), \quad H\Lambda^1(\Omega) = H(\mathrm{curl}, \Omega), \quad H\Lambda^2(\Omega) = H(\mathrm{div}, \Omega), \quad \text{and} \quad H\Lambda^3(\Omega) = L^2(\Omega),$$

as well as

$$d^0 = \nabla, \quad d^1 = \mathrm{curl}, \quad \text{and} \quad d^2 = \mathrm{div}.$$

We also denote $d^k = 0$ when $k$ is negative or larger than 3. The corresponding traceless spaces are denoted $\dot{H}\Lambda^k(\Omega)$, and the discrete traceless spaces are given by $\dot{V}_h^k(\Omega) = V_h^k(\Omega) \cap \dot{H}\Lambda^k(\Omega)$.

We are now able to define the discrete harmonic fields. We denote

$$\mathcal{B}_h^k(\Omega) = d^{k-1}V_h^{k-1}(\Omega), \quad \mathcal{Z}_h^k(\Omega) = \{u \in V_h^k(\Omega) \mid d^k u = 0\},$$

and

$$\mathcal{K}_h^k(\Omega) = \mathcal{Z}_h^k(\Omega) \cap \mathcal{B}_h^k(\Omega)^\perp.$$

Similarly, we define for traceless spaces

$$\mathring{\mathcal{B}}_h^k(\Omega) = d^{k-1}\mathring{V}_h^{k-1}(\Omega), \quad \mathring{\mathcal{Z}}_h^k(\Omega) = \{u \in \mathring{V}_h^k(\Omega) \mid d^k u = 0\},$$

and

$$\mathring{\mathcal{K}}_h^k(\Omega) = \mathring{\mathcal{Z}}_h^k(\Omega) \cap \mathring{\mathcal{B}}_h^k(\Omega)^\perp.$$

From the equivalence of discrete and continuous De Rham cohomology (see e.g. [AFW10, Section 5.6]), we get

$$\mathcal{K}_h^0(\Omega) \cong \mathbb{R}, \quad \mathcal{K}_h^1(\Omega) \cong \mathbb{R}, \quad \mathcal{K}_h^2(\Omega) \cong 0 \text{ and } \mathcal{K}_h^3(\Omega) \cong 0.$$

Similarly, from the equivalence of discrete and continuous De Rham cohomology with boundary condition, and Poincaré duality, we get

$$\mathring{\mathcal{K}}_h^0(\Omega) \cong 0, \quad \mathring{\mathcal{K}}_h^1(\Omega) \cong 0, \quad \mathring{\mathcal{K}}_h^2(\Omega) \cong \mathbb{R} \text{ and } \mathring{\mathcal{K}}_h^3(\Omega) \cong \mathbb{R}.$$

**Remark 2.** In the continuous case, we defined the harmonic fields as $\mathcal{K}(\Omega) = H(\mathrm{curl}^0, \Omega) \cap H_0(\mathrm{div}^0, \Omega)$. By using $H_0(\mathrm{div}^0, \Omega) = \nabla H^1(\Omega)^\perp$ and $H(\mathrm{curl}^0, \Omega) = \mathrm{curl} H_0(\mathrm{curl}, \Omega)^\perp$ from Proposition 17 of Appendix A, we recover two expressions for $\mathcal{K}(\Omega)$ which are similar to the discrete ones. However, in the discrete case, the two spaces $\mathcal{K}_h^1(\Omega)$ and $\mathring{\mathcal{K}}_h^2(\Omega)$ do not coincide. Indeed, elements of $\mathcal{K}_h^1(\Omega)$ are curl free but only weakly divergence free and tangent to the boundary meaning that they are not orthogonal to $\nabla H^1(\Omega)$ but only to the finite dimensional space $\nabla V_h^0(\Omega)$. On the other hand, elements of $\mathring{\mathcal{K}}_h^2(\Omega)$ are divergence free and tangent to the boundary, but only weakly curl free, that is, not orthogonal to $\mathrm{curl} H_0(\mathrm{curl}^0, \Omega)$ but to $\mathrm{curl} \mathring{V}_h^1(\Omega)$.

From these definitions, it is now straightforward to find the discrete Hodge decompositions.

$$V_h^k(\Omega) = \mathcal{Z}_h^k(\Omega)^\perp \oplus \mathcal{K}_h^k(\Omega) \oplus \mathcal{B}_h^k(\Omega),$$

$$\mathring{V}_h^k(\Omega) = \mathring{\mathcal{Z}}_h^k(\Omega)^\perp \oplus \mathring{\mathcal{K}}_h^k(\Omega) \oplus \mathring{\mathcal{B}}_h^k(\Omega).$$

There is also an equivalent of the Poincaré inequality in the discrete case given in [AFW06][Lemma 5.11].

**Proposition 9.** There exists a constant $C$, independent of $h$, such that for all $u$ in $\mathcal{Z}_h^k(\Omega)^\perp$

$$\|u\|_{L^2} \le C \|d^k u\|_{L^2}.$$

We will also be using the following lemma for the proof of convergence of the harmonic field. This lemma is proven in [AFW06][Lemma 5.9] for harmonic forms.

**Lemma 4.** For all $u_h$ in $\mathcal{K}_h^1(\Omega)$ there exists $u$ in $\mathcal{K}(\Omega)$ such that $\|u\| \le \|u_h\|$ and

$$\|u_h - u\|_{L^2} \le \|u - \Pi_h^1 u\|_{L^2}.$$
## 4.3 Numerical convergence of the harmonic helicity

We begin by studying the approximation of the harmonic field. As we have seen in Section 2.2, there are two different variational formulations for the harmonic fields. Although they give the same fields in the continuous case, this will not be true at the discrete level. As we will see, the classical Poisson formulation will give a discrete harmonic field in $K_h^1(\Omega)$, and the mixed formulation will give a discrete harmonic field in $\check{K}_h^2(\Omega)$.

Before studying the convergence of the numerical solutions, we state their well-posedness.

**Proposition 10.** There exists a unique solution to the following variational formulation. Find $u_h \in V_h^{\text{aff}}(\Omega)$ such that, for all $v_h \in V_h^0(\Omega)$,

$$ \int_{\Omega} u_h \cdot \nabla v_h = 0. \qquad (30) $$

Furthermore, $B_h^{\text{curl}} = \nabla u_h$ is in $K_h^1(\Omega)$.

*Proof.* We denote by $u$ the solution to Eq. (9), and $B = \nabla u$ the normalized harmonic field of $\Omega$. We begin by noticing that there is a unique solution to the following variational problem. Find $\tilde{u}_h \in V_h^0(\Omega)$ such that for all $v_h \in V_h^0(\Omega)$

$$ \int_{\Omega} \nabla \tilde{u}_h \cdot \nabla v_h = - \int_{\Omega} \Pi_h^1 B \cdot v_h. \qquad (31) $$

Indeed, this problem is a classical discretization of a Poisson equation, so the Lax-Milgam theorem applies. Now, define $u_h = \tilde{u}_h + \Pi_h^0 u$. Since $\nabla u_h = \nabla \tilde{u}_h + \Pi_h^1 B$, and $\Pi_h^0 u$ is in $V_h^{\text{aff}}(\Omega)$, we notice that $u_h$ is solution to Eq. (30) if and only if $\tilde{u}_h$ is solution to Eq. (31). As a consequence, Eq. (30) is also well-posed.

We now prove that $\nabla u_h$ is in $V_h^1(\Omega)$. It is straightforward that interelement continuity will be verified on all surfaces which are not included in $\Sigma_h$. Therefore, we take $S \in \Delta^2(\mathcal{T}_h)$ included in $\Sigma_h$. We know that there exist two cells $T_1$ and $T_2$ in $\Delta^3(\mathcal{T}_h)$ such that $S \in \Delta^2(T_1) \cap \Delta^2(T_2)$. We then denote $S = S_1$ when seen as an element of $\Delta^2(T_1)$, and $S = S_2$ when seen as an element of $\Delta^2(T_2)$, and order $T_1$ and $T_2$ so that the exterior normal of $T_1$ on $S_1$ is $n_{\Sigma_h}$. Since $u_h$ is in $V_h^{\text{aff}}(\Omega)$, we have $u_h|_{S_2} - u_h|_{S_1} = 2\pi$, so that $\nabla u_h \times n_{\Sigma_h}|_{S_2} = \nabla u_h \times n_{\Sigma_h}|_{S_1}$. As a consequence, interelement continuity is verified across $S$, so that $B_h^{\text{curl}} = \nabla u_h$ is in $V_h^1(\Omega)$.

Finally, we prove that $B_h^{\text{curl}}$ is in $K_h^1(\Omega)$. The fact that $B_h^{\text{curl}}$ is in $B_h^1(\Omega)^\perp$ is given directly by Eq. (30). Now, we take $T$ in $\Delta^3(\mathcal{T}_h)$. Since $B_h^{\text{curl}}|_T = \nabla u_h|_T$, we have $B_h^{\text{curl}}|_T \in B_h^1(T)$. Now, using that $B_h^1(T)$ is a subset of $Z_h^1(T)$, and that $B_h^{\text{curl}}$ is in $Z_h^1(\Omega)$ if and only if $B_h^{\text{curl}}|_T$ is in $Z_h^1(T)$ for all $T$, we obtain the desired result. $\square$

**Proposition 11.** There exists a unique solution to the following variational formulation. Find $(B_h^{\text{div}}, u_h) \in \check{V}_h^2(\Omega) \times (V_h^3(\Omega) \cap L_0^2(\Omega))$ such that, for all $(\tau_h, v_h) \in \check{V}_h^2(\Omega) \times V_h^3(\Omega)$

$$ \begin{gathered} \int_{\Omega} (\operatorname{div} B_h^{\text{div}}) v_h = 0, \\ \int_{\Omega} B_h^{\text{div}} \cdot \tau_h + \int_{\Omega} u_h (\operatorname{div} \tau_h) = 2\pi \int_{\Sigma} \tau_h \cdot n_{\Sigma}. \end{gathered} \qquad (32) $$

Furthermore, $B_h^{\text{div}}$ is in $\check{K}_h^2(\Omega)$.

*Proof.* The well-posedness comes from the exact same arguments as in the continuous case, by replacing the Poincaré inequalities and Hodge decompositions by their discrete counterparts.

We now prove that $B_h^{\text{div}}$ is in $\check{K}_h^2(\Omega)$. The fact that $B_h^{\text{div}}$ is in $\check{V}_h^2(\Omega)$ comes directly from the first equation of (32). Now, taking $\tau_h = \operatorname{curl} \rho_h$ in the second equation for $\rho_h \in \check{V}_h^1(\Omega)$, we get

$$ \begin{align*} \int_{\Omega} B_h^{\text{div}} \cdot \operatorname{curl} \rho_h &= 2\pi \int_{\Sigma} \operatorname{curl} \rho_h \cdot n_{\Sigma} \\ &= 2\pi \int_{\gamma'} \rho_h \cdot t' \\ &= 0, \end{align*} $$

so that $B_h^{\text{div}}$ is in $\check{B}_h^2(\Omega)^\perp$.

We prove the two following approximation results.
**Proposition 12.** There exists a constant C independent of h such that

$$ \begin{aligned} \|B_h^{\text{curl}} - B\|_{L^2} &\le Ch^s \|B\|_{L^2}, \\ \|B_h^{\text{div}} - B\|_{L^2} &\le Ch^s \|B\|_{L^2}. \end{aligned} $$

*Proof.* We begin by proving the convergence for $B_h^{\text{curl}}$. From Lemma 4, we know there exists $\tilde{B}_h$ in $K(\Omega)$ such that

$$ \|\tilde{B}_h - B_h^{\text{curl}}\|_{L^2} \le \|\tilde{B}_h - \Pi_h^1 \tilde{B}_h\|_{L^2}. $$

From the continuity of the circulation along $\gamma$ with respect to the Hcurl norm, we have

$$ \begin{aligned} \left| \int_{\gamma} (\tilde{B}_h - B) \cdot t \right| &= \left| \int_{\gamma} (\tilde{B}_h - B_h^{\text{curl}}) \cdot t \right| \\ &\le C \left\| \tilde{B}_h - B_h^{\text{curl}} \right\|_{H^{\text{curl}}} \\ &\le C \left\| \tilde{B}_h - B_h^{\text{curl}} \right\|_{L^2} \\ &\le C \left\| \tilde{B}_h - \Pi_h^1 \tilde{B}_h \right\|_{L^2}. \end{aligned} $$

Since $B$ and $\tilde{B}_h$ are in $K(\Omega)$, we get from [Val19, Lemma 5] that

$$ \begin{aligned} \|\tilde{B}_h - B\|_{L^2} &\le C \left( \|\text{div } \tilde{B}_h - \text{div } B\|_{L^2} + \|\text{curl } \tilde{B}_h - \text{curl } B\|_{L^2} + \left| \int_{\gamma} (\tilde{B}_h - B) \cdot t \right| \right) \\ &\le C \left\| \tilde{B}_h - \Pi_h^1 \tilde{B}_h \right\|_{L^2}. \end{aligned} $$

We then get by triangle inequality

$$ \begin{aligned} \|B_h^{\text{curl}} - B\|_{L^2} &\le \|B_h^{\text{curl}} - \tilde{B}_h\|_{L^2} + \|\tilde{B}_h - B\|_{L^2} \\ &\le C \left\| \tilde{B}_h - \Pi_h^1 \tilde{B}_h \right\|_{L^2}. \end{aligned} $$

Finally, using Proposition 8 and the continuous injection of $H(\text{curl}, \Omega) \cap H_0(\text{div}, \Omega)$ into $H^s(\Omega)$ for some $1/2 < s \le 1$, we get

$$ \begin{aligned} \|B_h^{\text{curl}} - B\|_{L^2} &\le Ch^s \|\tilde{B}_h\|_{H^s} \\ &\le Ch^s \|\tilde{B}_h\|_{L^2} \\ &\le Ch^s \|B\|_{L^2}. \end{aligned} $$

We now prove the convergence of $B_h^{\text{div}}$ to $B$. Using Eqs. (15) and (32), we have

$$ \begin{aligned} \int_{\Omega} B \cdot \tau + \int_{\Omega} u (\text{div } \tau_h) &= \int_{\Sigma} \tau_h \cdot n_{\Sigma}, \\ \int_{\Omega} B_h^{\text{div}} \cdot \tau + \int_{\Omega} u_h (\text{div } \tau_h) &= \int_{\Sigma} \tau_h \cdot n_{\Sigma}, \end{aligned} $$

for all $\tau_h$ in $\tilde{V}_h^2(\Omega)$. As a consequence, for all $\tau_h$ in $\tilde{Z}_h^2(\Omega)$, we get

$$ \int_{\Omega} (B_h^{\text{div}} - B) \cdot \tau_h = 0. $$

Since both $B_h^{\text{div}}$ and $\Pi_h^2 B$ are in $\tilde{Z}_h^2(\Omega)$, we have

$$ \int_{\Omega} (B_h^{\text{div}} - B) \cdot (B_h^{\text{div}} - B) = \int_{\Omega} (B_h^{\text{div}} - B) \cdot (\Pi_h^2 B - B), $$

so that

$$ \|B_h^{\text{div}} - B\|_{L^2} \le \|\Pi_h^2 B - B\|_{L^2}. $$

Using again Proposition 8 and the fact that $\Omega$ is s-regular, we get the desired result

$$ \|B_h^{\text{div}} - B\|_{L^2} \le Ch^s \|B\|_{L^2}. $$

□
We now study the well-posedness and the convergence of the vector potentials of B. A first remark is that the vector potential orthogonal to B is well studied in the literature via the Hodge Laplacian (see for example [AFW06]). As such, we study the approximation of such vector potentials which we denote $A_h^1$, and recover $A_h^2$, the one which is circulation free, by the discrete counterpart of Eq. (20)

$$A_h^2 = A_h^1 - \frac{1}{2\pi} \left( \int_{\gamma'} A_h^1 \cdot t' \right) B_h^{\text{curl}}. \qquad (33)$$

**Proposition 13.** There exists a unique solution to the following variational formulation. Find $(A_h^1, u_h) \in V_h^1(\Omega) \times V_h^2(\Omega)$ such that, for all $(\tau_h, v_h) \in V_h^1(\Omega) \times V_h^2(\Omega)$,

$$\int_{\Omega} A_h^1 \cdot \tau_h = \int_{\Omega} \mathrm{curl} \tau_h \cdot u_h, \\ \int_{\Omega} \mathrm{curl} A_h^1 \cdot v_h + \int_{\Omega} (\mathrm{div} u_h) (\mathrm{div} v_h) = \int_{\Omega} B_h^{\mathrm{div}} \cdot v_h. \qquad (34)$$

Furthermore, we have $\mathrm{curl} A_h^1 = B_h^{\mathrm{div}}$.

*Proof.* The well posedness is proven in the exact same way as in the continuous case Proposition 2, by replacing the Hodge decomposition and Poincaré inequality by their discrete counterparts. The fact that $\mathrm{curl} A_h^1 = B_h^{\mathrm{div}}$ is also proven in a similar way as in Proposition 2. Using $K_h^2(\Omega) = 0$, and the discrete Hodge decomposition (see Section 4.2), we have

$$V_h^2(\Omega) = Z_h^2(\Omega)^{\perp} \oplus B_h^2(\Omega),$$

so that $u_h = u_h^{\nabla} + u_h^{\mathrm{curl}}$ with $\int_{\Omega} u_h^{\nabla} v = 0$ for all $v_h$ in $V_h^2(\Omega)$ with $\mathrm{div} v_h = 0$, and $\mathrm{div} u_h^{\mathrm{curl}} = 0$. Therefore, taking $v_h = u_h^{\mathrm{curl}}$ in the second equation of Eq. (34), we get

$$\int_{\Omega} (\mathrm{div} u_h) (\mathrm{div} u_h^{\nabla}) = 0,$$

so that, since $\mathrm{div} u_h^{\mathrm{curl}} = 0$, we obtain $\mathrm{div} u_h^{\nabla} = \mathrm{div} u_h = 0$. From this, we get $\int_{\Omega} A_h^1 \cdot v_h = \int_{\Omega} B_h^{\mathrm{div}} \cdot v_h$ for all $v_h \in V_h^2(\Omega)$, and $B_h^{\mathrm{div}} \in V_h^2(\Omega)$ implies the desired equality. $\square$

**Proposition 14.** There exists a constant C independent of h such that

$$\|A_h^1 - A^1\|_{L^2} \leq C h^s \|B\|_{L^2}.$$

*Proof.* Define $(\tilde{A}_h^1, \tilde{u}_h)$ as the solution to the following variational problem. Find $(\tilde{A}_h^1, \tilde{u}_h) \in V_h^1(\Omega) \times V_h^2(\Omega)$ such that, for all $(\tau_h, v_h) \in V_h^1(\Omega) \times V_h^2(\Omega)$,

$$\int_{\Omega} \tilde{A}_h^1 \cdot \tau_h = \int_{\Omega} \mathrm{curl} \tau_h \cdot \tilde{u}_h, \\ \int_{\Omega} \mathrm{curl} \tilde{A}_h^1 \cdot v_h + \int_{\Omega} (\mathrm{div} \tilde{u}_h) (\mathrm{div} v_h) = \int_{\Omega} B \cdot v_h.$$

The well-posedness of this equation is obtained in the same way as for $A_h^1$. Furthermore, by continuity of the resolvant, and independence of the discrete Poincaré inequality constants in h, we get

$$\| \tilde{A}_h^1 - A_h^1 \|_{L^2} \leq C \| B - B_h^{\mathrm{div}} \|_{L^2}.$$

Furthermore, we obtain from of [AFW06, Theorem 7.9] and $\|B\|_{H^s} \leq C \|B\|_{L^2}$ that

$$\| \tilde{A}_h^1 - A^1 \|_{L^2} \leq C h^s \| B \|_{L^2}.$$

We then get the desired result from a triangle inequality and Proposition 12. $\square$

As we have seen in Section 2.3, there are two ways of computing the harmonic helicity of $\Omega$. The first one is done by computing circulations of $A^1$, and the second by taking the $L^2$ product of B and $A^2$. As we shall see in the following results, we can also recover the numerical harmonic helicity in two similar ways. To do this, we first need to establish convergence of the circulation of $A_h^1$, which is given by the following lemma.
**Lemma 5.** We have

$$ \left| \int_{\gamma'} (A_h^1 - A^1) \cdot t' \right| \leq Ch^s \|B\|_{L^2}. $$

*Proof.* From Proposition 13, we have

$$ \begin{aligned} \|A_h^1 - A^1\|_{H_{\text{curl}}}^2 &= \|A_h^1 - A^1\|_{L^2}^2 + \|\text{curl } A_h^1 - \text{curl } A^1\|_{L^2}^2 \\ &= \|A_h^1 - A^1\|_{L^2}^2 + \|B_h^{\text{div}} - B\|_{L^2}^2, \end{aligned} $$

so that by Propositions 12 and 14, we get

$$ \|A_h^1 - A^1\|_{H_{\text{curl}}} \leq Ch^s \|B\|_{L^2}. $$

Now, since curl $A_h^1 = B_h^{\text{div}}$, which is in $H_0(\text{div}, \Omega)$, the circulation of $A_h^1$ along $\gamma'$ is well-defined by Eq. (13), and this circulation is bounded by the $H(\text{curl}, \Omega)$ norm. We therefore have

$$ \left| \int_{\gamma'} (A_h^1 - A^1) \cdot t \right| \leq Ch^s \|B\|_{L^2}. $$

**Corollary 2.** Defining $A_h^2$ as in Eq. (33), we get

$$ \|A_h^2 - A^2\|_{L^2} \leq Ch^s \|B\|_{L^2}. $$

*Proof.* We get this result by writing

$$ A_h^2 - A^2 = A_h^1 - A^1 - \frac{1}{2\pi} \left( \int_{\gamma'} (A_h^1 - A^1) \cdot t \right) B_h^{\text{curl}} + \frac{1}{2\pi} \left( \int_{\gamma'} A^1 \cdot t \right) (B - B_h^{\text{curl}}). $$

We then obtain the desired result using Propositions 12 and 14 and Lemma 5.

**Theorem 2** (Convergence of harmonic helicity). We have

$$ |\mathcal{H}(\Omega) - \mathcal{H}(B_h^{\text{div}})| \leq Ch^s \|B\|_{L^2}^2, $$

where $\mathcal{H}(B_h^{\text{div}})$ can be either computed as

$$ \mathcal{H}(B_h^{\text{div}}) = - \left( \int_{\gamma'} A_h^1 \cdot t' \right) \left( \int_{\gamma} A_h^1 \cdot t \right), $$

or as

$$ \mathcal{H}(B_h^{\text{div}}) = \int_{\Omega} B_h^{\text{div}} \cdot A_h^2. $$

*Proof.* We use the formulation of the helicity through $A_h^2$, the other one being a simple consequence of the Bevir-Gray formula. To obtain the desired estimate, we write

$$ \int_{\Omega} B_h^{\text{div}} \cdot A_h^2 - \mathcal{H}(\Omega) = \int_{\Omega} (B_h^{\text{div}} - B) \cdot A_h^2 + \int_{\Omega} B \cdot (A_h^2 - A^2). $$

The conclusion then follows from Proposition 12 and Corollary 2.

**Remark 3.** Although we have proven convergence of the harmonic helicity, convergence of the shape gradient is a harder task. Indeed, such a result is related to traces of the numerical solutions. From the estimates obtained in this section, we only have convergence of $B_h^{\text{curl}}$ and $A_h^2$ in $H(\text{curl}, \Omega)$, and of $B_h^{\text{div}}$ in $H(\text{div}, \Omega)$. Hence, the tangential traces of both $A_h^2$ and $B_h^{\text{curl}}$ converges in $H^{-1/2}(\partial\Omega)$. Nevertheless, as we do not expect the numerical solutions to have any higher regularity than the ones prescribed by the discretized functional spaces$^4$, this is not sufficient to ensure convergence toward the shape gradient Eq. (23). This means that we cannot simply use Sobolev estimates inside $\Omega$, and need to work on the boundary directly in order to obtain, for example, convergence in $L^2(\partial\Omega)^3$ (or at least on the tangential part). Such a convergence result is therefore non-trivial with usual techniques, and this problem remains open.

<footnote>$^4$In the continuous setting, the reasoning works because $A^2$ and $B$ are in $H(\text{curl}, \Omega) \cap H_0(\text{div}, \Omega)$. Nevertheless, the divergence of $A_h^2$ and $B_h^{\text{curl}}$ is not in $L^2(\Omega)$ for Nedelec elements</footnote>
# 5 Numerical implementation and results

## 5.1 Specificity of simulations for stellarators

For both historical and practical considerations, it is frequently advantageous for the surfaces under examination to exhibit specific symmetries. In particular, a majority of stellarators are invariant under discrete rotations along the Oz axis, with an angle of $2\pi/N_p$, where $N_p$ takes values of 3 (as in the case of NCSX [Zar+01]) or 5 (as observed in W7X [War+17]). Additionally, stellarators commonly exhibit invariance under the continuous symmetry known as stellarator symmetry, as discussed in detail in [IPW20, Section 12.3] and [DH98]. In practice, these surfaces are represented by a set of coefficients $(R_{m,n}, (Z_{m,n})$ for $m \in \mathbb{N}$ and $n \in \mathbb{Z}$ which define the functions

$$R(u, v) = \sum_{m \ge 0} \sum_{n \in \mathbb{Z}} R_{m,n} \cos(2\pi(mu + nv)), \quad (35)$$

$$Z(u, v) = \sum_{m \ge 0} \sum_{n \in \mathbb{Z}} Z_{m,n} \sin(2\pi(mu + nv)). \quad (36)$$

Note the absence of sin terms for $R$ and cos terms for $Z$ due to the imposition of stellarator symmetry. The surface is then parametrized in cylindrical coordinates by $(R(u, v), \frac{2\pi v}{N_p}, Z(u, v))$ where $N_p$ stands for the discrete symmetry imposed.

For the numerical simulation, we truncate the number of Fourier components in (35) and (36).

## 5.2 Implementation

Using the symmetries of the system, we only work with one section of the stellarator, characterized by the coefficients $(R_{m,n}, Z_{m,n})$. Our first step is to use Gmsh [GR09] to mesh the interior of a polyhedral approximation of the surface.

Then we use the finite element library FENiCSx ([Scr+22b; Scr+22a; Aln+14]) to assemble the finite element problems Eqs. (30), (32) and (34) to obtain $B_h^{\text{div}}$, $B_h^{\text{curl}}$ and $A_h^1$. We implemented both first or second order FEEC elements defined in Section 4.1 with adequate periodic conditions at the cuts of our section. To solve the linear system required to get $B_h^{\text{curl}}$, that is a Poisson equation, we use the solver MUMPS: MUltifrontal Massively Parallel sparse direct Solver [Ame+01; Ame+06]. $B_h^{\text{div}}$ and $A_h$ are more expensive and complex to solve as both are defined in mixed formulations. Both are solved using the iterative Krylov method GMRES [SS86] applied after a careful preconditioning. More precisely, we use the block diagonal preconditioner proposed by Arnold et al. in [AFW06, Section 10.2] for Hodge-Laplacian problems. The preconditioning problems are solved using MUMPS. We believe that using an Auxiliary-space Maxwell (AMS) Solver [HX07] instead<fnref n="5" /> of a direct solver for the preconditioner would improve the efficiency and scalability of our method when using more than a few million degrees of freedom.

## 5.3 Numerical tests

Once we have computed $B_h^{\text{div}}$, $B_h^{\text{curl}}$ and $A_h$ that are represented in Fig. 2, we compute the harmonic helicity and its shape gradient. As was stated in Corollary 1, we know that axisymmetric tori have zero harmonic helicity. Numerically, we have found in this case a value of $10^{-9}$. This was computed with a major and minor radius of 1 and 0.1 respectively, second order elements, and $h$ equal to 0.025. For more interesting shapes, as NCSX's plasma, we do not have a reference. As mentioned in Section 5.1, we have a continuous description of the shape (i.e. not a polyhedral one), thus we perform two tests.

1. On the left side of Fig. 3, we perform a better and better approximation of the continuous shape. To this aim we provide finer and finer grid description of the surface to the mesher. Hence, variation of the obtained magnetic helicity are due both to variations of the domain and to the finite element approximation error.

2. On the right side of Fig. 3, we fix at the beginning a polyhedral shape that is an approximation of the continuous surface described by the set of coefficients $(R_{m,n}, Z_{m,n})$ and use finer and finer meshes of this fixed polyhedral shape as $h$ goes to 0. Hence, variations of the magnetic helicity are only due to the numerical approximation of the fields. This is the situation described by the theoretical analysis in Section 4. The reference solution was obtained on a high-performance computing cluster using a mesh

<footnote n="5">⁵We had issues porting AMS hypre [FY02] solvers to the background sparse linear library Petsc [Bal+98]</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2: One section (i.e. one third) of NCSX plasma. In the upper plot, we represent the function $u_h$ of Eq. (30). Its gradient $B_h^{\text{curl}}$ is shown in the middle figure. Note that the boundary conditions on the left and right cuts are a jump $2\pi/3$, whereas we have Neumann boundary condition on the plasma surface. The bottom figure is a representation of $A_h^2$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 3: On the left, we provide better and better approximations of the surface to the mesher. We then plot the characteristic size of an element versus the harmonic helicity. On the right, we mesh with different characteristic size $h$ the same polyhedral approximation of the continuous surface and compare it to a reference solution obtained with second order elements and $h = 0.025$. The two lines are merely indicative and their slopes correspond to theoretical first order (blue) and second order (red) convergence rates.</figcaption></figure>

size of $h = 0.025$ with second-order finite elements. This resulted in linear systems with approximately 5 million degrees of freedom, which were solved in a matter of minutes.

Except for the reference solution, all simulations can be conveniently performed on a laptop with 32 GB of RAM.

The shape gradient is shown in Fig. 4. As mentioned in Remark 3, we have not proven that this numerical approximation converges toward the correct shape gradient. However, we performed numerical tests using finite differences on the Fourier coefficients of the surface, which were consistent with the numerical shape gradient computed with $A_h^2$ and $B_h^{\text{curl}}$.

## 5.4 Two optimization programs

As the harmonic helicity of NCSX plasma is a negative quantity$^6$, we are interested in minimizing it. Indeed, we recall that in a stellarator, one generates a magnetic field that is in first approximation a harmonic field in the plasma domain $\Omega$. As the stability of the plasma is related to the twisting of the magnetic field, we are trying to increase the absolute value of the harmonic helicity to search for interesting and hopefully more stable plasma configurations.

As proved in Proposition 4, the harmonic helicity scales as $H(\lambda\Omega) = \lambda H(\Omega)$ for $\lambda > 0$. For this reason we impose an upper bound either on the volume of $\Omega$ or on the area of $\partial\Omega$. Under one of these constraints, it is unclear whether an optimal shape exists. Hence, we add the condition that $\partial\Omega$ satisfies a uniform ball condition. This imposes a strong regularity constraint on the surface and ensures existence of an optimal shape among these regular shapes, we refer for example to [PRS24; Ger23a]. Numerically, we set the uniform ball condition by imposing a lower limit on the minimal curvature radius everywhere on the surface $\partial\Omega$. These three constraints, namely upper bound on the volume $\Omega$ or on the area of $\partial\Omega$ and lower bound on the minimal curvature radii on $\partial\Omega$, along with their corresponding shape gradients, can be efficiently computed using the smooth parametrization of the surface outlined in Eqs. (35) and (36) and performing quadrature across either $\partial\Omega$ or $\Omega$.

Subsequently, we introduce smooth non-linear costs that blow up to guarantee the fulfillment of the specified constraints. Once all the costs are assembled, we apply the Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization algorithm from the scipy library [Vir+20] on the Fourier components of the description of the surface.

<footnote>$^6$Using a planar symmetry, we could also choose to take a positive value, see Proposition 4.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 4: One section of NCSX plasma. The plot shows the shape gradient $2A_h^2 \cdot B_h^{\text{curl}}$ on the boundary. NCSX plasma configuration has a negative harmonic helicity. Improving our criterion implies following the opposite of the shape gradient.</figcaption></figure>

<table><caption>Table 1: harmonic helicity and geometric properties of the optimized shapes in Fig. 5. As we aim to compare the optimized shapes with the NCSX reference, for BPC we penalized the area of the surface $\partial\Omega$ above $25 m^2$, whereas for BVC we penalized the volume of $\Omega$ above $3 m^3$. In both cases, we targeted a minimal curvature radius of $2 cm$; thus a more regular shape than the reference NCSX plasma shape.</caption><thead><tr><th>Plasma shape</th><th>Perimeter (m²)</th><th>Volume (m³)</th><th>min curvature radius (cm)</th><th>Harmonic helicity</th></tr></thead><tbody><tr><td>NCSX (ref)</td><td>24.52</td><td>2.96</td><td>0.83</td><td>-0.0220</td></tr><tr><td>BPC</td><td>25.02</td><td>5.45</td><td>1.89</td><td>-0.6427</td></tr><tr><td>BVC</td><td>33.5</td><td>2.99</td><td>2.00</td><td>-0.0726</td></tr></tbody></table>

**Bounded Perimeter and Curvature (BPC)** We initially focus on optimization with a bounded perimeter, employing a soft threshold set at $25 m^2$ to ensure comparability with the NCSX plasma shape. Our goal also involves achieving a minimal inverse curvature radius of $50 m^{-1}$, surpassing the regularity of the NCSX reference. The results of this simulation are depicted in Fig. 5 and summarized in Table 1. Notably, we observe a remarkable 30-fold improvement in harmonic helicity. While the shape appears to collapse toward the Oz axis, this effect is restrained by the imposed curvature constraints. This behavior may be attributed to our normalization choice for the harmonic field; as proximity to the Oz axis increases, the magnitude rises, given our fixed toroidal circulation.

**Bounded Volume and Curvature (BVC)** We also perform the optimization using a constraint on the volume instead of the perimeter. We apply an upper bound at $3 m^3$ along with the same curvature constraints. This time, we achieve more than a 3-fold improvement, and the shape does not seem to collapse toward the Oz axis. The shape is more intricate than the initial plasma shape of NCSX.

# 6 Conclusion and perspectives

We introduced a new shape functional for toroidal domains capturing the linkage of the corresponding harmonic field. Then, using carefull differential form pullbacks, we have been able to compute its shape derivative. We also showed that this functional can be efficiently computed numerically. We illustrated this using finite elements exterior calculus and applied it to a stellarator device for magnetic confinement in nuclear fusion.

It seems to us that the pure optimization of harmonic helicity provides degenerate forms, which appear to be less usable for applications in magnetic confinement. However, given the very significant potential gain (factors of 3 and 30 in our two simulations), studying multi-objective optimization with more conventional costs could yield interesting results.
<figure><img src="image_5.png" /><figcaption>Figure 5: The results of the two optimization problems, BPC (shown in green above) and BVC (shown in blue below). The associated costs are reported in Table 1. Additionally, the reference NCSX plasma is plotted in transparent orange.</figcaption></figure>
*Proof.* This is a simple consequence of the Hodge decomposition given in [AFW06, Section 2], in the cases $k=1$ and $k=2$ respectively. The reason no harmonic term appears in Eq. (39) is that the second De Rham cohomology space vanishes in $\Omega$, and thus, that the set of harmonic two forms is trivial. $\square$

**Proposition 16.** There exists $C$ such that for all $V_1 \in H(\mathrm{curl}^0, \Omega)^\perp \cap H(\mathrm{curl}, \Omega)$, $V_2 \in H(\mathrm{div}^0, \Omega)^\perp \cap H(\mathrm{div}, \Omega)$

$$ \|V_1\| \le C \|\mathrm{curl} V_1\|, $$

$$ \|V_2\| \le C \|\mathrm{div} V_2\|. $$

*Proof.* This is given by [AFW06, Eq. (2.17)] in the cases $k=1$ and $k=2$. $\square$

**Proposition 17.** We have the following orthogonality relations in $L^2(\Omega)^3$

$$ \nabla H^1(\Omega)^\perp = H_0(\mathrm{div}^0, \Omega), \qquad (40) $$

$$ \mathrm{curl} H(\mathrm{curl}, \Omega)^\perp = H_0(\mathrm{curl}^0, \Omega), \qquad (41) $$

$$ \nabla H_0^1(\Omega)^\perp = H(\mathrm{div}^0, \Omega), \qquad (42) $$

$$ \mathrm{curl} H_0(\mathrm{curl}, \Omega)^\perp = H(\mathrm{curl}^0, \Omega). \qquad (43) $$

Furthermore, all these linear subspaces are $L^2$ closed, so that the relations are still correct by taking the orthogonal to the other side.

*Proof.* Eqs. (40) to (43) are given by [AFW06, Eqs. (2.15) and (2.16)] in the cases $k=1$ and $k=2$. The fact that all the subspaces are $L^2$ closed is given by the continuity of $\mathrm{div}: H(\mathrm{div}, \Omega) \to L^2(\Omega)$ and $\mathrm{curl}: H(\mathrm{curl}, \Omega) \to L^2(\Omega)^3$ on the right-hand sides of Eqs. (40) to (43), and by [AFW06, Th. 2.3.] for the left-hand side. $\square$

## References

[Aln+14] M. S. Alnæs et al. “Unified form language: A domain-specific language for weak formulations of partial differential equations”. In: *ACM Transactions on Mathematical Software* 40.2 (2014), 9:1–9:37. DOI: 10.1145/2566630.

[Alo+18] A. Alonso-Rodríguez et al. “Finite Element Approximation of the Spectrum of the Curl Operator in a Multiply Connected Domain”. In: *Foundations of Computational Mathematics* 18.6 (2018), pp. 1493–1533. DOI: 10.1007/s10208-018-9373-4.

[Ame+01] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L’Excellent. “A fully asynchronous multifrontal solver using distributed dynamic scheduling”. In: *SIAM Journal on Matrix Analysis and Applications* 23.1 (2001), pp. 15–41.

[Ame+06] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, and S. Pralet. “Hybrid scheduling for the parallel solution of linear systems”. In: *Parallel Computing* 32.2 (2006), pp. 136–156.

[Amr+98] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault. “Vector potential in three-dimensional nonsmooth domains”. In: *Mathematical Methods in The Applied Sciences - MATH METH APPL SCI* 21 (1998), pp. 823–864. DOI: 10.1002/(SICI)1099-1476(199806)21:93.0.C0;2-B.

[AFW10] D. Arnold, R. Falk, and R. Winther. “Finite element exterior calculus: from Hodge theory to numerical stability”. In: *Bulletin of the American Mathematical Society* 47.2 (2010), pp. 281–354. DOI: 10.1090/S0273-0979-10-01278-4.

[AFW06] D. N. Arnold, R. S. Falk, and R. Winther. “Finite element exterior calculus, homological techniques, and applications”. In: *Acta Numerica* 15 (2006), pp. 1–155. DOI: 10.1017/S0962492906210018.

[Arn66] V. I. Arnold. “On the topology of three-dimensional steady flows of an ideal fluid”. In: *Journal of Applied Mathematics and Mechanics* 30.1 (1966), pp. 223–226. DOI: 10.1016/0021-8928(66)90070-0.

[Arn14] V. I. Arnold. “The asymptotic Hopf invariant and its applications”. In: *Vladimir I. Arnold - collected works: Hydrodynamics, bifurcation theory, and algebraic geometry 1965-1972*. Ed. by A. B. Givental et al. Berlin, Heidelberg: Springer Berlin Heidelberg, 2014, pp. 357–375. DOI: 10.1007/978-3-642-31031-7_32.
[AK21] V. I. Arnold and B. A. Khesin. *Topological Methods in Hydrodynamics*. Vol. 125. Applied Mathematical Sciences. Cham: Springer International Publishing, 2021. DOI: 10.1007/978-3-030-74278-2.

[Bal+98] S. Balay, W. Gropp, L. C. McInnes, and B. F. Smith. "PETSc, the portable, extensible toolkit for scientific computation". In: *Argonne National Laboratory* 2.17 (1998).

[BG80] M. K. Bevir and J. W. Gray. "Relaxation, flux consumption and quasi steady state pinches". In: *Proceedings of the Reversed-Field. Pinch Theory Workshop. Los Alamos, 1980*.

[BCS02] A. Buffa, M. Costabel, and D. Sheen. "On traces for $H(\text{curl}, \Omega)$ in Lipschitz domains". In: *Journal of Mathematical Analysis and Applications* 276.2 (2002), pp. 845-867. DOI: 10.1016/S0022-247X(02)00455-9.

[CDG01] J. Cantarella, D. DeTurck, and H. Gluck. "The Biot-Savart operator for application to knot theory, fluid dynamics, and plasma physics". In: *Journal of Mathematical Physics* 42.2 (2001). DOI: 10.1063/1.1329659.

[Can+00] J. Cantarella, D. DeTurck, H. Gluck, and M. Teytel. "Isoperimetric problems for the helicity of vector fields and the Biot-Savart and curl operators". In: *Journal of Mathematical Physics* 41.8 (2000), pp. 5615-5641. DOI: 10.1063/1.533429.

[Can+99] J. Cantarella, D. Deturck, H. Gluck, and M. Teytel. "Influence of Geometry and Topology on Helicity". In: *Magnetic Helicity in Space and Laboratory Plasmas*. American Geophysical Union (AGU), 1999, pp. 17-24. DOI: 10.1029/GM111p0017.

[DH98] R. Dewar and S. Hudson. "Stellarator symmetry". In: *Physica D: Nonlinear Phenomena* 112.1-2 (1998), pp. 275-280. DOI: 10.1016/S0167-2789(97)00216-9.

[EP23] A. Enciso and D. Peralta-Salas. "Non-existence of axisymmetric optimal domains with smooth boundary for the first curl eigenvalue". In: *ANNALI SCUOLA NORMALE SUPERIORE - CLASSE DI SCIENZE* (2023), pp. 311-327. DOI: 10.2422/2036-2145.202010_008.

[EG18] A. Ern and J.-L. Guermond. "Analysis of the edge finite element approximation of the Maxwell equations with low regularity solutions". In: *Computers & Mathematics with Applications* 75.3 (2018), pp. 918-932. DOI: 10.1016/j.camwa.2017.10.017.

[FY02] R. D. Falgout and U. M. Yang. "hypre: A Library of High Performance Preconditioners". In: *Computational Science — ICCS 2002*. Ed. by P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, and J. J. Dongarra. Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 2002, pp. 632-641. DOI: 10.1007/3-540-47789-6_66.

[Ger23a] W. Gerner. *Existence of optimal domains for the helicity maximisation problem among domains satisfying a uniform ball condition*. 2023. DOI: 10.48550/arXiv.2305.13642. arXiv: 2305.13642 [math-ph]. URL: http://arxiv.org/abs/2305.13642. preprint.

[Ger23b] W. Gerner. "Isoperimetric problem for the first curl eigenvalue". In: *Journal of Mathematical Analysis and Applications* 519.2 (2023), p. 126808. DOI: 10.1016/j.jmaa.2022.126808.

[GR09] C. Geuzaine and J.-F. Remacle. "Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities". In: *International Journal for Numerical Methods in Engineering* 79.11 (2009), pp. 1309-1331. DOI: 10.1002/nme.2579.

[HP18] A. Henrot and M. Pierre. *Shape variation and optimization: a geometrical analysis*. Vol. 28. EMS Tracts in Mathematics. Zürich, Switzerland: European Math. Soc (EMS), 2018.

[HX07] R. Hiptmair and J. Xu. "Nodal Auxiliary Space Preconditioning in $H(\text{curl})$ and $H(\text{div})$ Spaces". In: *SIAM Journal on Numerical Analysis* 45.6 (2007), pp. 2483-2509. DOI: 10.1137/060660588.

[IPW20] L.-M. Imbert-Gerard, E. J. Paul, and A. M. Wright. *An Introduction to Stellarators: From magnetic fields to symmetries and optimization*. arXiv:1908.05360 [physics]. 2020. DOI: 10.48550/arXiv.1908.05360.

[LRV15] E. Lara, R. Rodríguez, and P. Venegas. "Spectral approximation of the curl operator in multiply connected domains". In: *Discrete and Continuous Dynamical Systems - S* 9.1 (Mon Nov 30 19:00:00 EST 2015), pp. 235-253. DOI: 10.3934/dcdss.2016.9.235.

[MV19] D. MacTaggart and A. Valli. "Magnetic helicity in multiply connected domains". In: *Journal of Plasma Physics* 85.5 (2019). Publisher: Cambridge University Press, p. 775850501. DOI: 10.1017/S0022377819000576.
[Mof69] H. K. Moffatt. "The degree of knottedness of tangled vortex lines". In: *Journal of Fluid Mechanics* 35.1 (1969), pp. 117-129. DOI: 10.1017/S0022112069000991.

[Mon23] S. Montiel. *The Isoperimetric Problem for the Curl Operator*. 2023. arXiv: 2307.09556 [math]. URL: http://arxiv.org/abs/2307.09556. preprint.

[Ned80] J. C. Nedelec. "Mixed finite elements in R3". In: *Numerische Mathematik* 35.3 (1980), pp. 315-341. DOI: 10.1007/BF01396415.

[Pau+18] E. J. Paul, M. Landreman, A. Bader, and W. Dorland. "An adjoint method for gradient-based optimization of stellarator coil shapes". In: *Nuclear Fusion* 58.7 (2018), p. 076015. DOI: 10.1088/1741-4326/aac1c7.

[PRS22] Y. Privat, R. Robin, and M. Sigalotti. "Optimal shape of stellarators for magnetic confinement fusion". In: *Journal de Mathématiques Pures et Appliquées* 163 (2022), pp. 231-264. DOI: 10.1016/j.matpur.2022.05.005.

[PRS24] Y. Privat, R. Robin, and M. Sigalotti. "Existence of surfaces optimizing geometric and PDE shape functionals under reach constraint". In: *Interfaces and Free Boundaries* (2024). DOI: 10.4171/ifb/523.

[RT77] P. A. Raviart and J. M. Thomas. "A mixed finite element method for 2-nd order elliptic problems". In: *Mathematical Aspects of Finite Element Methods*. Ed. by I. Galligani and E. Magenes. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer, 1977, pp. 292-315. DOI: 10.1007/BFb0064470.

[SS86] Y. Saad and M. H. Schultz. "GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems". In: *SIAM Journal on Scientific and Statistical Computing* 7.3 (1986), pp. 856-869. DOI: 10.1137/0907058. eprint: https://doi.org/10.1137/0907058.

[Scr+22a] M. W. Scroggs, I. A. Baratta, C. N. Richardson, and G. N. Wells. "Basix: a runtime finite element basis evaluation library". In: *Journal of Open Source Software* 7.73 (2022), p. 3982. DOI: 10.21105/joss.03982.

[Scr+22b] M. W. Scroggs, J. S. Dokken, C. N. Richardson, and G. N. Wells. "Construction of Arbitrary Order Finite Element Degree-of-Freedom Maps on Polygonal and Polyhedral Cell Meshes". In: *ACM Transactions on Mathematical Software* 48.2 (2022), 18:1-18:23. DOI: 10.1145/3524456.

[Val19] A. Valli. "A variational interpretation of the Biot-Savart operator and the helicity of a bounded domain". In: *Journal of Mathematical Physics* 60.2 (2019), p. 021503. DOI: 10.1063/1.5024197.

[Vir+20] P. Virtanen et al. "SciPy 1.0: Fundamental algorithms for scientific computing in python". In: *Nature Methods* 17 (2020), pp. 261-272. DOI: 10.1038/s41592-019-0686-2.

[War+17] F. Warmer et al. "From W7-X to a HELIAS fusion power plant: On engineering considerations for next-step stellarator devices". In: *Fusion Engineering and Design*. Proceedings of the 29th Symposium on Fusion Technology (SOFT-29) Prague, Czech Republic, September 5-9, 2016 123 (2017), pp. 47-53. DOI: 10.1016/j.fusengdes.2017.05.034.

[Zar+01] M. C. Zarnstorff et al. "Physics of the compact advanced stellarator NCSX". In: *Plasma Physics and Controlled Fusion* 43 (12A 2001), A237-A249. DOI: 10.1088/0741-3335/43/12a/318.