A new family of boundary-domain integral equations for the diffusion equation with variable coefficient in unbounded domains

A system of Boundary-Domain Integral Equations is derived from the mixed (Dirichlet-Neumann) boundary value problem for the diffusion equation in inhomogeneous media defined on an unbounded domain. This paper extends the work introduced in [ 25 ] to unbounded domains. Mapping properties of parametrix-based potentials on weighted Sobolev spaces are analysed. Equivalence between the original boundary value problem and the system of BDIEs is shown. Uniqueness of solution of the BDIEs is proved using Fredholm Alternative and compactness arguments adapted to weigthed Sobolev spaces.


1.
Introduction. Boundary Domain Integral Equations appear naturally when applying the Boundary Integral Method to boundary value problems with variable coefficient. This class of boundary value problems has a wide range of applications in Physics or Engineering, such as, heat transfer in non-homogeneous media [27], motion of laminar fluids with variable viscosity [5], or even in the acoustic scattering by inhomogeneous anisotropic obstacle [6].
The popularity of the Boundary Integral Method is due to the reduction of the discretisation domain. For example, if the boundary value problem (BVP) is defined on a three dimensional domain, then, the boundary integral method reduces the BVP to an equivalent system of boundary integral equations (BIEs) defined only on the boundary of the domain. However, this requires an explicit fundamental solution of the partial differential equation appearing in the BVP. Although these fundamental solutions may exist, they might not always be available explicitly for PDEs with variable coefficients. To overcome this obstacle, one can construct a parametrix using the known fundamental solution. A discussion on fundamental solution existence theorems, algorithms for constructing fundamental solutions and parametrices is available in [24]; for classical examples of derivation of Boundary Domain Integral Equations refer to [7] for the diffusion equation with variable coefficient in bounded domains in R 3 ; [25] for the same problem applying a different parametrix; [26] for the Dirichlet problem in R 2 and [22] for the mixed problem for the compressible Stokes system, as an example of derivation of BDIEs from a PDE system.
The introduction of a parametrix for BVPs with variable coefficient leads to a system of integral equations not only defined on the boundary but also in the domain. Still, one can transform domain integrals into boundary integrals applying the methods shown in [1]. These methods help to preserve the reduction of dimension while also remove singularities appearing in the domain integrals.
The approximation of numerical solutions of BDIEs is a relevant problem nowadays. In particular, the very recent article [2] focuses on the solution of the analogous mixed BVP presented in this paper in R 2 . In [3], the authors show that it is possible to obtain linear convergence with respect to the number of quadrature curves, and in some cases, exponential convergence. Analogous research in 3D shows the successful implementation of fast algorithms to obtain the solution of boundary domain integral equations, see [27,14,28]. Furthermore, the authors [4] show the application of the Boundary Domain Integral Equation Method to the study of inverse problems with variable coefficients.
A parametrix is not unique, see discussion on [25,Section 1]. The study of different parametrices is adventageous to construct parametrices for PDE systems. Moreover, numerical methods may work with one parametrix more efficiently than with another. However, before attempting numerical experiments, results on the existence and uniqueness of solution need to be established and that is the purpose of this paper.
In this paper, we extend the results presented in [25] to unbounded domains which employ a different parametrix from the one used in [9].
In unbounded domains, the mixed problem is set in weighted Sobolev spaces to allow constant functions in unbounded domains to be possible solutions of the problem. Hence, all the mapping properties of the parametrix based potential operators are shown in weighted Sobolev spaces.
An analysis of the uniqueness of the BDIES is performed by studying the Fredholm properties of the matrix operator which defines the system. Unlike for the case of bounded domains, the Rellich compactness embeding theorem, see [18,Theorem 3.27], is not available for Sobolev spaces defined over unbounded domains. Nevertheless, we present a lemma to reduce the remainder operator to two operators: one invertible and one compact. Therefore, we can still benefit from the Fredholm Alternative theory to prove uniqueness of the solution.

Weighted Sobolev spaces.
Let Ω = Ω + be an unbounded exterior connected domain of R 3 . Let Ω − := R 3 Ω + the complementary (bounded) subset of Ω. The boundary S := ∂Ω is simply connected, closed and infinitely differentiable, S ∈ C ∞ . Furthermore, S := S N ∪ S D where both S N and S D are non-empty, connected disjoint submanifolds of S. The border of these two submanifolds is also infinitely differentiable: With regards to function spaces that we employ in this paper, D(Ω) := C ∞ comp (Ω) denotes the space of test functions, C ∞ functions with compact support inside of Ω. The space D * (Ω) denotes the space of distributions or generalised functions. The space D(Ω) is the set of restrictions to Ω of functions from D(R 3 ). We also use Sobolev spaces H k (Ω) with k ∈ Z; Bessel potential spaces on the boundary of the domain H s (∂Ω), where s ∈ R (see e.g. [18,16]  To ensure uniquely solvability of the BVPs in exterior domains, we will use weighted Sobolev spaces with weight ω(x) = (1 + |x| 2 ) 1/2 , (see e.g., [9]). Let be the weighted Lebesgue space and H 1 (Ω) the following weighted Sobolev (Beppo-Levi) space constructed using the L 2 (ω −1 ; Ω) space endowed with the corresponding norm . Taking into account that D(Ω) is dense in H 1 (Ω) it is easy to prove that D(Ω) is dense in H 1 (Ω). For further details, cf. [9, p.3] and more references therein.
If Ω is unbounded, then the seminorm If Ω is a bounded subdomain of an unbounded domain Ω and g ∈ H 1 (Ω), then g ∈ H 1 (Ω ).
For any generalised function g in H −1 (Ω), we have the following representation property, see [9, Section 2], g = 3 i=1 ∂ i g i + g 0 , g i ∈ L 2 (R 3 ) and are zero outside the domain Ω, whereas g 0 ∈ L 2 (ω; Ω). Consequently, D(Ω) is dense in H −1 (Ω) and 3. Traces, conormal derivatives and Green identities. We consider the following differential operator where a(x) ∈ C 2 , a(x) > 0, is a variable coefficient. It is easy to see that if a ≡ 1 then, the operator A becomes the Laplace operator ∆.
Here and thereafter, we will assume the following condition on the coefficient a(x).
Condition 3.1. The coefficient a(x) belongs to the space L ∞ (Ω). Furthermore, there exist two positive constants, C 1 and C 2 , such that: (3.2) Condition 3.1 is necessary so that the operator A acting on u ∈ H 1 (Ω) is well defined in the weak sense. Hence, we define the operator A in the weak sense as where For a scalar function u ∈ H 1 (Ω) in virtue of the trace theorem it follows that γ ± u ∈ H 1/2 (S) where the trace operators from Ω ± to S are denoted by γ ± respectively. Consequently, if u ∈ H 1 (Ω), then u ∈ H 1 (Ω) and it follows that γ ± u ∈ H 1/2 (S), (see, e.g., [18,19]). For u ∈ H s (Ω); s > 3/2, we can define by T ± the conormal derivative operator acting on S understood in the classical sense: where n(x) is the exterior unit normal vector to the domain Ω at a point x ∈ S. However, for u ∈ H 1 (Ω) (as well as for u ∈ H 1 (Ω)), the classical co-normal derivative operator may not exist in the trace sense. This issue is overcome by introducing the following function space for the operator A, (cf. [9]) . Now, if a distribution u ∈ H 1,0 (Ω; A), we can appropriately define the conormal derivative T + u ∈ H −1/2 (S) using the Green's formula, cf. [18,9], where γ + −1 : H 1/2 (S) → H 1 (Ω) is a continuous right inverse to the trace operator γ + : H 1 (Ω) −→ H 1/2 (S) while the brackets u, v S represent the duality brackets of the spaces H 1/2 (S) and H −1/2 (S) which coincide with the scalar product in L 2 (S) when u, v ∈ L 2 (S).
The operator T + : H 1,0 (Ω; A) −→ H −1/2 (S) is bounded and gives a continuous extension on H 1,0 (Ω; A) of the classical co-normal derivative operator (3.5). We remark that when a ≡ 1, the operator T + becomes the continuous extension on H 1,0 (Ω; ∆) of the classical normal derivative operator T + ∆ u = ∂ n u := n · ∇u. In a similar manner as in the proof [18,Lemma 4.3] or [10, Lemma 3.2], the first Green identity holds for a distribution u ∈ H 1,0 (Ω; A) Applying the identity (3.8) to u, v ∈ H 1,0 (Ω; A), exchanging roles of u and v, and then subtracting the one from the other, we arrive to the following second Green identity, see e.g. [18] 4. Boundary value problem. Now that we have shown that if u ∈ H 1,0 (Ω; A), then its trace and its conormal derivative are well defined, it is possible to formulate the mixed problem for the operator A for which we aim to derive an equivalent of system of boundary-domain integral equations (BDIEs).
Mixed problem. Find u ∈ H 1,0 (Ω; A) such that The previous BVP can be represented with the following operator equation The following result is well known and it has been proven [9, Appendix A] by using variational settings and the Lax Milgram lemma.
. It is clear that hypotheses of the Theorem 4.1 are satisfied under the assumption of Condition 3.1. Hence, the mixed BVP problem (4.1)-(4.3) is uniquely solvable.

5.
Parametrices and remainders. We define a parametrix (Levi function) P (x, y) for the differential operator A differentiating with respect to x, as a function on two variables that satisfies where δ(.) is the Dirac distribution and the term R(x, y) is a weakly singular distribution, i.e. O(|x − y| −2 ), so-called remainder. A given operator A may have more than one parametrix. For example, the parametrix was employed in [21,7], for the operator A, given in (3.1), where is the fundamental solution of the Laplace operator. The remainder corresponding to the parametrix P y is In this paper, we consider the parametrix P x used in [25,23], where analogous results to the ones presented in the upcoming sections have been obtained in bounded domains with smooth and Lipschitz boundary. The parametrix P x is defined as follows: which leads to the corresponding remainder Due to the smoothness of the variable coefficient a(x), both remainders R x and R y are weakly singular, i.e., which has been used to derive analogous results to those in this paper, in [9]. The parametrix P y has been widely analysed in the literature, see [21,20,14,7,8]. The difference between both parametrices relies on the dependence from the variable of the coefficient a(x) or a(y). Clearly, choosing a parametrix involving a(y) simplifies the expression of the remainder as the coefficient a(y) acts as a constant when differentiating with respect to x which is the variable of differentiation of the operator A. However, for some PDE problems, it is not always possible to obtain a parametrix that depends exclusively on a(y) and not on a(x). This is the case of the Stokes system, see [22]. Hence, the usefulness of the analysis of the family of parametrices depending on a(x).
6. Volume and surface potentials. Boundary-domain integral equations are usually formulated in terms of parametrix-based surface and volume potential operators. In this section, the surface and volume potentials based on the parametrix P x are introduced. We analyse their mapping properties in weighted Sobolev spaces. Additional boundedness conditions are often imposed on the variable coefficient a(x) in order to prove the boundedness properties of the potential operators.
Condition 6.1. We will assume the following condition from now on unless stated otherwise: 2) and (6.1), then where the constants k 1 and k 2 do not depend on g ∈ H 1 (Ω). This implies that the functions a and 1/a behave now as multipliers in the space H 1 (Ω). Furthermore, as long as a ∈ C 1 (S), then ∂ n a is also a multiplier.
The volume parametrix-based Newton-type potential and the remainder potential are respectively defined, for y ∈ R 3 , as

BOUNDARY-DOMAIN INTEGRAL EQUATIONS IN UNBOUNDED DOMAINS 5103
The parametrix-based single layer and double layer surface potentials are defined for y ∈ R 3 : y / ∈ S, as We also define the following pseudo-differential operators associated with direct values of the single and double layer potentials and with their conormal derivatives, for y ∈ S, The operators P, R, V, W, V, W, W and L can be expressed in terms the volume and surface potentials and operators associated with the Laplace operator, as follows The symbols with the subscript ∆ denote the analogous operators for the constant coefficient case, a ≡ 1. Furthermore, by the Lyapunov-Tauber theorem (cf. [15,13] and more references therein), 3), (7.5), and (7.6), follow from the parametrix relation with the fundamental solution (5.2) A similar argument would apply for the operators V and V as they all share the same integral kernel. The remainder relation (7.4) follows from expanding the expression (5.3) by applying the product rule.
To obtain relation (7.7) and (7.8), we need to apply the product rule to the kernel T + P (x, y) Multiplying by ρ and integrating over S, we see that the first term in (7.12) leads to harmonic single layer potential term in (7.2) and the second term coincides with the harmonic double layer potential term. Relations (7.9) and (7.10) directly follow from applying the conormal derivative operator T + at both sides of (7.5) and (7.7). These relations can be exploited to obtain mapping properties of the parametrix based surface and volume potentials taking into account those mapping properties already known for the analogous surface and volume potentials constructed with the fundamental solution of the Laplace equation.
One of the main differences with respect to the bounded domain case is that the integrands of the operators V , W , P and R and their corresponding direct values and conormal derivatives do not always belong to L 1 . In these cases, the integrals should be understood as the corresponding duality forms (or their their limits of these forms for the infinitely smooth functions, existing due to the density in corresponding Sobolev spaces).  Proof. Let us prove first the mapping property (8.1). From Theorem 7.1, we have that V g ∈ H 1 (Ω) for some g ∈ H −1/2 (S). Hence, it suffices to prove that V g ∈ L 2 (ω; Ω).
Differentiating using the product rule, we can write Taking into account relation (7.5) and applying (8.3) to h = V ∆ (g/a), we get By virtue of the mapping property for the operator V provided by Theorem 7.1, the last term belongs to L 2 (ω; Ω) due to the fact that V ∆ (g/a) = V g ∈ H 1 (Ω), and thus its derivatives belong to L 2 (ω; Ω). The term ∇a acts as a multiplier in the space L 2 (ω; Ω) due to Condtion 6.1. On the other hand, the term a∆V ∆ (g/a) vanishes on Ω since V ∆ (·) is the single layer potential for the Laplace equation, i.e., V ∆ (g/a) is a harmonic function. This, completes the proof for the operator V . The proof for the operator W follows from a similar argument.
Condition 8.1. In addition to Condition 3.1 and Condition 6.1, we will sometimes need the following condition: Remark 9. Note as well that due to Condition 3.1 and the continuity of the function ln a, the components of ∇(ln a) and ∆(ln a) are bounded as well.
Theorem 9.1. The following operators are continuous under Condition 6.1, Proof. Let g ∈ H −1 (R 3 ). Then, by virtue of the relation (7.3) Pg = P ∆ (g/a). Since Condition 6.1 holds, (g/a) ∈ H −1 (R 3 ) and thefore the continuity of the operator P follows from the continuity of P ∆ : , which at the same time implies the continuity of the operator (9.3), see [9, Theorem 4.1] and more references therein. Let us prove now the continuity of the operator R. Due to the second condition in (6.1), the components of ∇a ∈ L 2 (R 3 ) behave as multipliers in the space L 2 (ω −1 ; R 3 ). Let g ∈ L 2 (ω −1 ; R 3 ), then the relation (7.4) applies and gives Here, ∇ ln a is multipliers under Condition 6.1 in the space H −1 (R 3 ).
As Rg ∈ H 1 (Ω), we only need to prove that ∆Rg(y) ∈ L 2 (ω; Ω). Using the relation Taking the trace and the conormal derivative of (10.2), we obtain integral representation formulae for the trace and traction of u respectively: For some distributions f, Ψ and Ψ, we consider a more indirect integral relation associated with the third Green identity (10.2) u + Ru − V Ψ + W Φ = Pf, in Ω.
Proof. To prove that u ∈ H 1,0 (Ω; A), taking into account that by hypothesis u ∈ H 1 (Ω), so there is only left to prove that Au ∈ L 2 (ω; Ω). Firstly we write the operator A as follows: It is easy to see that the second term belongs to L 2 (ω; Ω). Keeping in mind Remark 7 and the fact that u ∈ H 1 (Ω), then we can conclude that the term u∇a ∈ H 1 (Ω) since due to the second condition in (6.1) ∇a is a multiplier in the space H 1 (Ω) and therefore ∇(u∇a) ∈ L 2 (ω; Ω). Now, we only need to prove that ∆(au) ∈ L 2 (ω; Ω). To prove this we look at the relation (10.5) and we put u as the subject of the formula. Then, we use the potential relations (7.3), (7.5) and (7.7) In virtue of the Theorem 9.2, Ru ∈ L 2 (ω; Ω). Moreover, the terms in previous expression depending on V ∆ or W ∆ are harmonic functions and P ∆ is the newtonian potential for the Laplacian, i.e. ∆P ∆ ( f a ) = f a . Consequently, applying the Laplacian operator in both sides of (10.7), we obtain: Thus, ∆u ∈ L 2 (ω; Ω) from where it immediately follows that ∆(au) ∈ L 2 (ω; Ω). Hence u ∈ H 1,0 (Ω; A). The rest of the proof is equivalent to [25,Lemma 5.1].
The proof of the following statement is the counterpart of [25, Lemma 5.2] for exterior domains. The proof follows from the invertibility of the operator V ∆ , see [18,Corollary 8.13]. (10.9) then Ψ * (y) = 0.
Proof. Take the trace of (10.9) and relation (7.5), to obtain VΨ * (y) = V ∆ Ψ * a (y) = 0, y ∈ S. 11. BDIES. In this section, we will derive a system of boundary domain integral equations formally segregated from the solution u of the BVP (4.1)-(4.3), following a similar approach as in [7,Section 5]. Consequently, we introduce Φ 0 ∈ H 1/2 (S) and Ψ 0 ∈ H −1/2 (S) as continuous fixed extensions to S of the functions φ 0 ∈ H 1/2 (S D ) and ψ 0 ∈ H −1/2 (S N ). Moreover, let φ ∈ H 1/2 (S N ) and ψ ∈ H −1/2 (S D ) be arbitrary functions formally segregated from u. Then, make in the three third Green identities (10.2)-(10.4) to obtain the following BDIES (M12) In what follows, we will denote by X the vector of unknown functions . We will denote by M 12 the matrix operator that defines the system (M 12): 4) and by F 12 the right hand side of the system F 12 = [ F 0 , γ + F 0 − Φ 0 ] . Using this notation, the system (M12) can be rewritten in terms of matrix notation as M 12 X = F 12 .
If Condition 6.1 and Condition 8.5 hold, then, due to the mapping properties of the potentials, F 12 ∈ F 12 ⊂ Y 12 , while operators M 12 : H → F 12 and M 12 : X → Y 12 are continuous. Here, we denote The following result shows that the BDIES (M12) is equivalent to the original mixed BVP (4.1)-(4.3).

Representation theorems and invertibility.
In this section, we aim to prove the invertibility of the operator M 12 : H → F 12 by showing first that the arbitrary right hand side F 12 from the respective spaces can be represented in terms of the parametrix-based potentials and using then the equivalence theorems.
The following result is the counterpart of [9, Lemma 7.1] for the new parametrix P x (x, y). The analogous result for bounded domains can be found in [7, Lemma 3.5]. Proof. Let us assume that such functions f * and Ψ * , satisfying (12.1), exist. Then, we aim to find expressions of these functions in terms of F * . Applying the potential relations (7.5), (7.3) to the equation (12.1), we obtain Applying the Laplace operator at both sides of the equation (12.2), we get On the other hand, we can rewrite equation (12.2) as V ∆ Ψ * a (y) = Q(y), y ∈ Ω, (12.4) where Q(y) := F * (y) − P ∆ (∆F * ) (y). (12.5) Now, we take the trace of (12.4) It is well known that the direct value operator of the single layer potential for the Laplace equation V ∆ : H −1/2 (S) −→ H 1/2 (S) is invertible (cf. e.g. [18,Corollary 8.13]). Hence, we obtain the following expresion for Ψ * : Relations (12.3) and (12.7) imply the uniqueness of the couple (f * , Ψ * ). Now, we just simply need to prove that the pair (f * , Ψ * ) given by (12.7) and (12.3) satisfies (12.1). For this purpose, let us note that the single layer potential operator, V ∆ (Ψ * /a) with Ψ * given by (12.7), as well as Q(y) given by (12.5) are both harmonic functions. Since Q(y) and V ∆ (Ψ * /a) are two harmonic functions that coincide on the boundary due to (12.6), then they must be identical in the whole Ω due to the uniqueness of solution to the Dirichlet problem for the Laplace equation, see [ Then there exists a unique triplet (f * , is a linear an bounded operator and (F 0 , F 1 ) are given by Proof. Taking Φ * = γ + F 0 −F 1 and applying the previous lemma to F * = F 0 +W Φ * we prove existence of the representation (13.1) and (13.2). The uniqueness follows from the homogenenous case when F 0 = F 1 = 0. Then, (13.2) implies Φ * = 0 and consequently, by (13.1) and Lemma 12.1, we get Ψ * = 0 and f * = 0.
We are ready to prove one of the main results for the invertibility of the matrix operator of the BDIES (M12).
Theorem 13.1. If Conditions (6.1) and (8.5) hold, then the following operator is continuous and continuously invertible: Proof. In order to prove the invertibility of the operator M 12 : H −→ F 12 , we apply the Corollary 13 to any right-hand side F 12 ∈ F 12 of the equation M 12 U = F 12 . Thus, F 12 can be uniquely represented as (f * , Ψ * , Φ * ) = C * F 12 as in (13.1)-(13.2) where C * : F 12 −→ L 2 (ω; Ω) × H −1/2 (S) × H 1/2 (S) is continuous. In virtue of the equivalence theorem for the system (M12), Theorem 11.1, and the invertibility theorem for the boundary value problem with mixed boundary conditions, Theorem 4.1, the matrix equation M 12 U = F 12 has a solution U = (M 12 ) −1 F 12 where the operator (M 12 ) −1 , is given by expressions where (f * , Ψ * , Φ * ) = C * F 12 . Consequently, the operator (M 12 ) −1 is a continuous right inverse to the operator (13.3). Moreover, the operator (M 12 ) −1 results to be a double sided inverse in virtue of the injectivity implied by Theorem 11.1.
14. Fredholm properties and invertibility. In this section, we are going to benefit from the compactness properties of the operator R to prove invertibility of the operator M 12 : X → Y 12 . This invertibility result is more general than the one presented in the previous section. The price to pay is imposing an additional condition on the variable coefficient. Unlike as in the bounded case, see similar to [9, Section 7.2], the Rellich compact embedding theorem, see e.g. [18,Theorem 3.27], cannot be applied as Ω is a bounded domain. Still, we can overcome this obstacle by decomposing the operator R into the sum of two operators: one which can be made arbitrarily small and the other one will be compact. Then, we shall simply make use of the Fredholm alternative to prove the invertibility of the matrix operator that defines the Proof. Let B(0, r) be the ball centered at 0 with radius r big enough such that S ⊂ B r . Furthermore, let χ ∈ D(R 3 ) be a cut-off function such that χ = 1 in S ⊂ B r , χ = 0 in R 3 B 2r and 0 ≤ χ(x) ≤ 1 in R 3 . Let us define by R c g := R(χg), R s g := R((1 − χ)g).
We will prove first that the norm of R s can be made infinitely small. Let g ∈ H 1 (Ω), then Consequently, we have the following estimate: Using the previous estimate is easy to see that when r →+∞ the norm R s g H 1 (Ω) tends to 0 due to Conditions (14.1). Hence, the norm of the operator R s can be made arbitrarily small.
To prove the compactness of the operator R c g := R(χg), we recall that supp(χ) ⊂ B(0, 2r). Then, one can express R c g := R Ωr ([χg| Ωr ]) where the operator R is defined now over Ω r := Ω ∩ B 2r which is a bounded domain. As the restriction operator | Ωr : H 1 (Ω) −→ H 1 (Ω r ) is continuous, in virtue of Theorem 9.1, the operator R c g : L 2 (Ω r ) −→ H 1 (Ω r ) is also continuous. Due to the boundedness of Ω r , we have H 1 (Ω r ) = H 1 (Ω r ) and thus the compactness of R c g follows from the Rellich Theorem (see [18,Theorem 3.27]) applied to the embedding L 2 (Ω r ) ⊂ H 1 (Ω r ). is Fredholm with zero index.
Proof. Using the previous Lemma, we have R = R s + R c so R s < 1 hence I + R s is invertible. On the other hand R c is compact and hence I + R s a compact perturbation of the operator I + R, from where it follows the result. Then, U will also solve the following extended system Furthermore,every solution of the system (15.2) will solve the equationM 12 0 U = F. The system (15.2) can be written also in matrix form as M 12 0 U = F where F denotes the right hand side and M 12 0 is defined as We note that the three diagonal operators: is a compact perturbation of M 12 . Consequently, M 12 is Fredholm with index zero. In addition, as the operator M 12 is one to one, we conclude that it is also continuously invertible. 16. Conclusions. A new parametrix for the diffusion equation in non homogeneous media (with variable coefficient) has been analysed in this paper. Mapping properties of the corresponding parametrix based surface and volume potentials have been shown in corresponding weigthed Sobolev spaces depending on several regularity and decay conditions on the variable coefficient a(x). A BDIES for the original BVP has been obtained. Results of equivalence between the BDIES and the BVP have been shown along with the invertibility of the matrix operator defining the BDIES using Fredholm alternative arguments overcoming the technicalities that unbounded domains present. Now, we have obtained an analogous system to the BDIES (M12) of [9] with a new family of parametrices which is uniquely solvable. Hence, further investigation about the numerical advantages of using one family of parametrices over another will follow.
Further generalised can be obtained by relaxing the smoothness of the boundary to Lipschitz domains. In this case, one needs the generalised canonical conormal derivative operator defined in [19,20]. Another possible generalisation could consider relaxing the smoothness of the coefficient, see [20].