Spectral approximation of the curl operator in multiply connected domains

A numerical scheme based on Nedelec finite elements has been 
recently introduced to solve the eigenvalue problem for the curl 
operator in simply connected domains. This topological assumption is not 
just a technicality, since the eigenvalue problem is ill-posed on 
multiply connected domains, in the sense that its spectrum is the whole 
complex plane. However, additional constraints can be added to the 
eigenvalue problem in order to recover a well-posed problem with a 
discrete spectrum. Vanishing circulations on each non-bounding cycle in 
the complement of the domain have been chosen as additional constraints 
in this paper. A mixed weak formulation including a Lagrange multiplier 
(that turns out to vanish) is introduced and shown to be well-posed. 
This formulation is discretized by Nedelec elements, while standard 
finite elements are used for the Lagrange multiplier. Spectral 
convergence is proved as well as a priori error estimates. It is also 
shown how to implement this finite element discretization taking care of 
these additional constraints. Finally, a numerical test to assess the 
performance of the proposed methods is reported.

A magnetic field satisfying the above equation with a constant λ is called a linear force-free field or also a free-decay field. In the theory of fusion plasma, for instance, such a field is the final state that makes the energy a minimum in order to leave the plasma in equilibrium.
From the mathematical point of view, to find a linear force-free field corresponds to solving an eigenvalue problem for the curl operator: In a bounded domain Ω, the natural boundary conditions are either u · n = 0 or u×n = 0. However, because of the Stokes theorem, the latter implies curl u·n = 0 on ∂Ω and, hence, for λ = 0, u · n = 1 λ curl u · n = 0 too. Therefore, u would vanish identically on ∂Ω and, hence, it is possible to prove that u should actually vanish on the whole domain Ω (see [17,Lemma 3]).
Thus we are in principle led to the following spectral problem for the curl operator: Find λ ∈ C and u ≡ 0 such that curl u = λu in Ω, div u = 0 in Ω, u · n = 0 on ∂Ω.
The second equation (which for λ = 0 is a consequence of the first one) rules out the trivial solutions λ = 0, u = ∇ϕ. The numerical approximation of this problem was analyzed in [14] for a simply connected domain Ω. However, when Ω is multiply connected, the set of eigenvalues of this problem is the whole complex plane (see [18,Theorem 2]).
The eigenvalue problem for the curl operator in multiply connected domains was recently analyzed in [9], where it was shown that additional constraints related to the homology of the domain have to be added for the problem to have a discrete spectrum. We consider in this paper one of the choices proposed in that reference. We introduce a variational mixed formulation of the resulting problem and a finite element discretization. We prove that the numerical approximation provides an optimal-order spectral approximation. We also discuss how to implement this numerical method and report some results for a numerical test which allows us to assess its performance.
We fix a unit normal n j on each Σ j and denote its two faces by Σ + j and Σ − j , with n j being the 'outer' normal to ∂Ω 0 along Σ + j . For any ψ ∈ H 1 (Ω 0 ), we denote by [ [ψ] ] Σj := ψ| Σ − j − ψ| Σ + j the jump of ψ through Σ j along n j . We denote by γ j the curves ∂Σ j and by t j the corresponding tangent unit vector oriented counterclockwise with respect to Σ + j . The set {γ j } J j=1 is a family of independent non-bounding cycles in the complement Ω of Ω (i.e., the union of the cycles of any non-empty subfamily cannot be the boundary of a surface contained in Ω ). A similar family of independent nonbounding cycles γ j J j=1 can be given for Ω, each cycle γ j being the boundary of a cutting surface Σ j of Ω . We denote by t j a corresponding tangent unit vector (see Figure 1). In order to obtain a well posed eigenvalue problem for the curl operator on a multiple connected domain, one additional constraint per cutting surface must be added. Whenever one can associate one cutting surface Σ j of Ω with one cutting surface Σ j of Ω (as in Figure 1), there are two alternatives for this additional constraint (see [9]): either γj u · t j = 0 or γ j u · t j = 0. We focus on the first one which, according to the Stokes Theorem, can be equivalently written as follows: Since for an eigenfunction of the curl operator with eigenvalue λ = 0, u = 1 λ curl u, we also have Σj u · n j = 0, 1 ≤ j ≤ J, which lead us to the following eigenvalue problem, whose analysis and numerical approximation is our goal: Let us remark that the last equation above makes sense since, as will be shown below (cf. (2)), the first three equations imply that u · n j ∈ L 2 (Σ j ).
3. Function spaces. In this section, we introduce some function spaces appropriate for setting and analyzing a convenient variational formulation of Problem 1. First, we recall the definitions of some classical spaces that will be used in the sequel: The spaces H(div, Ω) and H(curl, Ω) are respectively endowed with the norms defined by Here and thereafter, the symbol ⊥ is used to denote L 2 (Ω) 3 -orthogonality. We will also use the fractional Sobolev spaces H s (Ω) (0 < s < 1) endowed with the norms · s,Ω , which are well known to satisfy Let us remark that, according to [1,Proposition 3.7], there exists s > 1 2 such that H(curl, Ω) the inclusion being continuous.
Let Ω 0 := Ω\ J j=1 Σ j be the cut domain, with Σ j the cutting surfaces and n j the corresponding unit normal vectors as defined above. Let ·, · Σj denote the duality SPECTRAL APPROXIMATION OF THE CURL 239 pairing between H 1/2 (Σ j ) and H 1/2 (Σ j ). The following Green's identity has been proved in [1, Lemma 3.10].
Next, let us consider the space of the so-called harmonic Neumann fields: This is a finite-dimensional space, its dimension being equal to the number of cutting surfaces, as shown in the following lemma whose proof is essentially contained in [7, Lemma 1.3].
We can use the space of the harmonic Neumann fields to write a convenient direct decomposition of H(curl 0 , Ω).  Next step is to define another three function spaces which will play a central role in the forthcoming analysis. Let and V := X ∩ Z, endowed with · curl,Ω .
The following result follows immediately from the Helmholtz decomposition (1), the definition (5) of X and Lemma 3.4.
Moreover, we have the following characterization of the functions in X .
Proof. According to the definition (5), v ∈ X if and only if v ∈ H 0 (div 0 , Ω) (i.e., v satisfies the first two equations of the lemma) and Ω v · w = 0 for all w ∈ K T . Now, from Lemma 3.3, the latter is true if and only if Ω v · ∇φ j = 0, 1 ≤ j ≤ J, which in turn, by virtue of Lemma 3.1, is equivalent to Thus, we conclude the claimed equivalence.
Characterizations of the functions in Z and V follow immediately from this lemma.
In the following section we will introduce a variational formulation of Problem 1, which will be used for the theoretical analysis of this problem as well as for its finite element discretization. In this formulation, the eigenfunctions will be sought in the space Z. In what follows, we will establish several additional properties of this space that will be used in the sequel.
Proof. Let v ∈ D(Ω) 3 . According to Corollary 3.7, we only have to prove that curl v · n j , 1 Σj = 0, 1 ≤ j ≤ J, and this is a consequence of the Stokes theorem and the fact that, since v vanishes on Γ ⊃ γ j := ∂Σ j , The following commuting property will be the basis for the spectral characterization of the problem. It has been proved in [18,Theorem 1] in a more general context (see also [11,Prop. 2.3]). For the sake of completeness, we include an elementary proof.
Proof. The proof is based in a classical property, that in our case reads as follows: We have to prove that L also vanishes in Z. With this end, notice D( Therefore, l = − curll, so that by virtue of Lemma 3.11 it would be enough to show thatl ∈ Z to conclude the proof. To prove this, we will check that v =l satisfies the two properties from Corollary 3.7. First, since ∇ C ∞ (Ω) ⊂ C ∞ (Ω) 3 ∩ Z, we obtain from equation (8) that Then curll = −l ∈ H 0 (div 0 , Ω), and the first property from Corollary 3.7 is checked. For the second one, for each cutting surface Σ j , let U be an open connected set of R 3 containingΣ j , not intersectingΣ k , k = j, and such that U ∩ (Ω \ Σ j ) has two connected components, Ω Σ + j and Ω Σ − j , one at each side of Σ j as shown in Figure 2. Let χ ∈ C ∞ (Ω \ Σ j ) be any smooth function satisfying χ ≡ 1 in Ω Σ + j and χ ≡ 0 Then, because of Lemma 3.1 and equation Consequently, curll · n j , 1 Σj = 0, so that both properties of Corollary 3.7 are checked and, thus,l ∈ Z. As stated above, this allows us to conclude the proof.
To end this section, we will establish another density result that will be also used in the sequel. Let Φ be the subspace of smooth functions from the space Θ defined in (3), which are C ∞ up to its boundary. In general, these functions do not belong to C ∞ (Ω) because they may have (constant) jumps on each cutting surface.
Proof. We apply again the classical property used in the previous lemma. Let L be a bounded linear functional in Θ/C that vanishes on Φ/C. To conclude the density claimed in the lemma, it is enough to show that L vanishes on the whole Θ/C (cf. [8, (I.2.14)]).
Since Θ/C is a Hilbert space, there exists l ∈ Θ/C such that Hence, because of the Helmholtz decomposition (1), we have that ∇l ∈ H 0 (div 0 , Ω). Moreover, for each Σ j , let χ be as defined in the proof of the previous lemma and note that, by construction, χ ∈ Φ. Then, Lemma 3.1 yields Therefore, by virtue of Lemma 3.6, ∇l ∈ X and hence is L 2 (Ω) 3 -orthogonal to H(curl 0 , Ω) = ∇Θ (cf. Lemmas 3.5 and 3.2). Thus, for all ψ ∈ Θ/C and we conclude the proof. 4. Mixed variational formulation. The next step is to introduce a variational formulation of Problem 1. With this aim, thanks to Lemmas 3.6 and 3.5, we will impose the constraints in this eigenvalue problem by means of a Lagrange multiplier χ ∈ H(curl 0 , Ω). Moreover, the eigenfunction u will be sought in the space Z. This will ensure that the sesquilinear form on the right-hand side will be Hermitian, which will allow us to prove that the corresponding solution operator is self-adjoint. Thus, we are led to the following variational formulation.
The Babuška-Brezzi conditions for this mixed problem are easy to check. To do this, first note that thanks to Lemma 3.5 and (7), we have that v ∈ Z : Ω v ·η = 0 ∀η ∈ H(curl 0 , Ω) = Z ∩ X = V. • (inf-sup condition) there exists β > 0 such that Proof. The ellipticity in the kernel follows from the equivalence between · curl,Ω and curl · 0,Ω in the space V ⊂ X ∩ H(curl, Ω), which in turn follows from [1, Corollary 3.16] and Corollary 3.8. The inf-sup condition follows by taking v = η ∈ H(curl 0 , Ω) ⊂ Z.
Thanks to this lemma, we are in a position to define the solution operator corresponding to Problem 2: with w ∈ Z such that there exists χ ∈ H(curl 0 , Ω) satisfying By virtue of Lemma 4.2, this is a well posed problem (see, for instance, [8, Corollary I.4.1]). Moreover, clearly, T u = µu, with u ≡ 0 and µ = 0 if and only if (λ, u, 0) is a solution of Problem 2 with λ = 1/µ = 0. Therefore, our next step will be to obtain a spectral characterization of the operator T . With this end, first we establish the following additional regularity result. Consequently, T : Z → Z is compact.
Proof. Let f ∈ Z and w := T f ∈ Z, too. The same arguments used in the proof of Lemma 4.1 apply to the problem defining T and allow us to show that w ∈ X , χ ≡ 0 and curl(curl w − f ) = 0 in Ω. Then, on one side, w ∈ Z ∩ X = V (cf. (7)), so that T (Z) ⊂ V.
On the other hand, w ∈ Z ∩ X ⊂ H(curl, Ω) ∩ H 0 (div 0 , Ω) and curl w ∈ H(curl, Ω) ∩ H 0 (div 0 , Ω), too. Then, according to (2), there exists s > 1 2 such that H(curl, Ω) ∩ H 0 (div 0 , Ω) is continuously imbedded in H s (Ω) 3 . Hence, w ∈ H s (curl, Ω) and the estimate of the lemma holds true. Finally, we end the lemma from the compactness of the inclusion H s (curl, Ω) ∩ Z → Z, which in turn is a consequence of the fact that the inclusion H s (Ω) → L 2 (Ω) is compact. Proof. Let f , g ∈ Z, w := T f and v := T g. As shown in the proof of the previous lemma, w, v ∈ Z ∩ X , the corresponding Lagrange multipliers vanish and both, (curl w − f ) and (curl v − g), belong to H(curl 0 , Ω). Consequently, using that H(curl 0 , Ω) ⊥ X (cf. Lemma 3.5) and Lemma 3.11, we have that On the other hand, from the first equation of the problem defining T and Lemma 3.11 again, we have that which allow us to conclude that T : Z → Z is self-adjoint. • µ 0 = 0 is an infinite-multiplicity eigenvalue and its associated eigenspace is H(curl 0 , Ω); • {µ n } n∈N is a sequence of finite-multiplicity eigenvalues (repeated accordingly to their respective multiplicities) which converges to 0 and there exists a Hilbertian basis of associated eigenfunctions {u n } n∈N of V (i.e., such that T u n = µ n u n , n ∈ N).
Proof. From the definition of the operator T , its kernel is given by the last two equalities because of Lemmas 3.11 and 3.10, respectively. Then, since T (Z) ⊂ V (cf. Lemma 4.3) and Z = V ⊕ H(curl 0 , Ω) (cf. Lemma 3.9), the theorem follows from the two previous lemmas and the classical theory for selfadjoint compact operators (see, for instance, [4, Section 6.4]).
Taking into account the relation between the solutions of Problem 2 and the spectrum of T , we also have the following characterization of the solutions of Problem 2 and, because of Lemma 4.1, of Problem 1 too. Corollary 4.6. Problem 1 has a countable number of solutions (λ n , u n ), n ∈ N, λ n → ∞ and {u n } n∈N is a Hilbertian basis of V.
In the following section we will introduce a finite element discretization of Problem 2. In order to prove the convergence of the proposed numerical scheme we will use the following operator: where (w, ξ) ∈ Z × Θ/C is the solution of the following problem: By virtue of Lemma 3.2, H(curl 0 , Ω) = ∇Θ. Therefore, this problem is equivalent to the one used to define the operator T with χ = ∇ξ. Hence, G(f , g) = (T f , 0) and we have the following result.

5.
Finite element spectral approximation. In this section, we introduce a Galerkin approximation of Problem 2 and prove convergence and error estimates for the approximate eigenvalues and eigenfunctions.
With this end, we assume that Ω is a polyhedron and choose the cutting surfaces Σ j , 1 ≤ j ≤ J, also polyhedral. Let {T h } h>0 be a regular family of tetrahedral partitions ofΩ compatible with the cutting surfaces in the sense that each Σ j is a union of faces of tetrahedra T ∈ T h . Therefore, each T h can also be seen as a mesh of the cut domain Ω 0 . The mesh parameter h denotes the maximum diameter of all the tetrahedra of the mesh T h . From now on we will denote by C a generic constant, not necessarily the same at each occurrence but always independent of the mesh parameter h.
We use classical curl-conforming Nédélec finite elements: where, for k ∈ N, P k is the space of polynomials of degree not greater than k and P k the subspace of homogeneous polynomials of degree k. Let which by virtue of Corollary 3.7 can be written as follows: We denote by I N h the so-called Nédélec interpolation operator. We refer to [13,Section 5.5] for its precise definition and the properties of this interpolant that we will use in the sequel. This interpolant is well defined for functions in H s (curl, Ω) provided s > 1 2 , so that I N h : H s (curl, Ω) → N h is a bounded linear operator. Moreover, as it is shown in the following lemma, the Nédélec interpolant of functions from Z remains in this space.
Proof. Let I R h be the divergence-conforming Raviart-Thomas interpolation operator (see, for instance, [13,Section 5.4] for its definition and properties). This interpolant is well defined for functions in H s (Ω) 3 with s > 1 2 ([13, Lemma 5.15]). Let v ∈ H s (curl, Ω) ∩ Z with s > 1 2 . According to Corollary 3.7, curl v · n = 0 on Γ and Σj curl v · n = 0, 1 ≤ j ≤ J (note that the integral makes sense because curl v ∈ H s (Ω) 3 with s > 1 2 ). Therefore, curl I N h v · n = I R h curl v · n = 0 on Γ, where the first equality follows from [13,Lemma 5.40] and the second one from the fact that the Raviart-Thomas interpolant preserves vanishing normal components on the faces of the tetrahedra of the mesh (which in turn follows from the definition of this interpolant). Analogously, it also preserves the integral of these normal components, so that we have Thus, To discretize the Lagrange multiplier χ ∈ H(curl 0 , Ω), we use Lemma 3.2 to write χ = ∇ϕ with ϕ ∈ Θ, and approximate this space by The following discrete version of Lemma 3.2 has been proved in [3, Lemma 5.5].

Lemma 5.2. There holds
Now we are in a position to introduce a finite element discretization of Problem 2.
First of all note that, as in the continuous problem, for any solution (λ h , u h , ϕ h ) of Problem 3, ϕ h vanishes. In fact, since ∇Θ h ⊂ Z h (cf. Lemma 5.2), this follows by taking above v h = ∇ϕ h .
Taking arbitrary bases of the finite dimensional spaces Z h and Θ h /C, Problem 3 can be equivalently written as a degenerate generalized matrix eigenvalue problem in which both matrices are Hermitian but none is positive definite. However, we will show in the next section that it is also equivalent to another generalized matrix eigenvalue problem of smaller dimension which will be proved to be well posed.
In order to prove that the eigenvalues and eigenfunctions of Problem 2 are well approximated by those of Problem 3, we resort to the classical theory for mixed eigenvalue problems of the so-called type (Q1) reported in [12,Section 3]. With this aim, we have to check the following properties, which correspond to assumptions (3.12)-(3.16) from this reference.
• Babuška-Brezzi conditions for the continuous problem; they have been already established in Lemma 4.2. • Babuška-Brezzi conditions for the discrete problem, namely: -(ellipticity in the discrete kernel) there exists α * > 0 such that it has been proved in [1,Proposition 4.6]. -(discrete inf-sup condition) there exists β * > 0 such that it follows by taking v h = ∇ψ h , which according to Lemma 5.2 lies in Z h . • Density of the finite element spaces: for all (v, ψ) ∈ Z × Θ, it follows from the densities of C ∞ (Ω) 3 ∩ Z in Z (cf. Lemma 3.12) and Φ/C in Θ/C (cf. Lemma 3.13), the fact that for a smooth v ∈ Z its Nédélec interpolant I N h lies in Z h (cf. Lemma 5.1) and standard approximation properties of the Nédélec and the Lagrange interpolants (cf. [13,Theorem 5 An additional hypothesis assumed in the spectral approximation theory from [12,Section 3] is the compactness of the solution operator G defined in (9), which in our case has been already established in Lemma 4.7. Moreover, this theory also involves a formal adjoint operator G * , which in the present case is defined for (f , g) ∈ Z × Θ/C by G * (f , g) := (w * , ξ * ), with (w * , ξ * ) ∈ Z × Θ/C being the solution of the adjoint problem: However, according to Lemma 3.11, for f , v ∈ Z, Ω v · curlf = Ω curl v ·f and, then, in this case, the formal adjoint G * coincides with G. Therefore, all the hypotheses needed to apply [12, Theorem 3.1] hold true, which allows us to establish the main result of this paper. Let E h × {0} be the direct sum of the eigenspaces corresponding to λ is the distance from the continuous eigenspace E to the discrete space Z h and is the so-called gap between the continuous and the discrete eigenspaces.
To end this section we establish an appropriate estimate for the term γ h . With this end, recall that k ≥ 1 is the degree of the Nédélec finite elements and let s > 1 2 be a Sobolev exponent such that the inclusion (2) holds and is continuous. Proof. Let v ∈ E be such that v curl,Ω = 1 and µ := 1 λ , so that T v = µv. Then, from Lemma 4.3 it follows that v ∈ H s (curl, Ω) and Therefore, using again the standard error estimate for the Nédélec interpolant I N h v (cf. [13,Theorem 5.41(1)]) and taking into account that according to Lemma 5.1 The result follows from the definition of γ h and these two inequalities.
As a consequence of the two previous theorems we conclude that the eigenvalues and eigenfunctions of Problem 3 converge with optimal order to those of Problem 2.
6. Implementation issues. In this section we will briefly describe how to impose in N h the constraints defining Z h : curl v h · n = 0 on Γ and With this end, the results from [11,Section 5] and [14,Section 5] have been extended to multiply connected domains in [10, Section 5] yielding the following result.
where Γ := Γ \ J j=1 ∂Σ j and ∇ Γ denotes the surface gradient on Γ. The next step is to introduce a convenient basis of the space Θ h . Consider the standard finite element discretization of H 1 (Ω): On the other hand, for j = 1, . . . , J, let ϕ j ∈ Θ h be such that, for all nodes P , It is easy to check that The final step is to introduce a basis of Z h . With this aim, for simplicity, we assume that the boundary Γ is connected (otherwise, the same procedure should be repeated for each of its connected components).
Without loss of generality we order the nodal basis functions {ϕ j } L j=1 of L h so that the first K of them correspond to nodal values on Γ. Then, it is easy to check that Since surface gradients are determined up to an additive constant, we choose one vertex on Γ and drop out the basis function corresponding to this vertex (e.g., vertex number K).
Finally, let {Φ m } M m=1 be the canonical basis of N h ; without loss of generality, we assume that {Φ m } M m=M +1 are those related to the faces or edges on Γ. Then we have the following result whose proof can be found in [10, Proposition 5.1].
Let us remark that the matrices of the algebraic eigenvalue problem corresponding to the discrete mixed formulation can be easily obtained by static condensation from the matrices of the classical Nédélec and P k -continuous elements (see [10, Section 5.1] for more details). The resulting algebraic generalized eigenvalue problem has the form where the entries of u and ϕ are the components in the above given bases of u h and ϕ h , respectively. Both matrices above are Hermitian (actually, real symmetric), but none is positive definite. Since ϕ = 0, the above problem is equivalent to which, in turn, is equivalent to with a real symmetric and positive definite left-hand side matrix. This allows us to conclude that Problem 3 is well posed and that it has dim(Z h ) − dim(ker(C)) non-zero eigenvalues (repeated accordingly to their respective multiplicities).
7. Numerical results. We have implemented the method analyzed above for the lowest-possible order (k = 1) in a matlab code. We have applied it to compute the smallest positive eigenvalues in a toroidal domain as that shown in Figure 3, with r 1 = 1 and r 2 = 0.5, for which no analytical solution is available. We have used meshes T h with different levels of refinement; we identify each mesh by the corresponding number N h of tetrahedra. For each computed eigenvalue we have estimated the order of convergence and a more accurate value by means of a least-squares fitting of the model λ h,k ≈ λ ex + Ch t with h = N −1/3 h . Table 1 shows the obtained results for the four smallest eigenvalues. Note that the first two converge to a double-multiplicity eigenvalue of the continuous problem and the other two to another eigenvalue of multiplicity two.
A quadratic order of convergence can be seen from Table 1, which agrees with the theoretical results for the finite elements used (k = 1) and a smooth domain (s ≥ 1). Figure 4 shows a log-log plot of the computed errors for the eigenvalue λ 1 versus the number of tetrahedra. Once more, the quadratic order of convergence can be appreciated.