Asymptotic expansions for the conductivity problem with nearly touching inclusions with corner

We investigate the case of a medium with two inclusions or inhomogeneities with nearly touching corner singularities. We present two different asymptotic models to describe the phenomenon under speciﬁc geometrical assumptions. These asymptotic expansions are analysed and compared in a common framework. We conclude by a representation formula to characterise the detachment of the corners and we provide the possible extensions of the geometrical hypotheses.


V IR GINIE B O N N A IL L IE -N O Ë L C L A IR P OIG N A R D G R É G O R Y V I A L A S Y M P T O TIC E X PA N SIO N S F O R T H E C O N D U C TI V IT Y P R O B L E M W IT H N E A R LY T O U C HIN G IN C L U SIO N S W IT H C O R N E R D É V E L O P P E M E N T A S Y M P T O TIQ U E P O U R L E P R O B L È M E D E C AV IT É S U R U N D O M A IN E C O N T E N A N T D E S IN C L U SIO N S AV E C C OIN S P R O C H E S
Abstract. -We investigate the case of a medium with two inclusions or inhomogeneities with nearly touching corner singularities. We present two different asymptotic models to describe the phenomenon under specific geometrical assumptions. These asymptotic expansions are analysed and compared in a common framework. We conclude by a representation formula to characterise the detachment of the corners and we provide the possible extensions of the geometrical hypotheses.

Introduction
The aim of the paper is to investigate the asymptotic behavior of the solution to the conductivity (also called thermal) problem in a domain with two nearly touching inclusions with corner singularity. It is well-known that behavior of the electrostatic field near a corner of an inclusion depends on the angle of the corner, see [Gri85] for instance. For acute angles, from the domain of computation, the electric field remains bounded but as soon as the angle is larger than π, the electric field blows-up.
Problems with inclusions have given rise to an important literature, in various contexts, mainly motivated by mechanical and electrical engineering, or image processing. In the special case of smooth inclusions close to each other, or close to the boundary of the domain, we can mention for example the works [AKT06, BLY09, BT13, BTT18, KLY13,NP18]. For nonsmooth inclusions, we can refer to [MN86,NT18].
The configuration we are looking at here is quite trickier since we consider a domain with two inclusions with corner singularity at distance δ 1 from each other. Depending on the geometrical configuration, the corner asymptotics of the solution to the conductivity problem at δ = 0 may be dramatically different from the asymptotics of the solution at δ = 0. In this paper, we provide a common framework for two different methods to derive the asymptotic expansions of the solution at δ = 0. These expansions are mainly based upon the solution to the limit problem δ = 0 and the appropriate so-called profile terms. These profiles are solution to the conductivity problem in a sectorial domain with appropriate source terms. They appear naturally from the expansion, in the same vein as in [CCDV06].
Precisely, we are interested in the model problem where g is a given datum (trace of a smooth function), and the unbounded domain Π δ is defined by Π δ = R 2 \ Ω L δ ∪ Ω R δ , the domains Ω L δ and Ω R δ standing for the inclusions. The distance between Ω L δ and Ω R δ is 2δ, see Figure 1.1. Naturally, the formal limit case where the inclusions touch each other corresponds to The domain Π δ with the inclusions Ω L δ , Ω R δ , and the corresponding limit domain Π 0 .
Nevertheless, the asymptotic model is not unique since the dependence of the domains Ω L δ , Ω R δ with respect to the parameter δ is not explicitly given by the physical situation. We will investigate the following two choices: • Translation method. The size of the inclusions Ω L δ and Ω R δ is supposed to be independent of δ. Changing the value of the parameter δ merely corresponds to a translation of the inclusions along a given unit vector d. The subsequent asymptotic expansion is detailed in Section 3.1 where the domain Π δ will be denoted by • Π δ for notation consistency.
• Contraction method. The inclusions do not have a constant size with respect to δ, but are contractions of the limit case (i.e. δ = 0). This case is precisely addressed in Section 3.2. In particular, the geometrical setting will be made clear in Subsection 3.2.1, where, here again, the domain Π δ will be denoted Π δ . If the two presented asymptotic models coincide for δ = 0, both have pros and cons since they can handle different geometrical frameworks of inclusions. For comparison purposes, which yield one of the main goal of this paper, we restrict the geometrical configurations to a situation where both models might be considered. The precise framework is detailed in the following assumptions, but the reader can keep in mind the situation of By convention, we fix the horizontal axis along this direction. Moreover, Ω L 0 and Ω R 0 are supposed to be included into R − × R, and R + × R, respectively. Remark 1.2. -Replacing in Assumption (H1) one point by a finite number of points generates only technical difficulties. We do not address the case of inclusions without corners (for example two disks), which generates a cuspidal limit domain. Hypothesis (H2) prevents from this situation. Besides, the regularity of the boundaries far from the origin is not necessary and is only assumed here for simplicity.
Likewise, Hypothesis (H3) is technical and might be weakened. The relaxation of theses hypotheses will be discussed in Section 5.
Section 2 recalls the decomposition in regular and singular parts of the Laplacian equation (1.2) with adding a source term. Then we consider the two asymptotic models. First in Section 3, we present a formal asymptotic expansion where we do not take care of the corner singularities. Even though formal, this construction is useful for the complete analysis and enables to deal with difficulties step by step. We present an asymptotic model for the translation and contraction case and propose a unified formulation in Section 3.3. In Section 4, we take the singularities into account, we construct the corner profiles and we give the complete asymptotic expansions obtained by the two methods.
As far as we know, the comparison of these two asymptotic methods has not been addressed before. In this paper the two methods are studied in parallel and presented in a unified formulation which allows us to deduce a generic full asymptotic expansion of the solution of (1.1).
For δ > 0, the solution u δ of problem (1.1) and the limit solution u 0 both have singularities at the origin, but with different singular exponents because they are associated with different angles. In the geometry of Figure 1.1, u δ has singularities due to the re-entrant corners whereas u 0 has only singularities of the convex corners (of angle α ± , see Figure 1.1). The asymptotic procedures that will be described in the paper will allow to compare u δ and u 0 . Precisely, the first terms in the asymptotic expansion read where • The functions ψ and χ are smooth radial cutoff functions satisfying ψ(x) = 0 and χ(x) = 1, near 0.
• The coefficients c ± 0,0 are the first singular coefficients of u 0 , solution to (1.2). • The profiles K ± 0 and M are defined in the infinite domain P (see Figure 4.1) as the respective solutions to Problems (4.6)-(4.7).
Let us give some indication on the magnitude of the profiles with respect to δ: As we can see in (1.3), the singularities of u 0 are canceled through the cutoff function ψ (at distance O(δ) from the origin), and replaced by adapted counterparts through the profiles terms, rescaled to fit the δ-depending domain Π δ . It is worth noting whatever the method, the zero th -order terms of the expansion of u δ are given by (1.3). The differences in the expansions will appear at the next order terms. The goal of the paper is to push forward the expansion for each method to provide different ways to approximate u δ .

Corner asymptotics for the limit problem
Let us consider now the solution u to the generic problem: where f and h are smooth given data (for f = 0 and h = g the solution u is nothing but u 0 ). Since corners appear in the limit domain Π 0 , singularities arise for the solution u near the origin, preventing a full elliptic regularity. Precisely, if α ± denote the absolute value of the sector angles involved in Π 0 (see Figure 1.1) one has: Notice that according to (H3), α ± ∈ (0, π).
With Assumption (H4), the following decomposition (2.2) is valid. Otherwise, logarithmic singularities may appear in the corner asymptotics. We do not consider this case here for simplicity.
The standard theory of corner problems applies here (see [Dau88,Gri85,Kon67] for example). The potential u admits a splitting into regular and singular parts for any integer p: where, using polar coordinates (r, θ), • u p flat ∈ H p loc (Π 0 ), and for any multi-index β with |β| p, we have ∂ β u p flat (x) = O(r p−|β| ), near the corner 0, (2.3) • χ : R 2 → R + is a smooth radial cutoff function, with compact support and equal to 1 near 0, • we denote by • the coefficients c ± λ ∈ R depend on f and h and are called the singular coefficients of order λ of u in the upper or lower part of Π 0 . They can be numerically computed from the solution u via extraction formulae (see [Gri85,Chapter 8 • the singular functions s ± λ are defined for any λ ∈ Λ ± p by otherwise, see Figure 1.1 for the definition of θ ± 1 and θ ± 2 . Note that for λ = 0, the functions s ± 0 equal 1 for θ ± 1 < θ < θ ± 2 . • The last series in (2.2) comes from the data f and h, which do not necessarily vanish near the corner. The coefficients a ±,σ n are explicitly determined by the Taylor expansions of f and h. The functions a ±,σ n are defined for any positive integer n and any σ ∈ {1, 2} by Let us comment the dependence with respect to p. The bigger p, the more terms in the singular expansion, and the flatter the regular terms.
Before dealing with singularities, it is important to derive the asymptotic expansion in the ideal case where no singularity appears at any order. Actually, understanding how such a derivation occurs will facilitate the accurate asymptotic expansion with the singularities.

Asymptotic models without singularities
The goal of the section is to present the formal derivation of the asymptotic expansions for both translation and contraction cases, without accounting for the influence of corner singularities. We emphasize that throughout this section, we assume that no singularity appears in the derivation of the asymptotics, meaning that all the coefficients c ± λ and a ±,σ n of any expansion of type (2.2) equal zero. This implies in particular that the datum g is flat at any order near the origin. Let us mention that this assumption does not hold only for the limit term u 0 , but also for all terms which will be constructed along the asymptotic expansion (see Equation (3.16) for instance).
This ideal case has no chance to hold in general. Even the flatness of the data f and h are not enough to ensure that no singularity appears in the solution of one particular problem (2.1). Moreover during the procedure, the data will involve complex combinations of the previous terms. It is then not possible to exhibit a nontrivial case of such a situation (even if such a situation might theoretically happen). The present section mainly has a pedagogical role.
Under the assumption that no singularity appear, the expansion that will be derived below only involves coefficients with entire powers of δ. In the general case, singularities will introduce non entire powers of δ, and they will lead to a modification of the source terms of the following elementary problems, but the procedure will remain unchanged. Therefore even though very specific, this ideal framework is interesting for understanding the way the cascade of elementary problems is derived at any order to get the asymptotics. We derive such "ideal" asymptotics, in the two different geometrical frameworks. In addition, we unify both approach but introducing appropriate notations, making thus possible a comparison of the two methods. The difficulties arising from the corner singularities are postponed to Section 4 where the corner profiles are introduced.

Asymptotic model with translation of inclusions
This section is devoted to the construction of the asymptotic expansion for Problem (1.1) in the case where the inclusions are constant patterns translated from each other by a distance of 2δ. As mentioned previously, to distinguish the two geometrical frameworks, the notation will are enriched by a • for the translation case.

Geometrical setting
We assume here that the domains • Ω L δ and • Ω R δ are translations of constant size inclusions : It is classical to transform the problem into a domain independent of δ in order to make the parameter δ appear in the equations. In the present case, the "gap" between the inclusions can be seen as a thin layer. It is then natural to set the change of variables Through this change of variables, the domain Figure 3.1. Here, we obviously have Considering Problem (3.2), we postulate the two following Ansätze Note that the series are written in the sense of asymptotic expansions (i.e. truncate the series and make δ → 0, see [MNP00] for example), we do not expect (nor need) convergence as n → ∞. Inserting these Ansätze into Problem (3.2), we get the following cascades of elementary problems: with the convention • u ext = 0 and • u lay = 0 for < 0. In the first two equations of the problem, solved by • u lay n , the variable τ is actually only a parameter. The question of the solvability of this sequence of problems will be discussed in Section 3.3.

Asymptotic model with contraction of inclusions
As mentioned previously, the way each asymptotic expansion is derived is somehow arbitrary, and it depends mostly of the point of view of the modeling. In Section 3.1, we considered the case of two inclusions with δ-independent size, whose corners are distant from 2δ.
In this section we adopt another point of view, leading to the same limit geometry: we assume that the inclusions are contractions of the limit inclusion Ω L 0 ∪ Ω R 0 given by Figure 1.1. Let us first make precise the geometrical assumptions. As mentioned previously, to distinguish the two geometrical frameworks, the notation is enriched by a for the contraction case.

Geometrical setting
Identification of the thin layer. For δ small enough, we define here Ω R δ and Ω L δ as contractions in the normal direction ν of the respective limit inclusions Ω L 0 and Ω R 0 . A thin layer Π lay δ appears naturally, when passing from the domain Π δ to the limit domain Π 0 := Π 0 (see Figure 3.2).
We denote respectively by Π ext , Π lay δ and Π δ the 3 domains (see Figure 3.2). Let us introduce a suitable parameterization of ∂ Π δ . In the interior of the ball B ρ = {|x| < ρ }, the boundary ∂ Π δ is deduced by translation by ±δd from the right and left parts of ∂ Π 0 .
we choose a smooth, positive and δ-independent function b to parameterize ∂ Π δ : It is designed in such a way that the two definitions are consistent in the intersection { ρ 2 |x| ρ }. TOME 3 (2020) Remark 3.1 (Arbitrary choice for b near 0). -The choice of ρ is somehow arbitrary, but it does not change dramatically the derivation of the expansion, since the function u 0 is supposed to be flat near the origin 0, as stated by (2.2)-(2.3).
Remark 3.2 (The case of constant thickness). -In the simple geometrical framework where the plane sectors are equally distributed from both side of the horizontal axis (0, d), the simplest choice for b is the constant function defined by Once the identification and the parameterization of Π lay δ is performed, it is natural to introduce new variables (η, τ ), which respectively correspond to the normal and tangential variables along ∂Π 0 , in order to rewrite the Laplace operator in these new coordinates. We parameterize ∂Π 0 by the vector field Ψ defined on the torus T: By abuse of notation, we still denote by b and ν the functions b • Ψ and ν • Ψ defined on T. We then define the mapping Φ δ from (0, 1) × T into R 2 as Note that Φ δ applies only to the layer part (contrarily to • Φ δ which is defined in R 2 ). In local coordinates (η, τ ), the Euclidean metric (G ij ) i,j∈{1,...,2} is given by Denoting by |G| its determinant: 12 , Laplace's operator reads then and thus it can be expanded as where T k are differential operators independent of δ and κ denotes the curvature (extended by 0 at the origin). On the other hand, the normal derivative is rescaled ANNALES HENRI LEBESGUE by a factor 1/(δb): where B k are differential operators independent of δ.

Elementary problems for the contraction case
By analogy with the translation case, we denote by Γ trans the interface across which transmissions occur. In the contraction case, Γ trans is nothing but ∂Π 0 : We denote respectively by u ext,δ and u lay,δ the function u • Φ −1 δ in Π ext and Π lay . After this change of variables, the elementary problems read where ξ → g n (ξ) are obtained from the Taylor expansion of g • Φ −1 δ (ξ). Polynomial Ansätze similar to the translation case, lead to the following sequence of elementary problems: u ext n → 0 at infinity, Again, we postpone the discussion of the well-posedness of this sequence of equations to Section 3.3.

A unified formulation
In Sections 3.1 and 3.2 we have derived an asymptotic expansion for two different asymptotic families converging to the same limit problem. Actually, it is possible to describe a general framework containing the two procedures. Once again, we emphasize that in this section, we are dealing with the ideal cases where the singularities are omitted at any order. Section 4 is devoted to derive the sector profile problems to account for the singularities, that appear in the general case.
Geometrical setting. We consider the problems obtained after change of variables (solved by • u δ and u δ , respectively). We denote by Π the domain where the problem is set, it can be split into where Γ trans = ∂Π ext ∩ ∂Π lay . The layer Π lay has the form (η − , η + ) × I. Finally, we set Γ = ∂Π \ Γ trans .
• In the contraction case, η − = 0, η + = 1, I = ∂Π 0 , and Γ trans = ∂Π ext = {0} × I. For the sake of conciseness, for any function u defined in Π we naturally denote by u ext and u lay its restrictions to Π ext and Π lay respectively.
Layer problem. We consider the following equations (3.10) where the term T k u lay n+1−k is given explicitly by u lay for 0 n. Note that in the translation case, only the term T 2 = ∂ 2 τ does not vanish. The expression of k n depends on the geometrical framework: • In the translation case, k n = ∂ ν u ext n | η=η + . • In the contraction case, k n = b g n − n k=1 B k u lay n−k . Problem (3.10) is exactly obtained from (3.9) in the contraction case whereas the two following equations are missing in the translation case: for η = η + , ∂ ν u ext n = ∂ η u lay n+1 for η = η − . They will be taken into account in (3.16). In order to solve Problem (3.10), let us now define the operator R for any function j as

ANNALES HENRI LEBESGUE
Thanks to (3.10), one easily obtains the formal expressions ∂ η u lay n+1 (η, τ ) = k n (τ ) + R[j n ](η, τ ), (3.14) Exterior problem. In the outer domain Π ext the problem at the order n has the following form: (3.16) Find u ext n such that n → 0 at infinity, whereγ, φ n and g n are defined below.
• In the translation case, equations (3.12) together with (3.14)-(3.15) make appear the vector operatorγ defined bȳ and the data given by On the other hand, the Neumann datum g n is given by the Taylor expansion of g (see (3.3)) g n = • g n . • In the contraction case ∂Π ext \ Γ trans = ∅ thus no need to define g n . The operatorγ is scalar and given bȳ γ(u) = ∂ ν u| Γ trans , and the datum φ n is given by Sequential derivation of the terms. At order n = 0, the solution to Problem (3.10) is formally given by Formulae (3.14)-(3.15): where k 0 = ∂ ν u ext 0 | η=η + in the translation case and k 0 = b g 0 in the contraction case. The exterior Problem (3.16) is thus completely determined since • In the translation case, φ 0 = (0, 0), and no • g 0 is needed, • In the contraction case, φ 0 = 1 b ∂ η u lay 1 | η − = g 0 , which thus entirely determines a posteriori the terms u lay 0 and ∂ η u lay 1 . The general procedure works as follows (each situation is detailed for the sake of clarity).
Algorithm A (Translation case). -Assume that (u ext , u lay ) are known for any n − 1.
Thus, the datum φ n is well defined, see (3.17). • We can determine u ext n by solving Problem (3.16), since φ n is known, and g n only depends on the Taylor expansion of the Neumann datum g.
• Knowing u ext n , we compute u lay n with Formula (3.15), which is recalled here Algorithm B (Contraction case). -Assume that (u ext , u lay ) are known for any n − 1.
• We compute the following terms, involving u lay for n − 1 only.
• Problem (3.16) can be solved to compute u ext n since φ n is known, see (3.18). • Knowing u ext n , we compute u lay n with Formula (3.15), which is recalled here

ANNALES HENRI LEBESGUE
The above derivation process can be summarized as follows where J ext n and J lay n are the operators defined by the above problems (3.16) and (3.10) respectively. We will make use of J lay 0 as well, setting u lay −1 = 0 by convention.
Remark 3.3. -In the layer, the problem defining the terms u lay n is one-dimensional, and the tangential variable τ only acts as a parameter. The integration is only performed in the normal variable η. It is straightforward that the dependence on u lay n is polynomial in η, its degree being increased by 1 at each step. The dependence with respect to τ is more complex since it comes from traces on η = η ± of the terms u ext for < n.

Back to the physical domain
Remember that the above procedure is based on the assumption that all the coefficients c ± λ and a ±,σ n vanish in the expansions of type (2.2), not only for the limit term u 0 but also for u lay n with arbitrary n. This implies that at any order, all the terms u ext n and thus u lay n are flat near the origin, see Section 2. It is however important to transfer the obtained approximate solution in the original domain.
• Limit domain δ = 0. In the contraction case, the first term u ext 0 previously determined coincides with the limit term u 0 defined by (1.2).
In the translation case, the domain • Π ext , where u ext 0 is defined, is not connected. However the functionũ ext 0 defined in Π 0 by In the translation case, an approximation at order N is given by In the contraction case, approximations far from the corner have been derived thanks to the parameterization. In the vicinity of the corner, singularities should be dominant. However since we have assumed throughout Section 3 that no singularity appears, this means that at any order the functions are flat enough. Precisely, our hypothesis states that u δ is O(|x| p ) for any p ∈ N. This implies that 0 is a good approximation of u δ near the corner. Therefore an approximation of u δ at the order δ N is given by where the mapping Φ δ is defined in (3.5).

Handling singularities 4.1. Need for a specific procedure
As mentioned before, the procedure described in Section 3.3 fails to provide an asymptotic expansion of u δ when singularities appear. We detail here the reasons why.
First obstruction: lack of regularity for the sequential derivation. The formal derivation described in Section 3.3 fails since singularities appear in the solutions to (3.16). To explain this point, we distinguish the two configurations.
• Translation case. In order to solve Problem (3.16) in a variational framework, it is necessary that φ n has regularity H 1/2 ( In Algorithm A, the condition over the first component needs generically k n−1 ∈ H 1/2 ( Thus we need that has limited regularity (due to the presence of singularities s ± λ ), an expansion at any order is hopeless (if the angles α ± are close to π, the maximal order for the construction may be n = 1).
Second obstruction: terms are not necessarily flat near the origin. We assume here that the first obstruction occurs for the construction of u ext 2 . This arises if u ext 1 / ∈ H 2 (Π ext ). Precisely, according to (2.2) for p = 1, the following splitting holds for the exterior part with u 1,ext flat,1 ∈ H 2 (Π ext ), s ± λ / ∈ H 2 (Π 0 ), a ±,σ 1 ∈ C ∞ (Π 0 ). Therefore, if the coefficients c ± λ,1 vanish (again, this is not the general situation), we have u ext 1 ∈ H 2 (Π ext ), which seems sufficient for the procedure. However the terms a ±,σ 1,1 , which come from the Neumann traces, add new difficulties.
• Translation case. In order to build u lay 2 , it is necessary to derive u lay 0 (η, ·) = u ext 0 | Γ trans twice with respect to τ , which is impossible in general, since a +,σ 1,1 = a −,σ 1,1 . Introducing a truncation near the origin makes it possible to prevent this case, similarly to the contraction case. However this makes appear a non negligible error near 0.
• Contraction case. Here, u ext 2 can be built. Nevertheless, we do not end up with a suitable asymptotic expansion. Indeed, coming back to the physical domain (see Section 3.4), we get a correct approximation only away from the origin because of the extension by 0 near 0.
Third obstruction: error estimates fail. Should the above obstructions be overcome to define u 1 , a third issue would appear for the error estimates. Actually, the terms u δ − u [1] δ ψ(|x|/δ) cannot be used to obtain error estimates by computing the residual, since ∆u 1 is not well-defined in Π lay as a function.
How to overcome these obstructions? We have seen that the problem comes from the limited regularity of u ext 0 . We will tackle this issue by inserting only the regular part of the terms (the "flat" terms in expansions of type (2.2)) in the algorithmic procedure described in Section 3. The singularities will be handled in a different way (and then stay unaffected by the change of variable Φ δ ), with the introduction of profiles, see Section 4.2.

Heuristics for the introduction of the first profile
We have seen just above that the flat terms can be introduced in the procedure exposed in Section 3.3. We now explain how to take the singularities into account.
Let us go back to the approximation at the order 0. We set δ is given in Section 3.4. Here ψ is a smooth radial (1) cutoff function independent of δ and equal to 1 far from 0 and which vanishes near 0.
The remainder term r δ is well-defined in Π δ and satisfies the following transmission problem: where u 0 is the solution to (1.2). Note that in the translation case, the operators B k are zero. The term [∆; Ψ] generically refers to the commutator operator given by for smooth enough Ψ and V , and thus the term [∆; ψ( · /δ)] u 0 is of order δ −2 . In Problem (4.5), the main source terms involve the contribution of u 0 and g 0 both taken near the origin 0, since ψ( · /δ) equals 1 far from it. Splitting u 0 near the origin into regular and singular parts, see (4.2)-(4.3), we see that the principal part of the error is given by s ± 0 . It is natural to introduce the profile in the infinite domain P = lim δ→0 δ −1 Π δ , see Figure 4.1, defined as (4.6) (1) By abuse of notation, we use "radial" in the layer as well, but in the translation case, this corresponds to a dependence in η.
where s ± 0 have been extended by a constant on the support of ∇ψ. Moreover, the Neumann boundary condition also generates a profile denoted by M: The well-posedness and the behavior at infinity of K ± 0 and M are discussed in Section 4.4. These profiles make it possible to correct the expansion as follows Thanks to (4.6), the remainder R δ satisfies the following problem: The leading term of c + 0,0 s + 0 + c − 0,0 s − 0 − u 0 near the origin is supported by (a ± 1,i ) i=1,2 according to the splitting (4.2) of u 0 , and it is thus of order δ. The term vanishes inside a ball of radius O(1) centered at the origin, and thus its leading term is given by the behavior of the profiles K ± 0 and M at infinity, that is O(δ min(λ + 1 ,λ − 1 ) ), where λ ± 1 = π/α ± > 1. Therefore R δ has a smaller residual than r δ , and the expansion is improved.
In order to build the terms of order δ, it is necessary to account for three contributions • The next term in the splitting of u ext 0 which involves the singular functions a ±,σ 1 , see Section 2. • The Neumann source term (1 − ψ( · /δ)) g 0 near the origin.
• The termǔ 1 defined in (4.4), and in particular its first singular terms, involving the singular functions s ± 0 . The third contribution is addressed by the profile K ± 0 defined by (4.6), while the first two contributions will be handled by new specific profiles L ±,σ 1 defined as: The expansion at the next order reads then where c ± 1,0 are the first singular coefficients ofǔ ext 1 in the splitting near the origin. The remainder estimate comes from the first neglected term in the construction, that is Actually, to obtain the estimate in (4.10) requires a bit more work (see [CCDV06,Section 4.2] for more details on such estimates). Let us point out that the exponents λ ± 1 involved in this remainder are associated with the limit problem. They differ from singular exponents which arise in the profile problems, corresponding to the reentrant angles of the domain P. These additional exponents will actually not affect the δ-expansion since only the profile behavior at infinity is involved thanks to the rescaling x/δ, but the associated singularities impact the profiles themselves.
The above expansion has a priori an order of accuracy in O(δ min (λ + 1 , λ − 1 )). To go one step further, it is necessary to change the whole expansion, and to account for one more singular term in the splitting of u ext 0 , which modifies u 0,ext flat,0 and thusǔ 1 . So it is necessary to determine a priori the desired order of accuracy, so that u ext 0 be split with the appropriate number of singular functions, and then the expansion can be performed.

Structure of the complete asymptotic expansion
In the previous section, we have described the construction of the first terms in the expansion, taking into account the singularities coming from the corners of the domain. It is possible to go further and build a complete asymptotic expansion at any order. To do so, we need to fix a target order P , and expand every term defined in Π ext into singular and regular parts according to this choice. The contribution of all singular functions arising in such splittings will be handled through profiles (of type K, L, M) defined in the infinite domain P. The expansion of these profiles at infinity will give birth to other terms defined in Π ext , which in turn will be split into regular and singular parts.
where • the terms u P,ext µ are defined in the standard variable x inside Π ext , • the terms W ±,σ,P µ are defined on the infinite domain P, and are involved in the asymptotic expansion as functions of the scaled variable x δ , • the set of indices M P is given by ; p, q ± ∈ N 0 , and µ P , • for any compact subset ω ⊂ Π δ , the remainder R ext δ,P satisfies R ext δ,P H 1 (ω) = O(δ P ). In the layer, the expansion has a similar structure, but involves the normal dilation used to build the terms : The remainder R lay δ,P satisfies the following estimate on any compact subset ω ⊂ Π δ : √ δ R lay δ,P H 1 (ω) = O(δ P ). Remark 4.2. -If P is changed, the terms involved in the expansion changes as well, since the order of every splitting needs to be adapted accordingly.

Construction of profiles: existence and behavior at infinity
In the last two sections, we have made use of various profiles, which account for the singularities appearing at each order in the construction of the expansion. All the profiles (K ± λ , L ±,σ λ , M) satisfy a problem such as (4.11) The existence for Problem (4.11) relies on the space V = Z ∈ H 1 loc (P) ; ∇Z ∈ L 2 (P) and Z (1 + |X|) log(2 + |X|) ∈ L 2 (P) , and the bilinear form which is coercive on the quotient space V/R, see [AGG97,Néd01]. Within this variational framework, we can build the profiles and make precise their behavior at infinity:  ψ]s ± λ and g = 0 (which is needed to define the profile K ± λ ). We need to check that P [∆; ψ]s ± λ = 0.
To do so, we introduce the bounded domain P R = {X ∈ P ; |X| < R} (with R large enough to ensure ψ(X) = 1 for |X| R), and make use of an integration by parts where B R = {X ∈ ∂P ; |X| = R}, since ∂ ν s ± λ = 0 on ∂P, and ψ is identically equal to 1 in B R . A direct integration based on the definition of s ± λ implies that I R = 0.
• Relation (4.13) shows in particular that Z → 0 at infinity, which is crucial to deriving the asymptotic expansion, see Section 4.2. Note that Assumption (H4) allows us to avoid the use of the logarithmic singularity, which naturally appears in the two-dimensional Laplace-Neumann problem, see [CCDV06, Theorem 3.25].

Representation formula
From the asymptotic expansion, one can easily infer a representation formula far from the singularity. For any y far from the inclusions, let define the Neumann function N y as the solution to  and by definition ofǔ ext 1 and thanks to (4.2), one infers the following representation formula: u δ (y) − u 0 (y) = δ ∂Π 0 ∂ τ u 1,ext flat,0 ∂ τ N y dx + O(δ). Such a representation formula is useful to characterise (or to detect) the detachment of the two inclusions tied up by the corners.

Other possible geometrical frameworks
Let us now discuss some possible extensions. We have detailed the two asymptotic models (translation and contraction) in the special situation of assumptions (H1)-(H2)-(H3). Other geometric situations may be considered, for which the two above methods have to be slightly modified. We give in Figure 5.1 two examples for which assumption (H3) is violated. In Figure 5.1a, the domains Ω L δ and Ω R δ do not lie in separated half planes orthogonal to the horizontal direction. The contraction method can be applied directly but the translation method needs to be adapted. Indeed, a vertical strip cannot be inserted between the two inclusions and we need to use two slanted half-strips. The analysis remains very similar.
In Figure 5.1b, the domains Ω L and Ω R δ lie on the upper half plane. The translation method applies directly but the contraction needs adaptation in this case. Actually, the contraction method cannot be applied if the direction d remains to be horizontal. It is necessary to assume that the inclusions converge to each other at a point which lies under the origin 0.
In Figure 5.2, we present the situation of two nearly touching squares where both methods need more important adaptations. The limit situation shows a contact on a segment and not at a single point. Other profiles must be introduced to account for the geometry near its extremity points. Let us emphasize that hypothesis (H2) is crucial to our two asymptotic expansions. In particular regarding the case of two kissing balls, the limit problem is no more a corner problem, but a problem with a cusp. The structure of the singularities is then rather different in this case as presented in see [Dau96]. Our methods cannot be applied.