A discontinuous Poisson–Boltzmann equation with interfacial jump: homogenisation and residual error estimate

A nonlinear Poisson–Boltzmann equation with inhomogeneous Robin type boundary conditions at the interface between two materials is investigated. The model describes the electrostatic potential generated by a vector of ion concentrations in a periodic multiphase medium with dilute solid particles. The key issue stems from interfacial jumps, which necessitate discontinuous solutions to the problem. Based on variational techniques, we derive the homogenisation of the discontinuous problem and establish a rigorous residual error estimate up to the first-order correction.


Introduction
In this paper, we consider the steady-state problem of a nonlinear Poisson-Nernst-Planck (PNP) system, which describes multiple concentrations of charged particles (e.g. ions) subject to a selfconsistent electrostatic potential calculated from Poisson's equation. In particular, we shall investigate the PNP model on a multiphase medium. The prototypical multiphase medium in mind consists of an electrolyte medium (pores), which surrounds disjoint solid particles. Such models have numerous applications describing electro-kinetic phenomena in bio-molecular or electro-chemical models, photo-voltaic systems and semiconductors, see e.g. [1][2][3][4][5][6] and references therein. Our specific interests are motivated by models of Li-Ion batteries, see the relevant references [7][8][9][10].
We shall deal with the nonlinearity of the model within an analytic framework, in which the PNP system can be transformed into an equivalent scalar Poisson-Boltzmann (PB) equation, see e.g. [11][12][13]. This is possible, when reaction terms in the charged particle fluxes are omitted and the equations for the concentrations decouple since the charged particle concentrations are explicitly determined by the corresponding Boltzmann statistics. At the interface boundaries, this implies homogeneous Neumann conditions, which nevertheless allow for jumps in the concentrations of the species. For references applying linearisation of the PNP equations near the Boltzmann distribution see [14,15].
CONTACT Victor A. Kovtunenko victor.kovtunenko@uni-graz.at In order to be able to homogenise the equations over the entire pores-and-particles medium, a physically consistent continuation of the governing relations in the particles is needed. Here, we assume the Gouy-Chapman-Stern model for electric double layers (EDLs). [16] By this, charged species underlie transport through diffusion and electrostatic drift in the pores and pure diffusion in the solid phase. This model proposes a jump of the electrostatic field across the interface (a voltage drop) and a current prescribed at the interior boundary of the solid particles.
In our general setting, we assume that the current may have a jump across the interface as well. Allowing for jumps extends the modelling to describe e.g. a separator in super-capacitors.
In the following, we will derive a discontinuous formulation of the PB equation (valid both on the volume occupied by the solid particles and on the surrounding porous space). The key feature addressed in this manuscript is the imposed inhomogeneous jump conditions, which complement the PB equation (see (8)

below).
We emphasise that, from the point of view of partial differential equations, the considered PB equations are characterised by a strongly nonlinear term, which is unbounded and features exponential growth, see e.g. [17].
A first aim of this paper is to establish a proper variational setting of the problem, while a second part deals with its rigorous homogenisation. With respect to the later, the averaged effective coefficients of the limit problem represent the macroscopic behaviour of the EDL, which is of primary practical importance.
For reference concerning the classic homogenisation theories, we refer to [18][19][20][21][22][23]. The applied methods range from two-scale convergence (see e.g. [24]) over Gamma-convergence (see e.g. [25]) to unfolding (see [26]) and others. While formal methods of averaging are widely used in the literature, their verification in terms of residual error estimates is a hard task.
From the point of view of homogenisation, the principal difficulty of discontinuous problems concerns the nonstandard boundary conditions with jumps: firstly, related jump conditions are inherent for cracks. For models and methods used in crack problems, we refer to [27][28][29] and references therein. From a geometric viewpoint, cracks are open manifolds in the reference domain. Hence, classic Poincare-Friedrichs-Korn inequalities are valid in such situations. In contrast to cracks, the interfaces here are assumed to be closed manifolds disconnecting the reference domain. This fact requires discontinuous versions of Poincare-Friedrichs-Korn inequalities, which are then applied for semi-norm estimates.
Secondly, the interface boundary conditions are of Robin type. The homogenisation results known for linear problems with Robin (also called Fourier) conditions are crucially sensitive to the asymptotic rates of the involved homogenisation parameters. This issue concerns the coefficients in the boundary condition (cf. Lemma 1 below) and the volume fraction of solid particles in periodic cells (cf. Lemma 2 below), see e.g. [24,[30][31][32].
In the literature, on the one hand, the fully coupled PNP as well as decoupled PB equations are treated mostly under homogeneous boundary conditions of the Neumann type, see [11,12,14]. On the other hand, even inhomogeneous Neumann and also Robin boundary conditions are homogenised for the most part of linear equations, see e.g. [26,31]. In e.g. [4,6] only the pore model is described and extended in the solid phase inexplicitly.
Homogenisation of transmission problems with interface jumps can be found in works concerning models of diffusion in a partially fractured porous medium, models of heat conduction in composite materials or in models of absorption of a dissolved chemical in a fluid flowing through a porous medium. We refer to [33][34][35][36] for linear transmission problems as well as for linear problems under nonlinear [37] and contact [38] boundary conditions, and to [39] for a semi-linear case with bounded nonlinearity.
In comparison to the works mentioned above, the principal challenge of the underlying, discontinuous PB equation is the strongly nonlinear term together with nonstandard jump conditions of Robin type, which allow not only discontinuous fields but also discontinuous fluxes across a separator. In the present contribution, we homogenise it and derive the averaged limit problem. A further major result is a rigorous estimate of the residual error up to the first-order correction.
For these purposes, we develop a variational technique based on orthogonal Helmholtz decomposition following the lines of [21,23]. In a periodic cell, we decompose oscillating coefficients (describing the electric permittivity) by using the nontrivial kernel in the space of vector valued periodic functions, which is represented by sums of constant and divergence free (and thus, skew symmetric) vector fields (cf. Lemma 3). Employing solutions of appropriately defined discontinuous cell problems, we obtain a regular decomposition of the homogenisation problem (see Theorem 2).
A second result establishes the critical rates of the asymptotic behaviour with respect to a homogenisation parameter ε 0 + for coefficients in the inhomogeneous transmission condition: we find on the one side that the critical rate for the coefficient by interfacial jumps is 1 ε . This factor occurs in the discontinuous Poincare inequality (for the norm squared, cf. (21) below) and is thus relevant for a coercivity estimate, which in return contributes to the solvability of the discontinuous problem and the subsequent estimate of the homogenisation error.
On the other side, the critical rate for the flux prescribed at the interior boundary of solid particles is ε. At this rate, the interior boundary flux induces an additional potential, which distributes over the macroscopic domain in the homogenisation limit ε 0 + . If the asymptotic rate is lower than the critical one, then this flux vanishes in the limit. Otherwise, if the asymptotic rate is bigger, then the flux term diverges.
From the above description, we summarise the key points of this paper as follows: • the study of inhomogeneous interfacial conditions describing EDL; • the combination of the strongly nonlinear term of exponential type, jumps and Robin conditions; • a variational framework of the discontinuous problem; • the performing of the homogenisation procedure with rigorous error estimates; and • the identification of the critical asymptotic rates of the boundary coefficients.
Outline: In Sections 2.1-2.3, we first present the problem geometry, the physical and the mathematical model. Section 2.3 establishes moreover the equivalence of the steady state of the PNP system with the scalar Poisson-Boltzmann (PB) equation and the existence of a unique solution to the PB equation (see Theorem 1).
In Section 3, we consider the homogenisation problem and the residual error estimate. At first, we state three auxiliary Lemmata before stating the main homogenisation Theorem 2.
Finally, Section 4 provides a brief discussion of the obtained results.

Statment of the Problem
We start with the description of the geometry.

Geometry
Let ω denote the domain (which is an open set) occupied by solid particles of general shape (either single or multiple particles), which are located inside the unit cell ϒ = (0, 1) d ⊂ R d , d = 1, 2, 3. We assume that all particles ω ⊂ ϒ are disjunctively located as well as bounded away from the boundary ∂ϒ, i.e. ω ∩ ∂ϒ = ∅. We assume that the boundary ∂ω is Lipschitz continuous with outer normal vector ν = (ν 1 , . . . , ν d ) pointing away from the domain ω. Moreover, we distinguish the positive (outward orientated) surface ∂ω + and the negative (inward orientated) surface ∂ω − as the faces of the boundary ∂ω, when approaching the boundary ∂ω from outside, i.e. from ϒ \ ω or from the inside, i.e. from ω, respectively. For a two-dimensional example configuration see the illustration in Figure 1(a).
In the following, we consider a fixed, small homogenisation parameter ε ∈ R + and pave R d with periodic cells ϒ ε p indexed by p ∈ N. The periodic cells ϒ ε p are constructed from ϒ in the following ω ϒ/ω ∂ω − ∂ω + ν ϒ/∂ω (a) unit cell (b) periodic domain Ω/∂ω ε way: The position of every spatial point x = (x 1 , . . . , x d ) ∈ R d can be decomposed as into the integer-valued floor function coordinates x ε ∈ Z d and the fractional coordinates { x ε } ∈ ϒ. We shall then enumerate all possible integer vectors x ε by means of a natural ordering with the index p ∈ N. According to this index, we associate ε x ε with the pth cell ϒ ε p and ε{ x ε } = εy shall denote the local coordinates in all cells which correspond to y ∈ ϒ.
We will denote by ω ε p ⊂ ϒ ε p the respective solid particles obtained by means of the paving with { x ε } = y for y ∈ ω. We note that the rescaling does not change the unit outer normal vector ν. Evidently, the periodic mapping x → y, R d → ϒ, is surjective. This construction can be generalised to an arbitrary orthotope ϒ, see [26].
Let be the reference domain in R d with Lipschitz boundary ∂ and denote again the outer normal vector by ν. By reordering the index p, it is then possible to account for all solid particles ω ε p ⊂ with the index set p = 1, . . . , N ε , see [26,40]. We remark that N ε = O(ε −d ) as ε 0 + . By omitting solid particles which are 'too close' to the external boundary ∂ , we shall ensure a constant gap with the distance O(ε) between ∂ and all particles ω ε p . Thus, is divided into the multiple domains ω ε := ∪ N ε p=1 ω ε p corresponding to all the solid particles located periodically in the reference domain and the remaining porous space \ω ε .
In the following, we shall denote by ∂ω ε = ∪ N ε p=1 ∂ω ε p the union of boundaries ∂ω ε p and introduce the multiply-connected domains Moreover, for functions ξ , which are discontinuous over the interface ∂ω ε , we will denote the jump across the interface by Here, ∂ω + ε = ∪ N ε p=1 (∂ω ε p ) + summarises the positive faces (orientated towards the interior of the pore space \ω ε ), and ∂ω − ε = ∪ N ε p=1 (∂ω ε p ) − accounts for the negative faces (orientated towards the interior of the solid phase ω ε ).

Physical model
In the heterogeneous domain \∂ω ε , which consist of the particle volumes ω ε and the porous space \ω ε , we consider the electrostatic potential φ and (n + 1) components of concentrations of charged particles c = (c 0 , . . . , c n ) , n ≥ 1. The physical consistency requires positive concentrations c > 0.
At the external boundary ∂ , we shall impose Dirichlet boundary conditions φ = φ bath and c = c bath corresponding to a surrounding bath and given by constant values φ bath ∈ R and c bath = (c bath 0 , . . . , c bath n ) ∈ R n+1 + . We can then consider the normalised electrostatic potential φ − φ bath and concentrations c/c bath (i.e. c s /c s bath for all s = 0, . . . , n) and prescribe the following normalised Dirichlet conditions: In the following, all further relations will be formulated for the normalised potential and concentrations such that (1) holds.
Let z s ∈ R denote the electric charge of the sth species with concentration c s for s = 0, . . . , n. For the n + 1-components of charges particles, we shall assume the following charge-neutrality because of the normalisation supposed in (1). A necessary condition for (2) is min s∈{0,...,n} z s < 0 < max s∈{0,...,n} z s . The charge-neutrality is used to justify the Boltzmann statistics when extending the PNP equations to the solid phase with Ohm's law, see equations (5) and Proposition 2 below.
Moreover, the charge-neutrality assumption (2) implies also the following strong monotonicity property for a constant K > 0, which follows directly from Taylor expansion with respect to ( − z s ξ), see [11,Lemma 15]. We consider the following PNP steady-state system consisting of (n + 2) nonlinear, homogeneous equations: In both equations (4), D s ∈ L ∞ ( ) d×d , D s > 0, s = 0, . . . , n denote symmetric and positive definite diffusion matrices, which are in general discontinuous over ∂ω ε . In (4b), κ > 0 is the Boltzmann constant, and T > 0 is the temperature. We remark that the form of (4b) is based on assuming the Einstein relations for the mobilities. Moreover, equation (4a) models the effect of charges particles being included into the solid particles, which is well known, for instance, for Li + -ions, see e.g. [16].
In (5), A ∈ L ∞ (ϒ) d×d denotes the symmetric and positive definite matrix of the electric permittivity, which oscillates periodically over cells according to A ε (x) := A({ x ε }) and satisfies The entries of the permittivity matrix A are discontinuous functions in the cell ϒ across the interface ∂ω. A typical example considers piecewise constant A = σ ω I in ω and A = σ ϒ I in ϒ\ω, with material parameters σ ω > 0 and σ ϒ > 0, where I denotes here the identity matrix in R d×d . In the following, we denote by A ij , i, j = 1, . . . , d, the matrix entries of A. From a physical point of view, (5a) represents Ohm's law in the solid phase. Moreover, we remark that the equations on ω ε , i.e. (4a) for c and (5a) for φ are linear while the equations (4b) and (5b) on \ω ε form a coupled nonlinear problem on the porous space. The modelling of boundary conditions at the interfaces is a delicate issue. For the charge carries fluxes in (4), we assume homogeneous Neumann conditions In fact, imposing inhomogeneous conditions in (7) will lead to quasi-Fermi statistics instead of the Boltzmann statistics, see (16) below. In particular inter-species reaction terms would pose significant difficulties, which are out of the scope of the present work.
For the electrostatic potential in (5), we suppose the following inhomogeneous interfacial boundary conditions (see the related physics in [16]): Here α ∈ R + and g ∈ R are material parameters given at the interface. We note that by summing (8a) and (8b), we derive the relation implying that not only the electric potential φ but also fluxes ∇φ A ε ν are discontinuous functions across the interface ∂ω ε for g = 0.
The asymptotic weights 1 ε in front of [[φ]] and εg at the right-hand side of (8), which were already mentioned in the introduction, shall be discussed in detail during the below asymptotic analysis as ε 0 + . We emphasise that the interface conditions (8) couple the porous phase \ω ε with the solid phase ω ε by means of the jump in [[φ]]. In fact, the jump conditions (8) can be compared with the following two cases of simplified boundary conditions: First, if φ were continuous across ∂ω ε , i.e. [[φ]] = 0, then (8a) and (8b) would be decoupled into two usual Neumann boundary condition which do not represent the EDL. Second, if φ − were known on the solid phase boundary ∂ω − ε , then the model would be reduced to a model on the porous space \ω ε with the following inhomogeneous Robin (called also Fourier) boundary condition (see [13]) However, the subsequent homogenisation of this alternative model on the porous space \ω ε would nevertheless require a suitable continuation of φ + onto ω ε .

Mathematical model
In the following, we shall amend the state variables with the superscript ε in order to highlight the dependency on the cell size.
The following Proposition 2 states the crucial observation that introducing Boltzmann statistics allows to decouple the system of the homogeneous equations (11) and derive an equivalent scalar PB equation. Proposition 2: The system (10)- (12) it is equivalent to the following nonlinear Poisson-Boltzmann equation: find φ ε ∈ H 1 ( \∂ω ε ) such that together with the Boltzmann statistics determining c ε from φ ε , i.e.
We remark that the concentrations c ε in (16) are unique up to fixing the constant positive values within the solid particles ω ε .
By exploiting Proposition 2, we construct a solution (φ ε , c ε ) for the variational problem (10)-(12) from the scalar problem (15) for the potential φ ε . The n+1 concentrations c ε are afterwards explicitly determined by (16). Theorem 1: There exists the unique solution φ ε to the semilinear problem (15) satisfying the following residual estimate which is uniform with respect to ε > 0. Proof: We first emphasise that for the first two terms on the left-hand side of (20) the following discontinuous version of Poincare's inequality for homogeneous Dirichlet condition (15a) holds on the multiple domains \∂ω ε without interfaces ∂ω ε (see e.g. [34,36]): Therefore, the lower estimate (21) together with (3) ensures the coercivity of the operator of the problem (15b). The main difficulty of the existence proof arises from the unbounded, exponential growth of the nonlinear term in (15b). While classic existence theorems on quasilinear equations are thus not applicable here, the solution can nevertheless be constructed by a thresholding, see e.g. [17] and references therein for the details.
To derive the estimate (20), it suffices to insert φ = φ ε as the test-function in the variational equation (15b) and apply (3) in order to estimate below the nonlinear term at the left-hand side of (15b). Finally the right-hand side of (15b) can be estimated by means of the following trace theorem see [30] for the details. This completes the proof. We remark that in the following Section 3, we will refine the residual error estimate (20) by means of asymptotic analysis as ε 0 + and homogenisation.

Homogenisation and residual error estimate
We start the homogenisation procedure with three auxiliary cell problems. The first two cell problems serve to expand the inhomogeneous boundary traction g and the volume potential of the variational problem (15) from the porous space \ω ε onto the whole domain \∂ω ε . The third cell problem is needed to decompose the matrix A ε of oscillating coefficients in the cells with respect to small ε 0 + . This procedure will result in a regular asymptotic decomposition of the perturbation problem with a subsequent error estimate of the corrector term.
For a generic cell ϒ, we introduce the Sobolev space H 1 # (ϒ) of functions which can be extended periodically to H 1 (R d ). This requires matching traces on the opposite faces of ∂ϒ. Moreover, we shall denote by H 1 # (ϒ\∂ω) those periodic functions, which are discontinuous, i.e. allow jumps across the interface ∂ω.

Auxiliary results
We state the first auxiliary cell problem as follows: find L ∈ H 1 (ϒ\∂ω) such that In view of the homogenisation result stated in Theorem 2 in Section 3.2 below, the auxiliary problem (23) serves to expand the inhomogeneity of the boundary condition (8a) given by the material parameter g in terms of the weak formulation stated in (15b). The existence of a unique solution L in (23) follows via standard elliptic theory from the assumed properties (6) of A. With its help, we are able to prove the following result. ( 2 4 ) where l 1 : H 1 ( \∂ω ε ) → R is a linear form satisfying

Lemma 1: (The cell boundary-traction problem) For all test-functions
Proof: We apply the auxiliary cell problem (23). By inserting a constant test-function u, we calculate the average value Here, |∂ω| and |ϒ| denote the Hausdorff measures of the solid particle boundary ∂ω in R d−1 and of the cell ϒ in R d , respectively. Subtracting ϒ\∂ω L y u dy from (23), we rewrite it equivalently as where we have added to the residuum l(u) the trivial term In the following, we shall apply the discontinuous Poincare inequality and the Trace Theorem with K 2 > 0, which combine to the estimate By recalling that A ∈ L ∞ (ϒ) d×d and by applying Cauchy's inequality to the right-hand side of (27) and subsequently applying estimate (30) to L and u, we obtain the following estimate with K from (6) and K 3 from (30). For a proper test-function φ(x) with x = ε x ε + ε{ x ε }, we insert u(x, y) = φ(ε x ε + εy) into (27) and apply the periodic coordinate transformation y → x, ϒ → R d , by paving R d such that { x ε } = y (recall Section 2.1). After observing that dy → ε −d dx, dS y → ε 1−d dS x , ∇ y → ε∇ x , we also multiply (27) with the constant gε d and use (26) which is (24) with the following right-hand side term: where we denote L ε (x) . Similarly, (30) transforms into the uniform estimate with K 3 > 0 from (30). We note that the first line of (33) expresses the H 1 -norm by the standard homogeneity argument, see e.g. [22, Appendix, Lemma 1, p.370]. Therefore, the estimate (31) of l yields the following estimate of l 1 Here we used (6) and inequalities (30) for L and (33) for φ. Then, (34) follows (25) with the constant K = |g|(K + K 2 3 ) L H 1 (ϒ\∂ω) , which completes the proof. Remark 1: We remark that Lemma 1 justifies not only the a-priori estimate (22), but also refines it by specifying the limiting asymptotic term as ε 0 + , which consists of the constant potential |∂ω| |ϒ| g distributed uniformly over .
Another idea for the proof of Lemma 1 can be referred in [30,33]. The next auxiliary cell problem studies the asymptotic expansion of a volume force f ∈ H 1 ( \∂ω ε ), which is given on the porous space ϒ\ω surrounding the solid particle ω ⊂ ϒ. It will be applied in particular to the nonlinear term in (15b), i.e. we shall consider the specific volume force f (x) = − n s=0 z s exp (− z s κT φ 0 (x)) in Theorem 2 below. With x = ε x ε + ε{ x ε } (recall Section 2.1), the following unfolding operator is well defined, see [26]. For its modification near the boundaries ∂ of nonrectangular domains , see [40].
where l 2 : H 1 ( \∂ω ε ) → R is a linear form satisfying Proof: By inserting a constant test-function u into the auxiliary cell problem (35), we obtain the locally averaged value of M = M(x, y) Moreover, by using the average T ε f y , we can expand See [41] for the analysis of expansion (39) in terms of Fourier series. For fixed x the residual F(x, y) has zero average F y = 0 and estimates as due to the discontinuous Poincare inequality (30). By inserting (39) into (38), we calculate and thus derive by using again (39) After multiplying the identity (41) with u and integrating it over ϒ\∂ω, we subtract it from (35) and rewrite (35) equivalently as where we have added the trivial term ϒ\∂ω (M − M y ) u y dy = 0 and the residuum m(u) shortly denotes the right-hand side terms of (42).
For fixed x ∈ \∂ω ε , Cauchy's inequality yields for the first term on the right-hand side of (42) Thus, by applying the estimates (40) and (43) to F and the discontinuous Poincare inequality (30) to M and u, we estimate m(u) at the right-hand side of (42) as where K 4 = |ϒ\ω| |ϒ| 1 + |ϒ| |ϒ\ω| K 3 and by recalling K from (6) and K 3 from (30). Next, we substitute u = (T ε φ) as the test-function in (42) and use the property T ε f ·T ε φ = T ε (f φ) of the unfolding operator. After applying the periodic coordinate transformation y → x, { x ε } = y to (42) similar to the proof of Lemma 1, we arrive with T ε (f φ) → f φ and T ε φ → φ at (36) with where hence, where we have used (30) for M(x, · ) and (33) for f and φ. Thus, (46) implies the estimate (37) of the residual term l 2 given in (45) with This completes the proof.
The third cell problem considers the solutions of the following system of d linear equations: find a vector of periodic functions N = (N 1 , . . . , N d Here, H 1 # (ϒ\∂ω) denotes the space of periodic H 1 -functions and DN(y) ∈ R d×d for y ∈ ϒ\∂ω stands for the row-wise gradient matrix of the vector N, that is Moreover in (47), Dy = I ∈ R d×d yields the identity matrix. The solvability of (47) follows from the symmetry and positive definiteness assumption (6). The uniqueness of the solution N is provided due to the constraint N y = 0. Indeed, since N(y) + K with an arbitrary constant K solves also (47), the zero average condition is sufficient (and necessary) to ensure the uniqueness of the solution, see e.g. [21]. Finally, the solution is smooth locally in ϒ\∂ω. The system (47) is essential to determine the efficient coefficient matrix A 0 of the macroscopic model averaged over . In fact, following the lines of [21,23], we shall establish an orthogonal decomposition of Helmholtz type for the oscillating coefficients A ε .
The Helmholtz type decomposition is based on the left-hand side of (47) defining an inner product · , · in H 1 # (ϒ\∂ω). Due to [[y]] = 0, the variational equation (47) reads as N + y, u = 0 for all u ∈ H 1 # (ϒ\∂ω), which implies that N + y belongs to the kernel of this topological vector space. Thus, the fundamental theorem of vector calculus (the Helmholtz theorem, see e.g. [23]) permits the following representation as sum of a constant matrix A 0 and divergence free B(y) fields in R d×d : where B has zero average, i.e. Thus, we obtain the following lemma: Lemma 3: (The cell oscillating-coefficient problem) The constant matrix of effective coefficients is determined by averaging Moreover, A 0 is a symmetric and positive definite matrix with the entries: For the transformed solution vector N ε (x) also depends only on { x ε }, the following decomposition holds: The symmetry and positive definiteness of A 0 follow straightforward from the assumption in (6) Finally, we apply the periodic coordinate transformation y → x, ϒ → R d , with y = { x ε } to (48) and (55). With ∇ y → ε∇ x , we have for the row-wise gradient matrix D y N → εD x N ε and B → εB ε . Thus, we arrive at (51) and (54). The proof is completed.

The main theorem
Based on the Lemmata 1-3, we formulate the main homogenisation result: For smooth solutions φ 0 and N of (56) and (47), in the limit ε 0 + , the solution φ ε of (15) and the corrector term φ 1 := φ 0 + ε(∇φ 0 ) N ε to φ 0 satisfy the residual error estimate (improving (20)): Proof: First, we remark that the left-hand side of (57) defines a norm in H 1 ( \∂ω ε ) due to the lower estimate (21). Secondly, the unique solution φ 0 of (56) can be established by following the arguments given in the proof of Theorem 1. Moreover, the solution is smooth inside by standard arguments of local regularity of weak solutions, see e.g. [22,Appendix] and references therein.
Next, we prove the residual error estimate (57). Integrating (56) by parts on yields the strong formulation −div A 0 ∇φ 0 − |ϒ\ω| |ϒ| n s=0 z s e − z s By applying the Green formulas (13a) and (13b) in \ω ε and ω ε , respectively, we have for all and for all φ ∈ H 1 (ω ε ): By summing these two expressions and by using the continuity of ∇φ 0 across the interface ∂ω ε , we insert the strong formulation (58) into the above right-hand sides and rewrite problem (56) in the

Discussion
In the following, we shall summarise the main observations concerning the presented results.
• We observe that the first two terms on the right-hand side of (71) express the residual error near ∂ and at ∂ω ε . These terms are asymptotically of order O( √ ε) (as can be see by setting t 1 = O(ε −1/2 ) = t 2 in (75) and (76)) and thus constitute the leading order O(ε) in the residual error estimate (57). Therefore, by constructing corrector terms in form of the respective boundary layers, which is possible in practice for polyhedral domains with rational slopes of their flat boundaries with respect to the unit cell, the O(ε)-estimate (57) could be improved to the order O(ε 2 ).
• The factor 1 ε appears at the jump across interface ∂ω ε in the left-hand side of microscopic equation (15b). It is controlled by the coercivity condition (21). In the homogenisation limit, this term disappears in the macroscopic equation (56). Instead it rather contributes to the macroscopic coefficients in (50).
• The factor ε in front of the inhomogeneous material parameter g, which is prescribed at the solid phase boundary ∂ω − ε , presents the critical order. After averaging this factor guarantees the presence of the potential |∂ω| |ϒ| g distributed over the homogeneous domain in (56).