Div-Curl Problems and $\mathbf{H}^1$-regular Stream Functions in 3D Lipschitz Domains

We consider the problem of recovering the divergence-free velocity field ${\mathbf U}\in\mathbf{L}^2(\Omega)$ of a given vorticity ${\mathbf F}=\mathrm{curl}\,{\mathbf U}$ on a bounded Lipschitz domain $\Omega\subset\mathbb{R}^3$. To that end, we solve the"div-curl problem"for a given ${\mathbf F}\in{\mathbf H}^{-1}(\Omega)$. The solution is expressed in terms of a vector potential (or stream function) ${\mathbf A}\in{\mathbf H}^1(\Omega)$ such that ${\mathbf U}=\mathrm{curl}\,{\mathbf A}$. After discussing existence and uniqueness of solutions and associated vector potentials, we propose a well-posed construction for the stream function. A numerical method based on this construction is presented, and experiments confirm that the resulting approximations display higher regularity than those of another common approach.


INTRODUCTION
Let Ω ⊂ ℝ 3 be a bounded Lipschitz domain. Given a vorticity field ( ) ∈ ℝ 3 defined over Ω, we are interested in solving the problem of velocity recovery: This problem naturally arises in fluid mechanics when studying the vorticity formulation of the incompressible Navier-Stokes equations. Vortex methods, for example, are based on the vorticity formulation and require a solution of problem (1) at every time-step. 1 While our motivation lies in fluid dynamics, this "div-curl problem" also is interesting in its own right.
On the whole space ℝ 3 , this problem is a classical matter. Whenever is smooth and compactly supported, the unique solution of problem (1) that decays to zero at infinity is given by the Biot-Savart law. 2, Proposition 2. 16 However, the case where Ω is a bounded domain is significantly more challenging.
In numerical simulations of the incompressible Navier-Stokes equations, it is common to fulfil the constraint div = 0 only approximately, but it has recently been demonstrated that such a violation can cause significant instabilities. The importance for numerical methods to fulfil this constraint exactly was stressed by John et al. 3 One way of achieving this requirement is the introduction of a stream function, or vector potential: instead of solving problem (1) directly, one seeks an approximation ℎ of an auxiliary vector-field such that = curl . Because of the vector calculus identity div • curl ≡ 0, the velocity field ℎ = curl ℎ is always exactly divergence free. In particle methods, the particle positions ∈ Ω, = 1, … , , are updated by solving d d ( ) = , ( ) , = 1, … , using a time-stepping scheme. It makes sense to use a volume-preserving scheme for this problem. However, most of these schemes require a stream function and not the velocity as input, 4, Chapter VI. 9 and thus arises the desire to have a stream function of maximum regularity at hand. In this work we describe how to compute stream functions that are at least 1 (Ω)-regular-even on non-smooth domains-thereby improving on the regularity of approximations currently available in related literature. arXiv:2005.11764v2 [math.AP] 3 Feb 2021

Uniqueness. (Theorems 2 and 4)
If Ω is "handle-free", the solution ∈ 2 (Ω) of problem (1) can be made unique by prescribing its normal trace ⋅ ∈ − 1 2 (Γ). Moreover, if the prescribed boundary data fulfils ∫ Γ ⋅ d = 0 on each connected component Γ ⊂ Γ of the boundary, there exist conditions that uniquely determine a stream function ∈ 1 (Ω) such that = curl . 4. Construction of Solutions. (Section 5) The main novelty of this work lies in the explicit construction of 1 (Ω)-regular stream functions directly from the vorticity . Given a vorticity ∈ −1 (Ω) fulfilling the conditions of item 1 and boundary data ⋅ ∈ − 1 2 (Γ) fulfilling the conditions of item 2, this construction will yield a stream function ∈ 1 (Ω) such that = curl solves problem (1). If the domain is handle-free, the obtained solution will be the uniquely defined stream function ∈ 1 (Ω) from item 3. Numerical methods will be described in Section 7.
The results concerning existence and uniqueness of velocity fields follow from classical functional analytic arguments and are well-known, but covered for completeness. Existence of stream functions ∈ 1 (Ω) is due to Girault and Raviart 5,Theorem 3.4 . Their work, however, left unclear how to compute such a potential.

Problematic Approaches
A naive approach to the div-curl problem (1) relies on the observation that Based on this vector-calculus identity, it is tempting to solve three decoupled scalar Poisson problems -Δ = (curl ) , = 1, 2, 3, for the components of , say by prescribing the value of each one on the boundary. However, this approach is problematic: it is our aim to integrate , but instead this strategy asks that we differentiate first. Therefore, it needlessly requires to impose more regularity on the right-hand side. Moreover, there is no guarantee that its solution is divergence-free. Finally, since the tangential components of allow us to compute (curl ) ⋅ on the boundary, the boundary data must fulfil the compatibility condition (curl ) ⋅ = ⋅ . We will later see that the solutions of problem (1) are usually not 1 (Ω)-regular, and thus the classic existence and uniqueness results in 1 (Ω) for the scalar Poisson problems -Δ = (curl ) are not applicable either. Another straightforward approach assumes that ∈ 2 (Ω). One may then extend by zero to the whole space, yielding ∈ 2 (ℝ 3 ), and apply the Biot-Savart law to this extension. The normal trace ⋅ on Γ can then be prescribed by adding a suitable "potential flow". The main caveat of this strategy is that unless ⋅ = 0 on the boundary, the zero extensioñ will not be divergence-free, and in this case the Biot-Savart law fails to yield the correct result. We will later see that this approach can in fact be fixed by introducing a suitable correction on the boundary.

Our Results in Context
Clearly, for a given velocity field , the condition = curl alone does not uniquely determine : because of the vector calculus identity curl • (-∇) ≡ , any gradient may be added to without changing its curl. It is thus natural to enforce the so-called
These systems differ in the involved spaces and boundary conditions. For the -system we would like to prescribe ⋅ on Ω, while for the -system we actually do not care which boundary conditions are prescribed, as long as they ensure that the solution is unique, and-hopefully-as regular as possible.
Many results in the literature are concerned with only one of these systems, an overview of some results is given in Table 1. For example, the famous work of Amrouche et al. 6 is concerned with the -system and ∈ 2 (Ω). They propose tangential or normal boundary conditions for and numerical methods to approximate the resulting stream functions. However, even for perfectly smooth velocity fields , the resulting potentials will usually only have Sobolev regularity 1 2 (Ω), unless the domain is assumed to be more regular than just Lipschitz. In particular, functions from (curl; Ω) ∩ 0 (div; Ω) may develop quite strong singularities near corners of the domain, which makes it difficult to approximate them efficiently. 7, Figure 1.3 Note that these singularities can only occur in non-smooth domains: for 1,1 -domains one has T , N ∈ 1 (Ω), so in this case there is no need to look for more regular potentials. Almost all of the literature concerning numerics for the -system considers either T or N . 8, Chapter 6.1 Many authors also consider more general problems involving inhomogeneous material coefficients, 9 or the -setting, 10,11 which on the other hand again often requires higher regularity assumptions on the boundary. It has long been known that stream functions of regularity 1 (Ω) do exist, but so far conditions that uniquely characterise them have not been given. We are unaware of any previous approaches that allow us to efficiently compute ∈ 1 (Ω) for the case of general Lipschitz domains. Our work aims to close this gap. It proposes natural conditions which uniquely determine a vector potential 1 ∈ 1 (Ω) without explicitly involving boundary values of 1 . We believe it is because previous approaches do prescribe boundary conditions like T ⋅ = 0 or N × = that they yield less regular stream functions.
Our algorithm utilises the Newton operator: the bulk of the work lies in the explicit computation of a volume integral. Two companion scalar elliptic equations must also be solved on the boundary of the domain, but these are easily tackled using standard boundary element methods. While this construction solves both the -and -systems simultaneously, we first discuss existence and uniqueness of solutions for each of them individually.

Trace Spaces
We refer to McLean 14, Chapter 3 , Sauter and Schwab 15, Chapter 2 and Buffa et al. 16 for theory concerning extension of the traces to continuous and surjective mappings having continuous right-inverses (so-called lifting maps). The fact that the operators and are surjective is one of the main results of Buffa et al. 16 It allows for the derivation of Hodge decompositions, which we will also use later on. The surface differential operators: used in those definitions are defined by duality for all ∈ Here -∇ Γ and Γ are suitable extensions of the ordinary trace gradient and trace curl to operators with mappings -∇ Γ ∶ 1 2 (Γ) → T (curl Γ ; Γ). Finally, we will also make use of the Laplace-Beltrami operator -Δ Γ ∶= div Γ • (-∇ Γ ) ≡ curl Γ • Γ . This operator is known to be coercive on 1 (Γ)∕ℝ, the space of 1 (Γ)-regular traces with zero average on each Γ , = 0, … , 2 .

Geometry of the Domain
Throughout this article, we suppose that the Lipschitz domain Ω ⊂ ℝ 3 of interest is bounded and connected, and we will sometimes refer to the Betti numbers 1 and 2 . These numbers are related to the topological properties of the domain, see Figure 1. 1 is the genus of the domain, in other words the number of "handles". 2 is the number of 'holes' Θ , = 1, … , 2 in the domain. We define the exterior domain Θ 0 via: The domain's boundary thus always has 2 + 1 connected components Γ ∶= Θ , = 0, … , 2 . A domain with 1 = 0 is called handle-free, 1 hole-free if 2 = 0, and we say that the topology of Ω is trivial or simple if 1 = 2 = 0. We refer to Arnold and al. for more details. 17,Section 2 The geometric interpretation of the Betti numbers is best illustrated through an example. In the domain depicted in Figure 1,

FIGURE 1
A ring-shaped domain with three cubical holes is an example of a domain having non-trivial topology. 6, Section 3 There is one "handle" through it: 1 = 1. The red line is a representative of the equivalence class of non-bounding cycles. The three cubical inclusions ("holes") are not part of the domain: 2 = 3. The boundary Γ has four connected components: Γ 1 , Γ 2 , Γ 3 , and the domain's exterior boundary Γ 0 .
On the one hand, the value of the second Betti number 2 is relevant to questions regarding the existence results stated in item 1 and item 2 of Section 1. These existence theorems will make use of arbitrary but fixed functions ∈ ∞ 0 (ℝ 3 ), = 0, … , 2 , that act as indicators for the the connected components of the boundary: On the other hand, the value of 1 is crucial to the uniqueness results of item 3. For simplicity, we will restrict our attention to handle-free domains ( 1 = 0), that is domains for which every loop inside Ω is the boundary of a surface within Ω. The domain in Figure 1 is not handle-free, as the red loop is a representative of the equivalence class of non-bounding cycles. In that example, 1 = 1. Nevertheless, we will make some remarks on what changes in the following results need to be anticipated in order to recover uniqueness of solutions when 1 > 0.

Laplace and Newton Operator
The scalar Laplacian -Δ ∶= div • (-∇) is central to potential theory, and we assume that the reader is well aware of the classical existence and uniqueness results for boundary value problems related to this operator. Its vector-valued analogue is defined component-wise: This operator is also known as Hodge-Laplacian and can equivalently be written as: Let us denote by ( ) ∶= (4 | |) −1 the fundamental solution of the Laplacian -Δ. The Newton operator is defined on the space of compactly supported distributions  ′ (ℝ 3 ) via convolution as  ∶  ′ (ℝ 3 ) →  ′ (ℝ 3 ),  → ⋆ . In other words: . For a given ∈  ′ (ℝ 3 ), the distribution  is called the Newton potential of . This operator is an inverse for the Laplacian: 18 Moreover, because it is an operator of convolutional type, it commutes with differentiation. 18, Equation (4.2.5) Application of this operator always increases the Sobolev regularity of a distribution ∈ (ℝ 3 ) ∩  ′ (ℝ 3 ) by two, i. e.,  has the following mapping property and is continuous: 15 The Newton potential  of a compactly supported distribution ∈  ′ (ℝ 3 ) is regular outside supp and decays to zero at infinity: 19 Even more importantly, the following result characterises the Newton potential. 19, Chapter II, §3.1, Proposition 3 Lemma 1. Let ∈  ′ (ℝ 3 ) and ∈  ′ (ℝ 3 ). Then is the Newton potential of , that is =  , if and only if: This characterisation allows for the derivation of representation formulae for solutions of the Laplace equation on bounded domains, leading to boundary integral equations. This will appear at the end of this section. All of these results analogously hold for the vector-valued Newton operator , which is defined component-wise and distinguished by a bold-face font.

Decompositions of Helmholtz-Hodge Type
Decomposition theorems will play a central role in our analysis, so we collect the most important results here.
Note that in this decomposition neither curl nor -∇ are necessarily compactly supported. This makes the following result useful, for which we refer to the works of Bramble and Pasciak 12, Proposition 3.2 , Pasciak and Zhao 20, Lemma 2.2 , as well as Hiptmair and Pechstein 21, Remark 3 .

Trace Jumps and a Representation Formula
The trace operators introduced in the previous sections were all defined with respect to the domain Ω. One can instead also consider the corresponding traces with respect to the complementary domain Ω ∶= ℝ 3 ⧵ Ω. These one-sided traces exist whenever the restriction of a vector field ∈ 2 loc (ℝ 3 ) to the domains Ω and Ω is sufficiently smooth. If is smooth across Γ, then the one-sided traces coincide. Otherwise the difference of these traces is denoted by the jump operator ⋅ . For example, for one writes: The importance of the jump operator lies in the following representation formula.
Proof. The fact that = (-) directly follows from Lemma 1. Because div = 0 globally, we have -= curl(curl ) and curl(curl )| ℝ 3 ⧵Γ = . We thus obtain for all ∈ (ℝ 3 ): where we used that curl ∈ 2 loc (ℝ 3 ) since ∈ 1 loc (ℝ 3 ). Now, by definition of the rotated tangential trace: Applying the same methodology to the integral over ℝ 3 ⧵ Ω and using the definition of yields the desired result, because the fact that is smooth across Γ guarantees that = .

VELOCITY FIELDS
In this section we prove the existence of velocity fields solving (1) as claimed in Item 1. The abstract integrability condition is reformulated in Lemma 8. The uniqueness result for velocity fields presented in Item 3 is also covered.

Existence of Velocity Fields
has a solution ∈ 2 (Ω) if and only if lies in the space −1 (Ω) and fulfils the following integrability condition: Remark 1. A more general result was proven by Bramble and Pasciak, 12, Theorem 4.1 using entirely different techniques. We follow another route that sheds new light on the classical conditions imposed on and that helps in the construction of vector potentials later.
has closed range. 7, Box 3.1 The curl operator is symmetric and the dual of the mapping (27) is the operator curl | | | 2 (Ω) given in (24). Hence, Banach's closed range theorem yields That is, Evidently, problem (26)  . Thus, together with Lemma 6, the claim follows.
Proof of Theorem 1. Lemma 7 guarantees the existence of a ∈ 2 (Ω) such that curl = . This function does not necessarily fulfil div = 0. But in this case we let ∈ 1 0 (Ω) denote the unique solution to the Poisson problem and note that ∶= + ∇ solves the div-curl system (21).
The integrability condition (22) is most natural for the chosen method of proof. However, it is hard to verify in practice. For this reason, it is worthwhile considering equivalent alternative conditions. Lemma 8. Suppose that Ω ⊂ ℝ 3 is a bounded Lipschitz domain and let ∈ −1 (Ω). Together, the following conditions are equivalent to the integrability condition (22): If in particular ∈ 2 (Ω) and div = 0, condition (31b) is equivalent to: Remark 3. Notice that definition (8) guarantees that -∇ | | |Ω ∈ (Ω).
Remark 5. The normal trace ⋅ is well-defined because the solution of the div-curl system satisfies div = 0.
Proof. Let us first remark that the condition ⟨ , 1⟩ Γ = 0 is necessary. To see this, note that because div = 0, any solution ∈ 2 (Ω) of the div-curl system (21) must fulfil: Now let ∈ 2 (Ω) denote any solution of the div-curl system, whose existence is guaranteed by Theorem 1. Let ∈ 1 (Ω)∕ℝ be the unique solution of the Neumann problem: Then the function ∶= − ∇ fulfils the conditions of the theorem. To see that it is unique, let 1 , 2 ∈ 2 (Ω) denote two solutions of the div-curl system (21) that fulfil 1 ⋅ = 2 ⋅ = on the boundary Γ. Then their difference ∶= 1 − 2 solves In other words, is a so-called Neumann harmonic field. These functions form a space of dimension 1 . 6, Proposition 3.14 7, Section 4.3 By hypothesis 1 = 0, and thus = .
The above proof hints at what needs to be done in order to recover uniqueness in the case 1 ≠ 0. One needs to prescribe 1 functionals that determine the Neumann harmonic components of . A construction of these fields and corresponding functionals can be found in the work of Amrouche et al. 6

STREAM FUNCTIONS
We prove the existence result for stream functions of Item 2 in this section. The related uniqueness statement of Item 3 is proven in Theorem 4.

Existence of Stream Functions
The following theorem is a variant of a result by Girault and Raviart. 5, Theorem 3. 4 We give a different proof, which uses the Newton operator instead of Fourier transforms.

Theorem 3.
Let Ω ⊂ ℝ 3 denote a bounded Lipschitz domain. Then ∈ 2 (Ω) satisfies if and only if there exists a vector-field ∈ 1 loc (ℝ 3 ) such that Proof. "⇐" Note that the conditions (41a) and (41b) are exactly the integrability conditions (31a) and (32). Because of Lemma 7 and Lemma 8, these conditions are necessary to ensure the existence of a vector-field ∈ 2 (Ω) such that = curl in Ω. "⇒" In order to show sufficiency, the idea is to extend to ℝ 3 by "potential flows" matching ⋅ on Γ, then use the Newton operator.
We want to exploit the following scalar functions. For = 0, we let 0 ∈ 1 loc (Θ 0 ) denote the solution of the problem: whereas for = 1, … , 2 we define ∈ 1 (Ω )∕ℝ as the solution of: Because of condition (41b), it is well-known that these problems are well-posed.

Uniqueness of Stream Functions
For a given velocity field , a stream function can be constructed as in the proof of Theorem 3. In the following section, however, we construct a stream function directly from , so there the extended velocity field is a-priori unknown. The following result allows us to establish that the stream functions from Theorem 3 and Section 5 coincide if 1 = 0. This in turn will also allow us to establish higher regularity in Theorem 6. Note that one obtains uniqueness of , while neither explicitly referring to the extension of outside Ω, nor to the boundary values of on Γ. The proof also gives motivation for the constructions presented later.
Let us now make some remarks on the case 1 ≠ 0. We define Ω ∶= ℝ 3 ⧵ Ω as the complementary domain of Ω, and ∶= curl | Ω and ∶= curl | Ω . These functions are Neumann harmonic fields: Ultimately, the idea is to rely on the fact that in handle-free domains the space of Neumann harmonic fields only contains the zero element, and thus = and = . In the case 1 ≠ 0, however, neither Ω nor Ω are handle-free, and in fact we have 1 = 1 . The spaces of Neumann harmonic fields on Ω and Ω then each have dimension 1 . Buffa has derived the analogue of Lemma 4 for the case of Lipschitz polyhedra with 1 ≠ 0. 23 Because 1 = 1 , it contains an additional term from the 2 1 -dimensional space of harmonic tangential fields. Half of these components are fixed because of the condition curl | Ω = , the other half concerns the external harmonic fields. To ensure uniqueness of , one additionally needs to prescribe the Neumann harmonic components of ∶= curl | Ω .

CONSTRUCTION OF SOLUTIONS AND WELL-POSEDNESS
In this section, we provide a construction for a stream function ∈ 1 loc (ℝ 3 ) for the general case of a given vorticity field ∈ −1 (Ω). This construction may also be considered an alternative proof of the existence results of Theorems 1 and 3. The idea is to first find a suitable extensioñ ∈ −1 (ℝ 3 ) of , that additionally satisfies diṽ ∈ −1 (ℝ 3 ). The spurious divergence of̃ can then be cancelled out using a surface functional, and the problem can be solved by applying the Newton operator. In order for this approach to work it is crucial to make use of Lemma 6 and to interpret as a member of 0 (curl; Ω) ′ . Otherwise, one will usually only obtain diṽ ∈ −2 (ℝ 3 ) and the construction will fail.
In computational practice, one will often have ∈ 2 (Ω). Under this assumption we can simplify the construction, yielding an algorithm that is more easily implementable.

(69)
Then solves the div-curl system and ⋅ = on Γ, and is a stream function for . In the case of a handle-free domain, Theorems 2 and 4 guarantee that these functions are unique. Moreover, the solution continuously depends on and . In total we have therefore proven the following theorem.
These functions linearly and continuously depend on the data and , and we have: . (70)

Remark 6.
There are extensions of this result to domains with 1 ≠ 0, for example in Monk's book. 25,Theorem 3.47 However, this extension is lacking the statement ⋅ ∈ 2 (Γ) ⇐⇒ × ∈ 2 (Γ) which we need for our proof below. We thus refer to the original work of Costabel for 1 = 0, but one can expect that this equivalence also generalises to the case 1 ≠ 0. Under this assumption the following regularity result remains true if zero Neumann harmonic components are prescribed for = curl | ℝ 3 ⧵Ω .
This result can directly be applied to velocity fields solving the div-curl system (21). In the following, we show that it also implies higher regularity of the associated stream functions .
Then, the unique solution ∈ 2 (Ω) of the div-curl system (21) with ⋅ = on Γ fulfils ∈ 1 2 (Ω), and its uniquely determined stream function from Theorem 4 fulfils ∈ 3 2 loc (ℝ 3 ⧵ Γ). Proof. The regularity of is exactly Costabel's result Lemma 11. For the regularity of , we first consider the case = 0. Thus, let 0 denote the unique solution of the div-curl system (21) that satisfies 0 ⋅ = 0 on Γ, and let 0 ∈ 1 loc (ℝ 3 ) denote its associated stream function. An application of the representation formula for the vector Laplacian then yields: wherẽ ∈ 2 (ℝ 3 ) is 's zero extension as defined in (71). Clearly, from the mapping properties of the Newton operator, it follows that̃ ∈ 2 loc (ℝ 3 ). For the boundary term we note that from the construction of 0 in the proof of Theorem 3 it is clear that curl 0 = in ℝ 3 ⧵ Ω. This implies that and because of Lemma 11 this yields 0 ∈ 2 (Γ). The boundary term ′ 0 may thus alternatively be interpreted as a componentwise application of the scalar single layer potential operator  ′ to the components of 0 . For this operator the following mapping property is known: 15, Remark 3.1.18b and thus ′ 0 ∈ 3 2 loc (ℝ 3 ⧵ Γ). For general boundary data ⋅ = ∈ 2 (Γ), one then needs to solve the hypersingular boundary integral equation: for the unknown ∈ 1 2 (Γ)∕ℝ and set ∶= 0 + curl Γ , ∶= (̃ + ′ ). For the integral equation the following regularity result is known: 15, Theorem 3.2.3b ∈ 2 (Γ) ⇐⇒ ∈ 1 (Γ).

NUMERICAL APPROXIMATIONS
The case where ∈ (Ω) for some > 2, leads to a particularly simple discretisation using Raviart-Thomas elements. In practice this condition on is often not a real restriction. Here, we will only give a brief sketch of this scheme and its analysis. The case ∈ −1 (Ω) is technically more involved, and we restrict ourselves to giving some remarks on possible numerical realisations. At the end of this section we give a numerical example with ∈ ∞ (Ω) and ∈ ∞ (Ω). Even for such smooth data, the tangential T ∈ (curl; Ω) ∩ 0 (div; Ω) shows quite strong singularities. The newly proposed stream function ∈ 1 (Ω) is not smooth either, but displays increased regularity compared to T .

Meshes and Spaces
For our numerical approximations we will assume that the domain Ω is polyhedral, handle-free, and that a family { ℎ } ℎ>0 of shape-regular, quasi-uniform, tetrahedral meshes is available. We will write ∈  ℎ for the tetrahedra of such a mesh, the parameter ℎ > 0 refers to the average diameter of these tetrahedra. On these meshes we respectively define standard Lagrangian elements, Nédélec elements, Raviart-Thomas elements, and discontinuous elements of order ∈ ℕ: boundaries is "concentrated" in the term ‖ − ℎ ‖, and for ∈ 2 (Γ), as a consequence of the regularity result Theorem 6, we may at most expect convergence of order (ℎ 1 2 ). However, also note that (85) is an a-posteriori bound that is easily computable. The term ‖ − ℎ ‖, for example, can be made arbitrarily small by increasing the mesh resolution for ℎ . Because ℎ is the solution to a scalar equation on the boundary, adaptive refinement strategies are a comparatively cheap remedy to tackle the irregularity.

Remarks on the General Case ∈ −1 (Ω)
It is straightforward to discretise the general construction given in Section 5.1. An obvious choice would be a standard Galerkin method and finding ℎ ∈ ℎ,0 (Ω) such that ( ℎ , ℎ ) = ⟨ , ℎ ⟩ for all ℎ ∈ ℎ,0 (Ω). Ultimately, however, the rate of convergence for such a method would less depend on the regularity of itself, but more on that of its Riesz representative . If happens to be smooth, can turn out to be much less regular than the function it represents.
Note that the right side may be immediately replaced with ⟨ , ℎ ⟩, as ℎ ∈ 0 (curl; Ω) and ∈ 0 (curl; Ω) ′ by Lemma 6. This is not possible in the corresponding method for the tangential potential T , where the test functions do not have a vanishing tangential trace.
The true solution fulfils div N = 0 in Ω, it then suffices to cancel the spurious divergence of its zero extensioñ N . We have diṽ N = − ′ N ⋅ , so we may solve -Δ Γ = N ⋅ , set ∶= -∇ Γ , and obtain̂ ∶=̃ N + ′ ∈ 1 loc (ℝ 3 ). The numerical approximation N,ℎ will usually not be exactly divergence free, here the jumps of the normal traces on the internal faces will also need to be cancelled out in order to achieve 1 (Ω)-regularity. From here we may proceed analogously as before: for the correction of the normal trace of curl̂ we solve a boundary integral equation.

An Example Illustrating Higher Regularity
We consider the domain Ω ∶= (0, 1) 3 ⧵[0.1, 0.8] 3 and the smooth velocity field ∈ ∞ (Ω) associated to the vorticity ∈ ∞ (Ω) given by: The domain Ω was chosen to be asymmetric, non-smooth, non-convex, and topologically non-trivial ( 2 = 1), while at the same time being "easy" from the viewpoint of meshing. Neither for the tangential potential T , nor for the potential introduced in this work explicit expressions are known. Thus, the finite element method by Amrouche et al. has been implemented to compute T . 6, Equation (4.5) We use order = 2 for T,ℎ ∈ ℎ (Ω), such that the velocity field is recovered exactly. For the other stream function, notice that in this case the Newton potential: can be evaluated directly, without further discretisation. 28 The Laplace-Beltrami equation for ℎ and ℎ , as well as the hypersingular boundary integral equation for ℎ , are discretised as described above, with order = 2.
We remark that the method by Amrouche et al. requires to use the velocity field as input, while the approach discussed in this work only requires and the boundary data ⋅ . The numerical results along the line ( 1 , 0.5, 0.8) ⊤ are shown in Figure 2. In the interval 1 ∈ [0.1, 0.8] this line touches the internal boundary Γ 1 . One clearly sees that close to the corners at 1 = 0.1 and 1 = 0.8 neither solutions are smooth. But while the tangential potential T develops very steep gradients and jumps near 1 = 0.1, the new potential only exhibits a small kink, which suggests more regularity.

CONCLUSIONS AND OUTLOOK
In this work, we have established precise conditions under which a divergence-free velocity field ∈ 2 (Ω) can be recovered from its given curl and boundary data ⋅ . Additionally, minor complementary assumptions on the boundary data guarantees that this velocity field can be represented in terms of a stream function ∈ 1 (Ω), which can be explicitly constructed. This stream function is more regular than the tangential vector potential suggested by Amrouche et al. 6 The regularity result of Theorem 6 is sharp in several ways. Let us for example consider the case of a handle-free domain Ω with 1 = 0 and suppose that ≡ . It is then a classical result that the velocity field can be written in terms of the gradient of a scalar potential: = -∇ , where -Δ = 0 in Ω. Even if the given boundary data ⋅ is smooth, it is known from the regularity theory for the scalar Laplace equation that, on general Lipschitz domains, the highest regularity one can expect is ∈ 3 2 (Ω), which therefore only leads to ∈ 1 2 (Ω). For any > 0 there are indeed examples of domains Ω and boundary data ⋅ where ∉ 3 2 + (Ω). In this sense, the vector potential ∈ 3 2 (Ω) introduced in this work has the highest possible regularity one can expect for arbitrary Lipschitz domains.
However, an interesting question remains. Suppose that the given data and ⋅ are such that the velocity field does happen to have higher regularity, say ∈ (Ω) for some ≫ 1 2 . McIntosh and Costabel have proven that in this case, another vector potential ∈ 1+ (Ω) exists. 29,Corollary 4.7 In other words, there always exists a stream function that is more regular than its velocity field by one order. In the numerical example discussed in section 7.4, such a smooth vector potential is given by: However, the numerical experiments indicate that the vector potential proposed in this work is not smooth. Therefore, the problem of devising an algorithm to approximate reliably and efficiently the "smoothest possible" stream function remains open.