An Exterior Neumann Boundary-Value Problem for the Div-Curl System and Applications

: We investigate a generalization of the equation curl (cid:126) w = (cid:126) g to an arbitrary number n of dimensions, which is based on the well-known Moisil–Teodorescu differential operator. Explicit solutions are derived for a particular problem in bounded domains of R n using classical operators from Clifford analysis. In the physically signiﬁcant case n = 3, two explicit solutions to the div-curl system in exterior domains of R 3 are obtained following different constructions of hyper-conjugate harmonic pairs. One of the constructions hinges on the use of a radial integral operator introduced recently in the literature. An exterior Neumann boundary-value problem is considered for the div-curl system. That system is conveniently reduced to a Neumann boundary-value problem for the Laplace equation in exterior domains. Some results on its uniqueness and regularity are derived. Finally, some applications to the construction of solutions of the inhomogeneous Lamé–Navier equation in bounded and unbounded domains are discussed.


Introduction
In this work, we use Clifford analysis and, for the case n = 3, quaternionic analysis in the study of inhomogeneous Moisil-Teodorescu systems defined on bounded or unbounded domains of R n , n ≥ 3.For higher dimensions, our main goal is to provide an explicit general solution to the n-dimensional generalization of the equation curl w = g under certain conditions imposed on the data g over bounded Lipschitz domains.Thereafter, restricting ourselves to the three-dimensional case, we will analyze the div-curl system without boundary conditions on different classes of exterior domains.When the normal component of the vector field is also known, one speaks of a Neumann problem; one of the first papers to address this boundary-value problem was [1] using the fundamental theorem of vector calculus.Later, the results were extended to exterior domains in [2] under the condition that | div w| and | curl w| decay faster that 1/|x| 2 as |x| → ∞.
Historically, one of the first works on the application of quaternionic analysis for elliptic systems in unbounded domains was the article [3].In that work, weighted Banach spaces L 2,α , W k,2,α were employed in order to guarantee a good behavior at infinity of the Teodorescu transform and the Cauchy operator, which are some classical operators in quaternionic analysis.Unfortunately, the main disadvantage of considering the Teodorescu transform on unbounded domains in the usual form (that is, using the same Cauchy kernel E n (x) = −x/(σ n |x| n ) as in the bounded case) is that boundedness is lost in the classical integrable spaces L p (Ω).A different approach was presented in the articles [4][5][6].This approach relies on the use of a perturbation of the Teodorescu transform through the addition of a monogenic term whose singularity lies outside of the unbounded domain under consideration.More precisely, the studies in those papers employed modified Cauchy kernels.Unlike those works, we will employ the usual kernel in the present manuscript.Moreover, in spite of the fact that the integrability of the integral operators is not improved, the operators exhibit a good asymptotic behavior at infinity for our purposes.
To that end, we will employ results from [4][5][6] with modified Cauchy kernels.Indeed, they are still valid for the present case, only the integrability ranges will be different.
One of the main results presented in this work (Theorem 5) establishes that a weak solution of the div-curl system in exterior domains star-shaped with respect to (w.r.t.) infinity has the form where is harmonic in Ω − .In fact, it is this property together with the asymptotic condition lim |x|→∞ |x| 2 grad ψ 0 (x) = 0 which allow us to construct a general solution with no boundary conditions.The non-uniqueness of the solution is clear.Indeed, if we add the gradient of a harmonic function in Ω − to Equation (1), then it will have the same divergence and rotational over all Ω − .Moreover, the solution Equation (1) admits the following Helmholtz-type decomposition in exterior domains (Corollary 3): Here, with ψ 0 defined in Equation (2).In addition, we obtained a second Helmholtz-type decomposition for the class of strong local Lipschitz exterior domains in terms of the layer potentials (Corollary 2): where v 0 is as in Equation ( 4), and α 0 satisfies α 0 (x) − 2 P.V.
Note that in the solenoidal part of both the Helmholtz-type decompositions Equations ( 3) and (6), the following operator appears: curl Moreover, it coincides with the Biot-Savart operator defined over Ω − .In particular, when curl w = j represents the current density, we obtain the Biot-Savart law in electromagnetism.This important Biot-Savart operator is also part of a strategic decomposition of the Teodorescu transform Equations ( 52) and (53), in which much of our Clifford analysis is based.More precisely, the operator Equation ( 9) coincides with the vector part of the Teodorescu transform T 2,Ω applied to curl w, that is, Interestingly, the non-uniqueness of these general solutions Equations ( 3) and ( 6) allows us to transform the associated Neumann BVP for the div-curl system into a Neuman BVP for the Laplacian, with results available in the literature on its existence, regularity and uniqueness.This gives rise to the main result of Section 6, Theorem 7.
It is worth pointing out that each of the operators appearing in the expression of the general solution Equation ( 1) is important an operator from quaternionic analysis.In the present work, we will frequently employ an important decomposition of the Teodorescu operator used in [7][8][9] for bounded domains of R 3 .In turn, the radial operator appearing in the last term of Equation ( 1) was defined in [10,11] as a generalization to exterior domains of an important family of radial operators.We will talk briefly about this operator at the beginning of Section 5.In this work, we will express the general solution following a quaternionic approach.More precisely, we will embed the div-curl system in the algebra of quaternions and then project this quaternionic solution to a purely vector one without affecting the system.
The outline of the paper is as follows.In Section 2, we present the notation with basic theory of Clifford analysis as well as some facts about the regularity of the domain.In Section 3, a general weak solution of the generalization to n-dimensions of the equation curl w = g is provided using an appropriate embedding argument to the Clifford algebra Cl 0,n .In Section 4, we present the construction of hyper-conjugate harmonic pairs in unbounded domains in terms of certain layer potentials and give explicit formulas for a solution of the div-curl system without boundary conditions for exterior domains satisfying the strong local Lipschitz condition.In Section 5, we derive another explicit solution of the div-curl system in exterior domains now under the geometric condition that Ω − is star-shaped w.r.t.infinity.The construction of this second solution relies on the properties of a family of radial integral operators restricted to a family of harmonic functions with good behavior at infinity.In Section 6, we analyze in detail the regularity and asymptotic behavior of the solution of the div-curl system Equation (1).Thereafter, we adjust this solution to construct a weak solution of the Neumann BVP of the div-curl system in exterior domains.Due to the easy handling of the radial operator that appears in the construction of the general solution (Equation (1)) of the exterior div-curl system, we will thoroughly analyze and adjust this general solution instead of the other general solution Equation (3) found in this work (see Theorem 3 for more details).
In Section 7, we found an equivalence between the solutions of the inhomogeneous Lamé-Navier Equation (113) in elasticity and the solutions of a inhomogeneous div-curl system (Lemma 1).Later, we apply the results obtained in the previous sections and provide a weak general solution of the inhomogeneous Lamé-Navier Equation (113).Moreover, we give explicit solutions in appropriate interior or exterior domain in R 3 .We close the section showing that these weak solutions are in fact strong solutions through an embedding argument.

Clifford Algebras
Throughout, we will let δ ij be the Kronecker delta function.Let us consider the real Clifford algebra Cl 0,n generated by the elements e 0 , e 1 , e 2 , • • • , e n and e 0 = 1, together with the relation e i e j + e j e i = −2δ ij , for i, j = 1, 2, . . ., n.Then, a basis for Cl 0,n is the set Define conjugation in Cl 0,n as ab = ba, ∀a, b ∈ Cl 0,n .Denote the elements in the real Clifford algebra as a = ∑ A a A e A ∈ Cl 0,n and define the following projections: As a consequence, an arbitrary element a ∈ Cl 0,n can be written as Define the scalar, non-scalar, vector, paravector and non-paravector parts of Let us embed the (n + 1)-dimensional Euclidean space R n+1 in Cl 0,n by identifying each vector x i e i .Sc x * = x 0 and Vec x * = x = ∑ n i=1 e i x i ∈ R n denote the scalar and the vector part of any arbitrary paravector x * .In the sequel, Ω = Ω + ⊂ R n will be a domain with a sufficiently smooth boundary ∂Ω, its exterior domain will be denoted by Ω − := R 3 \ Ω, and the elements x ∈ Ω ± will be called vectors.We will say that w : Here, the coordinates w A are real-valued functions defined in Ω ± .In particular, paravector-valued functions are denoted by g(x) = ∑ n i=0 g i (x)e i .For further information on Clifford analysis, we refer to the monographs in [12][13][14].Meanwhile, for quaternionic analysis we suggest the monographs in [15,16].
Recall that the Moisil-Teodorescu differential operator is defined as where ∂ i represents the partial derivative operator with respect to the variable x i .We say ) be the subspace of integrable functions defined in Ω ± and taking values in A, with the usual norm given by In particular, in this work we will let A = R, R 3 , H, Cl 0,n , with H ∼ = Cl 0,2 the algebra of quaternions.We define the Clifford version of the Teodorescu transform as (see [14]) where σ n is the surface area of the unit sphere in R n , w ∈ L p (Ω ± , Cl 0,n ) and the Cauchy kernel E n is given by In the bounded case, the function T Ω is defined as T Ω : L p (Ω, Cl 0,n ) → W 1,p (Ω, Cl 0,n ), for each 1 < p < ∞ (see [13,14]).Meanwhile, T Ω − : L p (Ω − , Cl 0,n ) → W 1,p (Ω − , Cl 0,n ) in the unbounded case, for each 3/2 < p < ∞ (see [4,6]).Moreover, T Ω ± is a right inverse of the Moisil-Teodorescu operator D in Ω ± , that is, DT Ω ± = I in Ω ± .Furthermore, T Ω ± is monogenic in Ω ∓ .Finally, if w ∈ L p (∂Ω), then the Cauchy operator is defined by where η is the outward normal vector to ∂Ω.

Geometric Properties of the Domain
It is well known that many properties of Sobolev spaces depend on the regularity of the domain.In the present study, we will require some classical Sobolev embedding theorems and the trace theorem for bounded and unbounded domains.Following [17] [Par.4.9], we will impose the following condition on the geometry of the domain.
Strong local Lipschitz condition: We say that Ω ± satisfies the strong local Lipschitz condition if there exist positive numbers δ and M, a locally finite open covering U j of ∂Ω, a real-valued function f j of n − 1 variables for each j, such that the following conditions are satisfied: (i) There exists some R ∈ N with the property so that every collection consisting of R + 1 of the sets U j has an empty intersection.(ii) Let Ω ± δ = {x ∈ Ω ± : dist(x, ∂Ω) < δ}.For every pair of points x, y ∈ Ω ± δ such that |x − y| < δ, there exists j such that x, y ∈ V j = {x ∈ U j : dist(x, ∂U j ) > δ}. ( (iii) Each function f j satisfies a Lipschitz condition with constant M.
(iv) For some Cartesian coordinate system For the bounded domain Ω = Ω + , the above requirements reduce to the simpler condition that Ω has a locally Lipschitz boundary.
In the following, we will always specify the regularity assumptions on ∂Ω which will be required to employ standard results from the theory of Sobolev spaces.Note that if Ω − satisfies the strong local Lipschitz condition, then the Sobolev embedding theorem [17] [Th.4.12, Part II] assures that W m,p (Ω − ) ⊂ C 0 (Ω − ), for mp > n.In particular, the strong local Lipschitz condition implies the cone condition (see [17] [Par.4.6]).Again by the Sobolev embedding theorem [17] for mp > n.Here, C 0 B (Ω − ) is the space of bounded and continuous functions in Ω − .In particular, if w is harmonic in Ω − and |w(x)| → 0 as |x| → ∞, then max In addition, if w ∈ W m,p (Ω − ), then the trace of w is well defined; moreover, for all mp > n, it follows that tr w = w| ∂Ω ∈ L p (∂Ω).In summary, in order to ensure the existence of the trace of functions in unbounded domains, we will work with strong local Lipschitz domains.

Clifford Integral Operators
Following the notation of the decomposition used in [7] for the case n = 3, we denote the component operators of the Teodorescu transform T Ω acting over a Clifford-valued function w = Sc w + NSc w = w 0 + NSc w as follows: Here, the component operators are given by Observe that the scalar product in the identity Equation ( 25) is the product of vectors in R 2 n .Additionally, recall that Sc(ab) = Sc(ab) = a • b, for all a, b ∈ Cl 0,n .On the other hand, if n = 3 and w = w 0 + w is a quaternion-valued function, then the integrand of T 2,Ω reduces to the cross product between E n and w [7,8].A similar decomposition of Equation (24) was also used in [9] for the perturbed Teodorescu transform, whose analysis allowed to give the explicit form of a right inverse of curl + λ, with λ ∈ C.
The next proposition is a direct consequence of differentiating under the integral sign and of the standard identity See [7] [Prop.3.2] for details of the proof in the particular case n = 3.

An N-Dimensional Generalization of the Div-Curl System
In this section, we are interested in the analysis of an n-dimensional inhomogeneous Moisil-Teodorescu system, whose component equations give rise to an n-dimensional generalization of the div-curl system where w = ∑ n i=1 w i e i is a paravector-valued function with a vanishing scalar part and g = ∑ A g A e A is a Cl 0,n -valued function.More precisely, we will assume that Sc g = g 0 and NSc g = ∑ |A|=2 g A e A = ∑ i,j g i,j e i e j when n > 3 and g = ∑ 3 i=0 g i e i is paravector-valued when n = 3.By applying the operator D in both sides of Equation ( 32), we can check that the condition Sc D[NSc g] = NPa NSc D[NSc g] = 0 is necessary for Equation (32) to have a solution.On the other hand, observe that w solves Equation (32) if and only if w is a solution of where g 0 = Sc g and g ij is the coefficient of e i e j in the expression of g.Note that the left-hand side of Equation ( 33) is the additive inverse of the divergence of w.Meanwhile, the lefthand side of Equation ( 34) is the generalization of the curl operator to n-dimensions, which was previously studied in [19].It is worth mentioning that the main goal in that article was to find necessary and sufficient conditions to obtain a unique solution w ∈ W 1,p 0 (Ω, R n ), which depends continuously on NSc g.
In the present work, we will follow a different approach.More precisely, our approach will hinge on solving the embedding of Equations ( 33) and ( 34) into the Clifford structure provided by the equivalent Equation (32).The steps of our construction method are the following:

•
Step 1. Find a Cl 0,n -valued solution to Equation (32) using the Teodorescu transform, that is, Step 2. Proceed then to describe the kernel of the non-paravector component operator NPa T 2,Ω [NSc g] and restrict the right-hand side of Equation ( 32) to this class.

•
Step 3. In turn, our Cl 0,n -valued solution becomes a paravector-valued solution after the restriction in Step 2. Afterwards, we use the theory of hyper-conjugate harmonic pairs in order to construct a paravector and monogenic function whose scalar part coincides with the scalar part of w given in Step 1.
Finally, taking the scalar part, the paravector part and the non-paravector part of the above-mentioned equation, we reach the conclusions.
We say that Sc w = w 0 and NSc w are hyper-conjugate harmonic pairs in Ω ± if w = w 0 + NSc w is monogenic in Ω ± .We will illustrate now a way to generate monogenic paravector-valued functions when only the scalar part is known.The idea is the same as that used for the three-dimensional case in [8] [Cor.A.2].It is worth pointing out here that this is not the only procedure.colorredSee, for instance, the radial integral operator used in reference [7] [Prop.2.3] for star-shaped bounded domains in R 3 .
The three-dimensional singular Cauchy integral operator satisfies both the identity S 2 ∂Ω = I and the Plemelj-Sokhotski formulas: Here, y ∈ Ω ± , Ω + = Ω and Ω − = R n \ Ω.The notation n.t.-means that the limit must be taken non-tangential.The scalar component operator of S ∂Ω acting over scalar valued functions will be denoted by K 0 , and it is of particular interest for the scope of this work.More precisely, This scalar operator is well known in harmonic analysis, and it is fundamental in the classical Dirichlet problem.Moreover, whenever Ω be a bounded Lipschitz domain, then , where the value of (Ω) depends only of the Lipschitz character of ∂Ω.Returning to the construction of hyper-conjugate harmonic pairs, a natural way to construct them is through the Cauchy operator Equation ( 19), which generates monogenic functions.More precisely, , where (I + K 0 )ϕ 0 = tr w 0 is monogenic in Ω and Sc w = w 0 .In other words, w 0 and w = Vec F ∂Ω [2ϕ 0 ] are hyper-conjugate harmonic pairs.Proof.Observe that the Plemelj-Sokhotski formula Equation ( 36) describes the trace of the Cauchy operator.As a consequence tr By the maximum principle for harmonic functions, we conclude that Sc F ∂Ω [2ϕ 0 ] = w 0 in Ω, as desired.
As mentioned at the beginning of the present section, the necessary conditions for the equivalent system Equation (32) to have a solution coincide with the first two hypotheses on Ω in the following result.Meanwhile, the third boundary hypothesis imposed on NSc g is used to ensure that the solution has vanishing non-paravector part.
and NPa NSc [ηNSc g] = 0 on ∂Ω, then a weak solution w of the n-dimensional div-curl system in Equations ( 33) and (34) is given by where This solution is unique up to the gradient of a scalar harmonic function in Ω.Moreover, in the case n > 3, we have w Proof.In the proof, we will follow Steps 1, 2 and 3 described above.Using Gauss' theorem on Cl 0,n (see [14] [Rmk.A.2.23]), we obtain By hypothesis, NPa NSc D[NSc g] = 0 in Ω and NPa NSc [ηNSc g] = 0 on ∂Ω.It follows that NPa T 2,Ω [NSc g] = 0 in Ω, as desired.On the other hand, Proposition 2 guarantees that T 0,Ω [NSc g] is a scalar harmonic function in Ω.In turn, Proposition 3 implies that F ∂Ω [2ϕ 0 ] is monogenic and its scalar part is T 0,Ω [NSc g], where Equation ( 40) , and note that it is purely vectorial by virtue of NPa T 2,Ω [NPa g] ≡ 0. As a consequence, Equation (39) is reached.Moreover, This means that w satisfies the equivalent system, Equation (32).The fact that w belongs to the Sobolev space W 1,p (Ω, R 3 ) is a direct consequence of the properties of the Teodorescu and Cauchy operators.Finally, if n > 3, then [NSc g] 1 ≡ 0. In turn, this identity implies that E n (x − y) • NSc g = 0, which means that T 0,Ω vanishes in Ω.

Corollary 1.
Let Ω and NSc g be as in Theorem 2.Then, a right inverse of the n-dimensional generalized curl operator is Proof.Taking g 0 = 0 in Equation (39), we readily obtain the expression for the right inverse of the generalized curl operator.To see if R Ω + [NSc g] is divergence-free, we will use the alternative expression which is what we wanted to prove.
Before closing this section, we must point out that Theorem 2 and Corollary 1 generalize to n dimensions those results recorded as [8] [Th.A.1] and [8] [Cor.A.3], respectively, which were valid for bounded Lipschitz domains in R 3 .On the other hand, as illustrated in the last part of Theorem 2, this construction is mathematically much more interesting in the three-dimensional case, needless to say that this is the physically most relevant.In fact, if n = 3, then the assumptions on NSc g = g in Theorem 2 and Corollary 1 that involve the non-paravector part disappear, while the hypothesis Sc D[NSc g] = 0 becomes the irrevocable condition that g has zero divergence.Finally, we point out that the present work is not the first to make use of Clifford analysis and the construction of hyper-conjugate harmonic pairs to address inhomogeneous Moisil-Teodorescu systems.A recent work in which these tools were employed is [21].

Unbounded Domains
From now on, we will restrict our attention to the case n = 3, and analyze the classical div-curl system in unbounded domains Ω − = R 3 \ Ω.To that end, we require some hypothesis on ∂Ω to guarantee that the operator I − K 0 is invertible in L p (∂Ω, R).
In Section 3, we use the fact that the operator I + K 0 is invertible in L p (∂Ω) for all 2 − (Ω) < p < ∞, when ∂Ω is Lipschitz.We are interested now in the inversion of the operator I − K 0 in L p (∂Ω).To analyze in more detail the range of p for which this operator is invertible, let us define the boundary averaging operator A as Note also that K 0 does not interfere with the averaging process: This and [8] [Prop.3.3]) show that the operator I − K 0 sends L p 0 (∂Ω, R) into itself.Moreover, this operator has an inverse with Lipschitz and γ > 0.
The following result provides an alternative form to complete a scalar-valued harmonic function to a paravector-valued monogenic function in the exterior domain Ω − .

Proof.
Mimicking the analysis at the end of Section 2.2, we can assure tr w 0 ∈ L p (∂Ω), for p > 3/m.By using the Plemelj-Sokhotski formula in Ω − Equation (36), we obtain that is, tr where the last inequality comes from the continuity of the operator (I Using the asymptotic behavior of w 0 (x) = o(1) as |x| → ∞ and Equation (50), it follows that as |x| → ∞.Due to the the uniqueness of the Dirichlet problem in exterior domains, we readily conclude that Sc F ∂Ω [2 comes from Equation (50) and Cauchy's integral formula is used for exterior domains [14] [Th.7.14].

Teodorescu Transform over Unbounded Domains
To start with, note that the Teodorescu transform defined previously in Equation ( 17) reduces in the three-dimensional case to Furthermore, its decomposition is simplified to the expression In turn, the operators T 0,Ω − , T 1,Ω − , T 1,Ω − , defined previously in the n-dimensional case as Equations ( 25)- (27), respectively, are reduced to the following expressions, respectively: Note that the last identities in the above expressions of T 0,Ω − , T 1,Ω − and T 2,Ω − are derived from Proposition 1. Observe that T 0,Ω − , T 1,Ω − , T 2,Ω − are the divergence, gradient and curl of the Newton potential L Ω − [ϕ](x) = Ω − (ϕ(y)/|x − y|) dy, respectively.This potential operator (also known as volume potential) has been extensively studied in various works, such as [22] [Sec.2.2] and [23].In addition to the role that these component operators have in our construction of solutions for Equation (57), they also provide a lot of analytical information.For instance, T 2,Ω and T 2,Ω − are the Biot-Savart operators for bounded and unbounded domains, respectively.
We recall next some properties of the Teodorescu transform in classical L p spaces.As mentioned above, one of the disadvantages of using a kernel without any modification is that its integrability range is reduced.Indeed, if > 0 and Integrating over Ω yields which is finite for p > 3/2.Let w ∈ L p (Ω − ).By utilizing the Fubini-Tonelli theorem, we obtain From the fact that the kernel E 3 (x − •) and w belong to L p (Ω − ) for 3/2 < p < ∞, we readily obtain T Ω − [w] ∈ L p (Ω − ).

The Div-Curl System over Unbounded Domains
In this stage, we will give an explicit solution for the Equation (57) on unbounded domains of R 3 satisfying the strong Lipschitz condition with weaker topological constraints.To that end, we will recall some auxiliary results reported in [7,8,24].Fortunately, the operator theory needed for the quaternionic integral operators over unbounded domains is already well developed [4][5][6].The novelty now lies in the use of the monogenic completion method discussed in Proposition 4 via the single layer operators.
Let us consider the div-curl system without boundary conditions where ) and div g = 0 in Ω − .Note that the equivalence of the systems in Equations ( 32)-( 34) is readily verified when n = 3.Moreover, the system in Equation ( 57) is equivalent to Equation ( 58) Due to the action of the operator D to a vector-valued function, w is rewritten in quaternionic notation as D w = − div w + curl w.In the same way as for the bounded case, the mentioned equivalence will be the key in the analysis of the exterior div-curl system.
By Equations ( 61) and (62), we can conclude that T 0,Ω − [ g] ∈ W 2,p (Ω − , R).As a consequence, by Proposition 4, Using this and the decomposition Equation (52), we obtain is a purely vector solution of Equation ( 58), whence the conclusion of this result follows.
Note that Equation (59) can be rewritten as where . Define the single layer potential [25] as It is worth pointing out that the Cauchy operator evaluated in scalar functions possesses a decomposition in terms of the operators div and curl [8].More precisely, Using the last equation with ϕ 0 = 2α 0 and (I − K 0 )α 0 = tr T 0,Ω − [ g] as above, and replacing the second and third expressions of Equation (53) in Equation (59), we observe that the solution of the div-curl system can be rewritten in a way similar to the classic Helmholtz decomposition theorem.More precisely, we have the following result.
Corollary 2. Under the same hypothesis of Theorem 3, the solution, Equation (59), admits a Helmholtz-type decomposition of the form and Comparing the decomposition Equation (69) with the classical Helmholtz decomposition on the entire three-dimensional space [26] [p.166] [27] [Lem.3.1, 3.2], we readily observe that they adopt similar forms.Note that the vector field v is divergence-free.This follows from div which is a consequence of the first equation of Equation ( 53) and the proof of Theorem 3.

Div-Curl System in Exterior Domains
In this section, we derive another explicit solution to the div-curl system Equation (57), this time using another method to generate hyper-conjugate harmonic pairs.The cornerstone now is a radial integral operator defined on an infinite ray instead of the integral equation method provided by the layer potentials.
For the remainder of this manuscript and for the sake of convenience, we will suppose that Ω is star-shaped w.r.t. the origin.It is worth pointing out that if Ω is star-shaped w.r.t.any other point, then a simple translation would make it star-shaped w.r.t. the origin.The radial integral operator mentioned above was recently proposed and firstly analyzed in [10,11].There exists an important family of radial integral operators in star-shaped domains, which takes on the form where usually α > −1.Using standard relations such as ∂w 0 (tx)/∂t = x • grad w 0 (tx), one may readily verify the following relations: div I α = I α+1 div; grad . This family of operators plays an important role in the theory of special functions as well as in mathematical physics.
Another interesting application appears in quaternionic analysis when α = 0. Indeed, this radial integral operator generates harmonic functions for each u 0 which is a harmonic function in This means that the radial operator U Ω provides an explicit way to generate hyper-conjugate harmonic pairs in the star-shaped domain with respect to the origin when the scalar part is known.For convenience, we recall next the main result of [7].
If div g = 0 in Ω, and if g = g 0 + g ∈ L p (Ω) for 1 < p < ∞, then a general weak solution of the div-curl system is given by Moreover, this solution is unique up to the gradient of a harmonic function in Ω.
We now turn our attention to unbounded domains.To start with, note that a similar radial integral operator J α acting on functions defined on Ω − = R 3 \ Ω was defined in [10,11], for star-shaped domains Ω. Equivalently, Ω − will be star-shaped w.r.t.infinity, which means that any infinite ray emanating from x ∈ Ω − is entirely contained in Ω − .In other words, Ω − is star-shaped w.r.t.infinity if λx ∈ Ω − for all x ∈ Ω − and λ > 1.More precisely, the following operator preserves the above-mentioned properties of the operator I α when it is restricted to a class of functions with a suitable behavior at infinity: Define the class Note that this class of scalar-valued harmonic functions is non-empty in that E 3 (x) = −x/|x| 3 belongs to A ∞ (Ω − ) component-wise.The harmonicity in Ω − is straightforward in that it is monogenic in R 3 \ {0}.Meanwhile, the radiation condition at infinity readily follows.Indeed, observe that grad Let us define We will call U Ω − the exterior monogenic completion operator in light of the next result.
Proof.Beforehand, note that w 0 + U Ω − [w 0 ] satisfies Equation (80) if and only if w 0 and Using the fact that w 0 is harmonic in Ω − and some identities from the vector calculus, it follows that the following is satisfied for The action of the rotational operator to the integrand of However, lim t→∞ t 2 |x| 2 grad w 0 (tx) = 0 by that hypothesis.We conclude finally that Equation (84) reduces to − grad w 0 (x), as desired.
We must mention that the operator J α played a fundamental role in [10] [Th.2] to obtain the general solution of the biharmonic equation.It was also crucial to obtain the general solution of the div-curl system in exterior domains when the known data g 0 and g belong to the class of functions A ∞ (Ω − ) component-wise [10] [Th. 3].Our next result is more general in that we consider arbitrary integrable functions in Ω − and not only harmonic functions in the class A ∞ (Ω − ).
Theorem 5. Suppose that Ω is a bounded domain and Ω − is star-shaped w.r.t.infinity.Let 3/2 < p < ∞, g 0 + g ∈ L p (Ω − , H) and div g = 0 in Ω − .Then, a weak solution w of the div-curl system Equation (57) in Ω − is given by ( This general solution is unique up to the gradient of a harmonic function in Ω − .
Proof.The proof is similar to that of Theorem 3; only the generation of the monogenic function whose scalar part coincides with the operator T 0,Ω − [ g] changes.Consider it only remains to verify that the hypothesis of Proposition 5 holds.In other words, we will show that T 0,Ω − [ g] belongs to the family of functions in A ∞ (Ω − ): Using the estimation Equation (87) and letting R → ∞, we obtain grad which means that the harmonic function Finally, since this monogenic function has the same scalar part as the quaternionic solution satisfies the equivalent system Equation (58), which is what we wanted to establish.
In our derivation of the solution to the exterior div-curl problem Equation (57), we followed a path different from the classical works by Girault and Raviart [28].The present solution hinges on the exterior monogenic completion operator U Ω − defined in Equation (79) (which was firstly introduced in [10,11]), and on the properties derived in the present work for the component operators of Corollary 3.Under the hypothesis of Theorem 5, the solution Equation (85) admits a Helmholtztype decomposition where Moreover, div v * is harmonic in Ω − .
Proof.By Equation (53), we obtain Meanwhile, a simple computation shows that which yields Equation (91).On the other hand, due to T 0,Ω − [ g] being harmonic in Ω − .It only remains to prove that the second term of Equation ( 97) is also harmonic in the exterior domain, but this fact is derived from the fact that Unfortunately, unlike the Helmholtz-type decomposition given in Equation (69), the new decomposition, Equation (69), is not divergence-free in the exterior domain Ω − .Later, in Theorem 7, the regularity of the solution Equation (85) will be analyzed as well as its asymptotic behavior.

Neumann Boundary-Value Problems
In this stage, we will analyze an exterior Neumann boundary-value problem associated with the div-curl system Equation (57).More precisely, we will check that there exists a Helmholtz-type solution of the boundary-value problem which preserves the optimal behavior at infinity whenever g = g 0 + g belongs to L p (Ω − ).
Firstly, note that the normal trace of the solution Equation ( 85) is well defined.As a consequence, solving for the following exterior Neumann boundary-value problem gives: which is equivalent to solving the Neumann boundary-value problem for the Laplace equation in exterior domains Here, g = g 0 + g and is the general solution provided by Theorem 5.More precisely, w = H Ω − [g] + grad u 0 solves Equation (98) if and only if u 0 solves Equation (99).It is the non-uniqueness of the solution of the div-curl system without boundary conditions and the fact that the normal trace of H Ω − [g] is well defined that allow us to formulate the equivalent Neumann problem Equation (99).
Theorem 6 (Neudert and von Wahl [27] [Th.2.1]).Let Ω ⊂ R 3 be a bounded domain with a smooth boundary, let Ω − = R 3 \ Ω and assume ψ ∈ C 0 (∂Ω, R).Then, the Neumann boundaryvalue problem Before introducing the main theorem of this section, we will establish some crucial results.To start with, we will prove next that the composition of the exterior monogenic completion operator U Ω − with T 0,Ω − [ g] preserves the regularity and asymptotic behavior of the Teodorescu transform T Ω − .Proposition 6.Let Ω be a bounded domain and Ω − be star-shaped w.r.t.infinity.Let 3/2 < p < ∞ and g ∈ L p (Ω − , H).
. Without loss of generality, we will only calculate the asymptotic behavior of grad(x × grad T 0,Ω − [ g](x)) 1 .
Taking the e 1 component of Equation ( 102), note that grad x x 2 x 3 |x − y| 7 dy.

Theorem 7.
Let Ω be a bounded domain with smooth boundary, let Ω − be star-shaped w.r.t.infinity, and suppose that Ω − satisfies the strong local Lipschitz condition.Let 3 < p < ∞, g = g 0 + g ∈ L p (Ω − , H), div g = 0 in Ω − and ϕ 0 ∈ C 0 (∂Ω).Then, the exterior Neumann boundary-value problem Equation (98) has a unique solution where ) is the general solution described by Theorem 5 and, in turn, Proof.The first part of the proof is reduced to verifying that H Ω − [g]| ∂Ω • η ∈ C 0 (∂Ω) in light of Theorem 6 and the equivalence between the systems of Equations ( 98) and (99).Meanwhile, the uniqueness of solutions of Equation ( 98) is derived from the uniqueness of solutions of Equation ( 99).From Proposition 6, ).In turn, the Sobolev Imbedding Theorem [17] [Th.4.12, Part II] assures that W Note that grad u 0 (x) = O(|x| −2 ) as |x| → ∞ follows from Theorem 6.The fact that ) satisfies the same decay condition at infinity results from Proposition 6 (when we obtained that as |x| → ∞) and from the asymptotic behavior of the Teodorescu transform which reads To establish the last part of the proof, note that Vec is divergence-free in Ω in that the Teodorescu transform T Ω − is well defined over all R 3 and it is monogenic in Ω.From Equation (82), we have By [27] [Lem.2.2], it follows that u 0 (x) = O(|x| −2 ) as |x| → ∞, as needed.
In [27], [Th.3.2] was given an exhaustively classification of the asymptotic behavior of the solutions of the Neumann BVP Equation ( 98) under an appropriate functional setting.In that work, the authors used the solutions of the div-curl system in the entire threedimensional space and correct the boundary values by harmonic vector fields; the second part is similar to the equivalent BVP Equation (99) considered in this work.
Regularity of the solution: We can apply the Sobolev embedding theorem so that, if we require a higher regularity for the function g = g 0 + g, then the range in which the embedding is achieved is improved Observe that up to this step, we have still not used the extra geometric condition of the domain.If Ω − satisfies the strong local Lipschitz condition, then W 2,p (Ω − ) ⊂ C 0 (Ω − ) for p > 3/2 follows from [17] [Th.4.12, Part II].Therefore, the range of integrability in the hypotheses of Theorem 7 can be improved from where H Obviously, we can modify the functional framework of our Neumann boundaryvalue problem for the div-curl system in the context of weighted Sobolev spaces which, as shown in [29], gives a correct functional setting to the exterior Neumann problem for the Laplace equation.This approach was also used to analyzed the regularity of the Teodorescu transform in exterior domains [3].

Right Inverse of the Curl and Double Curl Operator
Let us define the subspace of divergence-free L p -functions Sol p (Ω − ) as the set Sol p (Ω − ) := { u ∈ L p (Ω − ) : div u = 0 in Ω − }.Taking g 0 ≡ 0 in the general solution Equation (85), we readily obtain a right inverse operator for the curl operator in exterior domains whose complement is a star-shaped domain, namely, Meanwhile, when Ω − satisfies the strong local Lipschitz condition, the operator is also a right inverse for the curl operator in Sol p (Ω − ) when (I Moreover, both operators are divergence-free invariant in Ω − .Due to the Helmholtz-type decomposition Equation (91), a right inverse operator of the curl curl operator in exterior domains of star-shaped domains is given by Similarly, we can obtain a right inverse operator for curl curl by taking g 0 = 0 in the Helmholtz-type decomposition Equation (69).Indeed, if (I is also a right inverse of the curl curl operator in exterior domains satisfying the strong local Lipschitz condition.As a consequence of this discussion, given g ∈ Sol p (Ω − ), there exists S 2,Ω −

Lamé-Navier Equation
In this section, we will apply the results obtained in the previous sections in order to provide an explicit solution to the well known Lamé-Navier problem in elasticity [30].Let us consider the inhomogeneous Lamé-Navier equation where λ and µ are known as the first and second Lamé parameters, respectively.The parameters on the right-hand side of Equation ( 113) have physical significance.For example, T 0 denotes the temperature field, f represents the body forces, and the residual strain res ij defines the vector field E res as follows: It is worth recalling that this system with the right-hand side of Equation (113) equal to zero was originally introduced by G. Lamé while he was studying the method of separation of variables for solving the wave equation in elliptic coordinates [31].Recently, several works have addressed the homogeneous Lamé-Navier equation using the tools of quaternionic analysis.For instance, [32,33].
Quaternion algebra was also used in [34] to give an extension of the classical Kolosov-Muskhelishvili formulas from elasticity to three dimension.This approach is also based on the classical harmonic potential representation due to Papkovich and Neuber, as well as on a monogenic representation.In the latter technique, the main tool is the decomposition of harmonic functions as the sum of a monogenic with an antimonogenic function in the quaternion setting.For the complete details of this decomposition, see [35].Here, we will proceed following a completely different path.In fact, we will show that the solutions of Equation ( 113) can be constructed by solving a specific div-curl system whose solutions are readily at hand with the theory developed herein.Lemma 1.Let Ω ± be a bounded or unbounded domain in R 3 .Let u ± satisfy the system with f ± = f ± 0 + f ± being a quaternionic solution of the inhomogeneous Moisil-Teodorescu system D f ± = G in Ω ± , and Then, u ± is a solution of the inhomogeneous Lamé-Navier Equation (113) in Ω ± , respectively.
Let f ± = f ± 0 + f ± be a quaternionic solution of the Moisil-Teodorescu system D f ± = G in Ω ± , respectively.Moreover, f ± = T Ω ± [ G] + H ± , where H ± is an arbitrary monogenic function in Ω ± .Equating the scalar parts of D f ± = G, we readily obtain div f ± = 0 in Ω ± , respectively.Now, equating the vector parts of D f ± = G, we have The conclusion of the result follows now by Equation (115).
Proof.By construction, div f = div F = 0 and div h = div H = 0 in Ω − , with F and H being the vector parts of the arbitrary monogenic functions in Ω and Ω − , respectively.It only remains to verify that f and h belong to L p (Ω + ) and L p (Ω − ), for 1 < p < ∞ and 3/2 < q < ∞, respectively.Since T Ω ± : W m,p (Ω ± ) → W m+1,p (Ω ± ) for m ≥ 0, then f ∈ W 1,p (Ω) and h ∈ W 1,p (Ω − ) for 1 < p < ∞ and 3/2 < p < ∞, respectively.The result readily follows from Lemma 1 and Theorem 4 for the star-shaped domain Ω + , and from Lemma 1 and Theorem 5 for the exterior domain Ω − .For the unbounded case scenario, the regularity of the solution comes from Corollary 5, due to the fact that u − = H Ω − [− f 0 /(λ + 2µ) − f /µ].In the bounded case, the proof of Proposition 6 yields that U Ω + = U Ω defined in Equation (74) belongs to W 2,p (Ω), as desired.where Here, f = f 0 + f = T Ω [ G] + F and h = h 0 + h = T Ω − [ G] + H are quaternionic solutions of the inhomogeneous Moisil-Teodorescu systems D f = G in Ω and Dh = G in Ω − , and F and H are arbitrary monogenic functions in L p (Ω ± ), respectively.These general solutions are unique up to the gradient of a harmonic function in Ω ± , respectively.