A New Family of Boundary-Domain Integral Equations for the Dirichlet Problem of the Diffusion Equation in Inhomogeneous Media with $H^{-1}(\Omega)$ Source Term on Lipschitz Domains

The interior Dirichlet boundary value problem for the diffusion equation in non-homogeneous media is reduced to a system of Boundary-Domain Integral Equations (BDIEs) employing the parametrix obtained in (Fresneda-Portillo, 2019) different from (Chkadua et. al 2009). We further extend the results obtained in (Fresneda-Portillo, 2019) for the mixed problem in a smooth domain with $L^{2}(\Omega)$ right hand side to Lipschitz domains and source term $f$ in the Sobolev space $H^{-1}(\Omega)$, where neither the classical nor the canonical co-normal derivatives are well defined. Equivalence between the system of BDIEs and the original BVP is proved along with their solvability and solution uniqueness in appropriate Sobolev spaces.


Introduction
The popularity of the boundary integral equation method (BIE) is owed to the reduction of dimension of a boundary value problem (BVP) with constant coefficients and homogeneous right hand side defined on a domain of R n . By applying the BIE method, one can reformulate the original BVP in terms of an equivalent integral equation defined exclusively on the boundary of the domain. This method has already been extensively studied for many boundary value problems, for instance: Laplace, Helmholtz, Stokes, Lamé, etc. [3,4,5]. This method requires an explicit formula for the fundamental solution of the PDE operator in the BVP which is not always available when the BVP has variable coefficients [2,6].
The overcome this issue, one can construct a parametrix (Levi function) [2,Section 3] for the PDE operator and use it to derive an equivalent system of Boundary-Domain Integral Equations following a similar approach as for the BIE method. However, the reduction of dimension no longer applies as volume integrals will appear in the new formulation as a result of the remainder term. This is also the case for non-homogeneous problems with constant coefficients, [3,Chapter 1 and 2].
In order to preserve the reduction of dimension, one can use the method of radial integration method (RIM) which allows to transform volume integrals into boundary only integrals [14]. This method has been successfully implemented to solve boundary-domain integral equations derived from BVPs with variable coefficients [15,13]. This method is also able to remove various singularities appearing in the domain integrals.
The recent developments on numerical approximation of the solution of BDIEs show that there are effective and fast algorithms able to compute the solution. For example: the collocation method [7,9] which, although leads to fully populated matrices, it can be further enhanced by using hierarchical matrix compression and adaptive methods as shown in [10] to reduce the computational cost. Localised approaches to reduce the matrix dimension and storage have also been developed [11,12] which lead to sparse matrices.
Moreover, reformulating the original BVP in the Boundary Domain Integral Equation form can be beneficial, for instance, in inverse problems with variable coefficients [17].
On the one hand, the family of weakly singular parametrices of the form P y (x, y) for the particular operator has been studied extensively studied [2,18,19,20]. Note that the superscript in P y (x, y) means that P y (x, y) is a function of the variable coefficient depending on y, this is P y (x, y) = P (x, y; a(y)) = −1 4πa(y)|x − y| .
On the other hand, is another parametrix for the same operator A. In this case, the parametrix depends on the variable coefficient a(x). This parametrix was introduced in [1] for the mixed problem in smooth 3D domains and in [21] for the mixed problem in Lipschitz domains. Some preliminary results for the mixed problem in exterior domains have also been obtained [22].
However, most of the numerical methods to solve BDIEs aforementioned are tested for the Dirichlet problem [10,7,9,15]. In order to compare the performance of parametrices P x (x, y) and P y (x, y), one needs first to prove the equivalence between the original Dirichlet BVP and the system of BDIEs as well as the uniqueness of solution (well-posedness) of the system of BDIEs what corresponds to the main purpose of this paper.
The study of new families of parametrices is helpful at the time of constructing parametrices for systems of PDEs as shown in [1, Section 1] for the Stokes system. In this case, the fundamental solution for the pressure does not present any relationship with the viscosity coefficient whereas the parametrix for the pressure depends on two variable viscosity coefficients: one depending on y and another depending on x, see also [27].
The parametrix preserves a strong relation with the fundamental solution of the corresponding PDE with constant coefficients. Using this relation, it is possible to establish further relations between the surface and volume potential type operators for the variable-coefficient case with their counterparts for the constant coefficient case, see, e.g. Different families of parametrices lead to different relations with their counterparts for the constant coefficient case. For the parametrices considered in this paper, these relations are rather simple, which makes it possible to obtain the mapping properties of the integral potentials in Sobolev spaces and prove the equivalence between the BDIE system and the BVP. After studying the Fredholm properties of the matrix operator which defines the system, its invertibility is proved, what implies the uniqueness of solution of the BDIE system.
In this paper, we extend the results obtained in [22,1,26] by considering the source term of the equation Au = f in the Sobolev space H −1 (Ω). This happens for example, when the source term f is the Dirac's delta distribution. The Dirac's delta is an example of distribution that does not belong to the space L 2 but belongs to H −1 and is used in many applications in Physics, Engineering and other mathematical problems [23,24,25].
This generalisation for the source term introduces an additional issue on the definition of the co-normal derivative which is needed to derive BDIEs.
The co-normal derivative operator is usually defined with the help of first Green identity, since the function derivatives do not generally exist on the boundary in the trace sense. However this definition is related to an extension of the PDE and its right-hand side from the domain Ω, where they are prescribed, to the boundary of the domain, where they are not. Since the extension is non-unique, the co-normal derivative appears to be non-unique operator, which is also non-linear in u unless a linear relation between u and the PDE right-hand side extension is enforced. This creates some difficulties particularly in formulating the BDIEs.
To overcome these technical issues, we introduce a subspace of H 1 (Ω) which is mapped by the PDE operator into the space H − 1 2 (Ω) for the right hand side [16]. This allows to define an internal co-normal derivative operator, which is unique, linear in u and coincides with the co-normal derivative in the trace sense if the latter does exist. This approach is applied to the formulation and analysis of direct segregated BDIEs equivalent to the stated Dirichlet BVP with a varaible coefficent and right hand side from H −1 (Ω).
Last but not the least, we generalise in this paper the results for the two-dimensional case and smooth boundary domains [26]. We shall consider the partial differential equation where the variable smooth coefficient a(x) ∈ C 2 (Ω) is such that 0 < a min a(x) a max < ∞, ∀x ∈ Ω, a min , a max ∈ R, u(x) is the unknown function and f is a given distribution on Ω. It is easy to see that if a ≡ 1 then, the operator A becomes ∆, the Laplace operator.
In what follows D(Ω) := C ∞ comp (Ω) denotes the space of Schwartz test functions, and D * (Ω) denotes the space of Schwartz distributions, H s (Ω), H s (∂Ω) denote the Bessel potential spaces, where s ∈ R (see e.g. [4,3] for more details). We recall that the spaces For u ∈ H 1 (Ω), the partial differential operator A is understood in the sense of distributions. Using the usual notation of distribution theory, equation (2.1) can be written as Using the differentiation properties of distributions, one can obtain the following identity Au, v Ω = − a∇u, ∇v Ω , ∀v ∈ D(Ω).
To simplify the notation, we introduce the operator E defined as follows which allow us to write equation (2.2) as 3) defines the continuous linear Let us also define the so-called aggregate operator[28, Section 3.1] of A, asǍ : where the bilinear functional E : For any u ∈ H 1 (Ω), the functionalǍu belongs to H −1 (Ω) and is an extension of the functional Au ∈ H −1 (Ω) whose domain is thus extended from H 1 (Ω) to the domain 3 Traces, co-normal derivatives and Green identities From the trace theorem for Lipschitz domains, we know that the trace of a scalar function For u ∈ H s (Ω), s > 3/2, we can define on ∂Ω the conormal derivative operators, T ± , in the classical sense where n + (x) is the exterior unit normal vector directed outwards the interior domain Ω at a point x ∈ ∂Ω. Respectively, n − (x) is the unit normal vector directed inwards the interior domain Ω at a point x ∈ ∂Ω. Sometimes, we will also use the notation T ± x u or T ± y u to emphasise the differentiation variable. When the variable of differentiation is obvious or is a dummy variable, we will simply use the notation T ± u.
It is well known that the classical co-normal derivative operator is generally not well [28,29]. Consequently, to correctly define a conormal derivative, one can draw on the first Green identity. This is the case for the generalised co-normal derivative and the canonical co-normal derivatives [28, Definition 3.1 and 3.6].
If u, v ∈ H 1 (Ω), u satisifying Au = r Ωf in Ω for somef ∈ H −1 (Ω), then, the first Green identity holds in the following form In order to appropriately define the canonical co-normal derivative[30, Definition 6.5], we introduce the following space.
Definition 2. Let s ∈ R and A * : H s (Ω) −→ D * (Ω) be a linear operator. For t ≥ − 1 2 , we introduce the space In this paper, A * will refer to either A or ∆ in the above definition. Also, we remark, whereÃu :=E(Au).
The canonical co-normal derivatives T + u is independent of (non-unique) choice of the operator γ −1 ; independent of the source termf , unlike to generalised co-normal derivative defined in (3.1); it is linear with respect to u and has the following continuous mapping and v ∈ H 1 (Ω), then, the first Green identity for the canonical co-normal derivative holds in the following form, (cf. [20,Theorem 2.13]).
Furthermore, if u ∈ H 1,− 1 2 (Ω; A) and v ∈ H 1 (Ω), then, the first Green identity for the canonical co-normal derivative holds in the following form[20, Theorem 2.13] In the particular case of u ∈ H 1,0 (Ω; A) and v ∈ H 1 (Ω), then, the first Green identity takes the form To obtain the second Green identity for u ∈ H 1,− 1 2 (Ω; A) and v ∈ H 1 (Ω), we use the first Green identity for the canonical co-normal derivative for u, i.e. identity (3.2) and subtract it from the first Green identity for the generalised co-normal derivative for v, this is, swapping u by v in formula (3.1). Hence, supposing that r Ω Av = r Ωf with f ∈ H −1 (Ω), we obtain the following second Green identity If u, v ∈ H 1,− 1 2 (Ω; A), then we arrive at the familiar form of the second Green identity for the canonical extension and canonical co-normal derivatives In the particular case, when u, v ∈ H 1,0 (Ω; A), the previous identity becomes (3.5)
To obtain a system of boundary-domain integral equations for the boundary value problem (4.1a)-(4.1b), we intend to use Boundary Integral Method (BIM) approach [5].
However, this method requires an explicit fundamental solution which is not always available when the PDE differential operator has variable coefficients, as it is the case for the operator A. To overcome this problem, one can introduce a parametrix [2,1,6].
Definition 4. A distribution P (x, y) in two variables x, y ∈ R 3 is said to be a parametrix or Levi function for a differential operator A x differentiating with respect to x, if the following identity is satisfied where δ(.) is the Dirac distribution and R(x, y) is remainder.
A parametrix for a given operator A x might not be unique. This is the case, for example, for the operator A. One parametrix [11,2] is given by where P ∆ (x − y) = −1 4π|x − y| is the fundamental solution of the Laplace equation. The remainder corresponding to the parametrix P y is given by In this paper, for the same operator A defined in (2.1), we will use another parametrix [1] P (x, y) : which leads to the corresponding remainder Note that the both remainders R x and R y are weakly singular, i.e., This is due to the smoothness of the variable coefficient a(·).

Volume potentials
Following the steps of the boundary-domain integral method [2], we will later on replace u by the parametrix (4.3) in the first Green identity (3.2). This will give an integral representation formula of the solution u in terms of surface and volume potential-type integral operators. In this section, we define these surface and volume integral operators and study their mapping properties which will be applied to prove the main results of this paper.
For the function ρ defined on Ω ⊂ R 3 , e.g., ρ ∈ D(Ω) the volume potential and the remainder potential operator, corresponding to parametrix (4.  Since Ω is a bounded domain, then, the compact embedding theorem for Sobolev spaces[4, Chapter 2] can be applied to the remainder operators (5.11)-(5.13) of the previous theorem to obtain the following corollary.

Surface potentials
The single layer potential operator and the double layer potential operator associated with the Laplace equation ∆u = 0, are defined as where T + ∆ is the normal derivative operator (i.e. T + with a(x) ≡ 1 in (3)). Similarly, one can defined the corresponding potentials parametrix-based for y ∈ R 3 and y / ∈ ∂Ω, as Due to (4.3) and the fact that the operators V and W can be also expressed in terms the surface potentials and operators associated with the Laplace operator, Since the mapping properties in Sobolev spaces of the single layer potential and double layer potential for the Laplace equation are well known [20,29], one can easily derive analogous mapping properties for the operators V and W as a consequence of the relations (5.16) and (5.17), along with Theorems 3.3-3.7 proved in [20]. These mapping properties are reflected in the following results which are included for completeness of the paper as some are key to prove the main results.
Let Ω be a bounded Lipschitz domain, let 1 2 < s < 3 2 . Then, the following operators are bounded, µr The following result follows from the relations (5.16)-(5.17) and the analogous jump properties for the harmonic surface potentials [20, Theorem 3.6].
6 Integral representation of the solution in terms of the surface and volume potentials In this section, we will obtain an integral representation formula for the solution u of the original BVP (4.1a)-(4.1b). These results will be useful to construct a system of BDIEs equivalent to the original Dirichlet BVP. We will follow a similar approach as in [20] but using the new parametrix (4.3).
(iii) The trace of u,i.e. γ + u, can be represented in terms of the surface and volume potentials as follows Proof. Let us prove item (i). First, we consider the first Green identity (3.2) with the roles of u and v interchanged In order to apply the first Green identity, we needed u ∈ H 1 (Ω) and v ∈ H 1,0 (Ω; A). Let us take v := P (x, y) as the parametrix. Let us remark that as P is the parametrix, then For every distribution u ∈ H 1 (Ω), AP, u Ω ∈ H 1 (Ω), and hence in L 2 (Ω) due to the mapping properties of the operator R given in Theorem 5.1. This implies that, as a distribution, P ∈ H 1,0 (Ω; A) Then, the identity (6.3) can now be reformulated in terms of the surface and volume integral operators as Since E(u, P ) = − AP, u Ω and AP = δ + R for being P a parametrix, we obtain what implies the result.
Let us prove now item [(ii)]. Since r Ω Au =f withf ∈ H −1 (Ω), then, we need to use the generalised second Green identity (3.4), again swapping the roles of u and v As before, make v = P (x, y) to obtain a new representation formula in terms of the parametrix-based surface and volume potentials Taking into account that AP = δ + R and rearranging terms, we obtain Item (iii) directly follows from item (ii) by taking the trace of (6.9), keeping in mind the jump property γ + W γ + u = − 1 2 γ + u + Wγ + u given by Corollary 2 and the mapping properties given in Corollary 3.
To derive the boundary-domain integral equation systems, we will use the integral representation formulas obtained in the previous theorem. However, we will substitute that both the trace and generalised conormal derivatives are independent from u. Hence, we will use the distributions Ψ and Φ in their place as unknowns alongside u, and consider the new boundary domain integral equation (6.10) We will show now that any triple (u, Ψ, Φ) satisfying the previous relation, solves the PDE (4.1a).
The following two statements are a generalisation of Lemma 9 and Lemma 10 in [21] to the case wheref ∈ H −1 (Ω).
The following Lemma is a direct consequence of the invertibility of the direct value of the single layer potential for the Laplace equation [4,Corollary 8.13]. A proof of the Lemma is available in [21].

BDIE system for the Dirichlet problem
We aim to obtain a segregated boundary-domain Let us represent the generalized co-normal derivative and the trace of the function u as T + (f , u) = ψ, γ + u = ϕ 0 , and we will regard the new unknown function ψ ∈ H − 1 2 (∂Ω) as formally segregated of u. Thus we will look for the couple (u, ψ) ∈ H 1 (Ω) × H − 1 2 (∂Ω). To obtain one of the possible boundary-domain integral equation systems we will use equation (6.10) in the domain Ω, and equation (6.2) on ∂Ω. Then we obtain the following system (A1) of two equations for two unknown functions, in Ω, (7.1a) Note that for ϕ 0 ∈ H 1 2 (∂Ω), we have the inclusion F 0 ∈ H 1 (Ω) iff ∈ H −1 (Ω) due to the mapping properties of the surface and volume potentials given in Theorem 5.1, Theorem 5.2 and Corollary 3.
The system (A1), given by (7.1a)-(7.1b) can be written in matrix notation as where U represents the vector containing the unknowns of the system, the right hand side vector is and the matrix operator A 1 is defined by: We note that the mapping properties of the operators involved in the matrix imply the continuity of the operator Let us prove that the Dirichlet boundary value problem (4.1) in Ω is equivalent to the system of the Boundary Domain Integral Equations (7.1a)-(7.1b). ii) If a couple (u, ψ) ⊤ ∈ H 1 (Ω) × H − 1 2 (∂Ω) solves the BDIE system (A1) then u solves the BVP and the functions ψ satisfy (7.3).
Let us now prove the invertibility of the operator A 1 Theorem 7.2. The operator is invertible.
Proof. To prove the invertibility, let A 1 0 be the matrix operator defined by  Hence, further investigation about the numerical advantages of using one family of parametrices over another will follow. Now, numerical methods can be applied when the source term belongs to H −1 (Ω) [25,23,24].
Further work will consist of extending the results presented in this paper to unbounded domains, non-smooth coefficients or other BVP problems with different boundary conditions as well as providing a localised version of the BDIEs (A1) presented in this paper, inspired by the works, [2,11].
We highlight again that analysing BDIEs for different parametrices, i.e. depending on the variable coefficient a(x) or a(y), is crucial to understand the analysis of BDIEs derived with parametrices that depend on the variable coefficient a(x) and a(y) at the same time, as it is the case for the Stokes system [1,27].