Analysis of a Mogi-type model describing surface deformations induced by a magma chamber embedded in an elastic half-space

Motivated by a vulcanological problem, we establish a sound mathematical approach for surface deformation effects generated by a magma chamber embedded into Earth's interior and exerting on it a uniform hydrostatic pressure. Modeling assumptions translate the problem into classical elasto-static system (homogeneous and isotropic) in an half-space with an embedded cavity. The boundary conditions are traction-free for the air/crust boundary and uniformly hydrostatic for the crust/chamber boundary. These are complemented with zero-displacement condition at infinity (with decay rate). After a short presentation of the model and of its geophysical interest, we establish the well-posedness of the problem and provide an appropriate integral formulation for its solution for cavity with general shape. Based on that, assuming that the chamber is centred at some fixed point $\bm{z}$ and has diameter $\varepsilon>0$, small with respect to the depth $d$, we derive rigorously the principal term in the asymptotic expansion for the surface deformation as $\varepsilon=r/d\to 0^+$. Such formula provides a rigorous proof of the Mogi point source model in the case of spherical cavities generalizing it to the case of cavities of arbitrary shape.


Introduction
Measurements of crustal deformations are crucial in studying and monitoring volcanoes activity. In this context, one of the main goals is to investigate the mechanics of volcanic systems in order to obtain models for which the predicted synthetic data best fit the observed ones. Indeed, if the model is accurate it gives good images of the Earth's interior. In particular, concerning volcanoes monitoring, synthetic data might be used to localize underground magma chambers and their volume changes. In this sense, a well-established and widely used model, in the geological literature, is the one proposed by Mogi in [33] where the magma chamber is modeled by a pressurized spherical cavity of radius r, buried in a homogeneous, isotropic, elastic half-space (R 3 − ) at depth d r. Specifically, Mogi provides an explicit formula of the leading term of the displacement field on the boundary (R 2 ) of the half-space when the ratio r/d goes to zero. This model has been further investigated by McTigue in [30] who extends the analysis by detecting the second order term in the asymptotic formula. For other geometries, we mention [16] where Davies considers the case of a pressurized cavity in the shape of prolate and oblate ellipsoids, restricted to radius ratio to depth sufficiently small, deriving, also in this case, an explicit formula of the leading order term.
We stress that the analysis carried out in [16,30,33] and the derivation of the formulas therein are informal. For this reason the principal aim of this article is to contribute with a rigorous mathematical analysis of such modeling with the crucial difference relative to the shape of the cavity; our analysis does not require this shape to be a sphere or an ellipsoid but a general subdomain of the half-space. More precisely, we derive a rigorous asymptotic expansion of the boundary displacement vector field due to the presence of a cavity C, buried in an isotropic, homogeneous and elastic half-space (R 3 − ), of the form where z is the center of the cavity, ε is a suitable scaling parameter (see Section 2) and Ω, shape of the cavity, is a bounded Lipschitz domain containing the origin. Furthermore, we follow the assumption in [16,33] supposing that the boundary of the cavity C is subjected to a constant pressure p. In details, denoting by u ε the displacement vector field, we have, for y ∈ R 2 , (1.1) u k ε (y) = ε 3 |Ω| p ∇ z N (k) (z, y) : MI + O(ε 4 ), for k = 1, 2, 3, as ε → 0 (see Theorem 4.1), where u k ε stands for the k-th component of the displacement vector. Here p ε 3 represents the total work exerted by the cavity on the half-space, ∇N (k) denotes the symmetric part of the gradient related to the k-th column vector of the Neumann function that is the solution to div(C ∇N(·, y)) = δ y I in R 3 − , ∂N(·, y)/∂ν = 0 on R 2 , with the decay conditions at infinity where ∂/∂ν is the conormal derivative operator. In addition, I is the identity matrix in R 3 and M is the fourth-order moment elastic tensor defined by for q, r = 1, 2, 3, where I is the symmetric identity tensor, C is the isotropic elasticity tensor and n is the outward unit normal vector to ∂Ω. Finally, the functions θ qr , with q, r = 1, 2, 3, are the solutions to Cn on ∂Ω, with the decay conditions at infinity The leading term in (1.1) contains the information on the location of the cavity, z, on its shape and on the work p ε 3 , hence it might be used to detect the unknown cavity if boundary measurements of the displacement field are given. In order to establish these results we use the approach introduced by Ammari and Kang, based on layer potentials techniques for bounded domains, see for example [4,5,6], and following the path outlined in [8]. We highlight that proving the existence and uniqueness results for unbounded domains with unbounded boundary is, in general, much more difficult with respect to the case of bounded or exterior domains. The main obstacle is the control of both solution decay and integrability on the boundary. Indeed, the usual approach is based on the employment of some weighted functional spaces, see for example [7]. Here, taking advantage of the explicit formula of the Neumann function N and rewriting the problem into an integral formulation via layer potential techniques, we are able to prove the existence and uniqueness of the displacement field generated by an arbitrary finite cavity C in a standard Sobolev setting. These results follow showing the invertibility of the integral operators, some of them with singular kernel, that come from the layer potential technique. After that, we derive the asymptotic expansion (1.1), as ε goes to zero, for a cavity of a generic shape.
The approach based on the determination of asymptotic expansions in the case of small inclusions or cavities in bounded domains, which goes back to Friedman and Vogelius [22], has been extensively and successfully used for the reconstruction, from boundary measurements, of location and geometrical features of unknown conductivity inhomogeneities in the framework of electrical impedance tomography (see, for example, [11] and [5] for an extended bibliography). In the last decade some of these results have been extended to elasticity in the case of bounded domains, see [4]. Specifically, starting from boundary measurements given by the couple potentials/currents or deformations/tractions, in the case, respectively, of electrical impedance tomography and linear elasticity, information about the conductivity profile and the elastic parameters of the medium have been inferred. It is well known that without any a priori assumptions on the features of the problem, the reconstruction procedures give poor quality results. This is due to the severe ill-posedness of the inverse boundary value problem modeling both the electrical impedance tomography [2] and the elasticity problems [3,34]. However, in certain situations one has some a priori information about the structure of the medium to be reconstructed. These additional details allow to restore the well-posedness of the problem and, in particular, to gain uniqueness and Lipschitz continuous dependence of inclusions or cavities from the boundary measurements. One way to proceed, for instance, is to consider the medium with a smooth background conductivity or elastic parameters except for a finite number of small inhomogeneities [4,5]. Therefore, by means of partial or complete asymptotic formulas of solutions to the conductivity/elastic problems and some efficient algorithms, information about the location and size of the inclusions can be reconstructed, see for example [4,5,24].
Despite the extensive literature in this field, we remark that the mathematical problem of this work represents an intriguing novelty because we have to deal with a pressurized cavity, that is, a hole with nonzero tractions on its boundary, buried in a half-space. These two peculiarities do not allow to reduce the boundary value problem to a classical one based on cavities (see, for example, problems in [4,5] and reference therein).
The paper is organized as follows. In Section 2, we illustrate the model of ground deformations and motivate the assumptions at the basis of the simplified linearized version. In Section 3, we recall some arguments about linear elasticity and layer potentials techniques and then we analyze the well-posedness of the direct problem via an integral representation formula for the displacement field. Section 4 is devoted to the proof of the main result regarding the asymptotic formula for the boundary displacement field. In addition, as a consequence of the asymptotic expansion, we obtain the classical Mogi's formula for spherical cavities.
Notation. -We denote scalar quantities in italic type, e.g. λ, µ, ν, points and vectors in bold italic type, e.g. x, y, z and u, v, w, matrices and second-order tensors in bold type, e.g. A, B, C, and fourth-order tensors in blackboard bold type, e.g. A, B, C. The The unit outer normal vector to a surface is represented by n. The transpose of a second-order tensor A is denoted by A T and its symmetric part by A = 1 2 A + A T . To indicate the inner product between two vectors u and v we use u · v = i u i v i whereas for second-order tensors A : B = i,j a ij b ij . The cross product of two vectors u and v is denoted by u×v. With |A| we mean the norm induced by the inner product between second-order tensors, that is, |A| = √ A : A.
Acknowledgements. -The authors are indebted with Maurizio Battaglia for drawing the attention to the problem and for some useful discussion on its volcanological significance. The authors thank the New York University in Abu Dhabi (EAU) for its kind hospitality that permitted a further development of the present research.

Description of the mathematical model
Monitoring of volcanoes activity, targeted to the forecasting of volcanic hazards and development of appropriate prevention strategies, is usually performed by combining different types of geophysical measurements. Ground deformations are among the most significant being directly available (in particular, with the development of Global Positioning Systems) and, at the same time, straightly interpretable in term of elastic behaviors of the Earth's crust adjacent to the magma chamber (see [10,18,37]). A well-established model has been proposed by Mogi, [33], following previous results (see description in [17,28,36]), based on the assumption that ground deformation effects are primarily generated by the presence of an underground magma chamber exerting a uniform pressure on the surrounding medium. Precisely, the model relies on three key founding schematizations.
1. Geometry of the model. -The earth's crust is an infinite half-space (with free air/crust surface located on the plane x 3 = 0) and the magma chamber, buried in the half-space, is assumed to be spherical with radius r and depth d such that r d.

Geophysics of the crust. -
The crust is a perfectly elastic body, isotropic and homogeneous, whose deformations are described by the linearized elastostatic equations, hence are completely characterized by the Lamé parameters µ, λ (or, equivalently, Poisson ratio ν and shear modulus µ). The free air/crust boundary is assumed to be a traction-free surface.
3. Crust-chamber interaction. -The cavity describing the magma chamber is assumed to be filled with an ideal incompressible fluid at equilibrium, so that the pressure p exerted on its boundary on the external elastic medium is hydrostatic and uniform.
Assuming that the center of the sphere is located at z = (z 1 , z 2 , z 3 ) with z 3 < 0, the displacement u = (u 1 , u 2 , u 3 ) at a surface point y = (y 1 , y 2 , 0) is given by in the limit ε := r/|z 3 | → 0 (see Fig. 2.1). A higher-order approximation has been proposed by McTigue [30] with the intent of providing a formal expansion able to cover the case of a spherical body with finite (but small) positive radius.
Being based on the assumption that the ratio radius/depth ε := r/|z 3 | is small, the Mogi model corresponds to the assumption that the magma chamber is wellapproximated by a single point producing a uniform pressure in the radial direction; as such, it is sometimes referred to as a point source model. However, even if the source is reduced to a single point, the model still records the spherical form of the cavity. Different geometrical form may lead to different deformation effects.
The Mogi model has been widely applied to real data of different volcanoes to infer approximate location and strength of the magma chamber. The main benefit of such strategy lies in the fact that it provides a simple formula with an explicit appearance of the basic physical parameters depth and total work (combining pressure and volume) and, thus, that it can be readily compared with real deformation data to provide explicit forecasts.
The simplicity of formulas (2.1) makes the application model viable, but it compensates only partially the intrinsic reductions of the approach. As a consequence, variations on the geometry of the magma chamber have been proposed to provide more realistic frameworks but with the target to furnish explicit formulas for ground deformations; in this way synthetic data can be compared with real data via appropriate inversion algorithms. Such attempt has focused on oblate and prolate ellipsoids [16,40], rectangular dislocations [35] and horizontal penny-shaped cracks [20]. For completeness, we mention also the attempts to study the heterogeneity of the crust, see [29] and references therein, the case of a non-flat crust surface [39,14] and the combination of elastic properties with gravitational effects and time-dependent processes modeling of the crust, see [9,13].
In all cases, refined descriptions have the inherent drawback of requiring a detailed knowledge of the crust elastic properties. In absence of reliable complete data and measurements, the risk of introducing an additional degree of freedom in the parameter choice is substantial. This observation partly supports the approach of the Mogi model which consists in keeping as far as possible the parameters choice limited and, consequently, the model simple.
Now, let us introduce in details the boundary value problem which emerges from the previous assumptions on the geometry of the model, geophysics of the crust and crust-chamber interaction. Denoting by R 3 − the (open) half-space described by the condition x 3 < 0, the domain D occupied by the Earth's crust is D : − , describing the magma chamber, is assumed to be an open set with bounded boundary ∂C. Hence, the boundary of D is composed by two disconnected components: the two-dimensional plane which constitutes the free air/crust border, and the set ∂C, corresponding to the crust/chamber edge. In this setting, we end up with the following boundary value problem which we consider with the asymptotic conditions at infinity The conjunction of conditions on ∂D and at infinity determines a displacementtraction problem. At this point, the model provides the displacement u of a generic finite cavity C. The next step is to deduce a corresponding point source model, in the spirit of the Mogi spherical one. To this aim, we assume the cavity C of the form where d, r > 0 are characteristic length-scales for depth and diameter of the cavity, its center dz belongs to R 3 − and its shape Ω is a bounded domain containing the origin. The Mogi model corresponds to Ω given by a sphere of radius 1.
Introducing the rescaling (x, u) → (x/d, u/r) and denoting the new variables again by x and u, the above problem takes the form (with unchanged far-field conditions) where ε = r/d, C ε := z + εΩ, D ε := R 3 − C ε and p is a "rescaled" pressure, ratio of the original pressure p and ε.
In the next section we study the well-posedness of the linear elastic model in the case of a finite pressurized cavity C and then we derive the asymptotic expansion of the solution when C := C ε := z + εΩ.

Integral representation and well-posedness of the direct problem
Since C := λI ⊗ I + 2µI, the elastostatic Lamé operator L for a homogeneous and isotropic elastic medium is given by where u represents the vector of the displacements and ∇u = 1 2 ∇u + ∇u T the strain tensor. With ∂u/∂ν we depict the conormal derivative on the boundary of a domain, that is, the traction vector, which has the expression ∂u ∂ν := (C ∇u)n = λ(div u)n + 2µ( ∇u)n or, equivalently, ∂u ∂ν = 2µ ∂u ∂n + λ(div u)n + µ(n × rot u).
Here, we analyze the linear elastostatic boundary value problem where C is the cavity and p is a constant representing the pressure. For the Lamé parameters, we consider the physical range 3λ + 2µ > 0 and µ > 0 which ensures positive definiteness of C. For the sequel, we recall that the positive definiteness of the tensor C implies the strong ellipticity which corresponds to the request µ > 0 and λ + 2µ > 0, see [23]. The aim of this section is to provide an integral representation formula and to establish the well-posedness of the problem. To do that, we consider three steps: -firstly, we recall Betti's formulas, definition and some properties of single and double layer potentials of linear elasticity; -then, we give the expression of the fundamental solution N of the half-space with null traction on the boundary, found by Mindlin in [31,32]; -finally, we represent the solution to (3.1) by an integral formula through the fundamental solution of the half-space.
All these objects will be used to prove the well-posedness of the problem (3.1).

3.1.
Preliminaries. -We recall Betti's formulas for the Lamé system which can be obtained by integration by parts, see for example [4,27]. Given a bounded Lipschitz domain C ⊂ R 3 and two vectors u,v ∈ R 3 , the first Betti formula is where the quadratic form Q associated to the Lamé system is From (3.2) it is straightforward to find the second Betti formula Formula (3.2) will be used to prove that the solution of (3.1) is unique, and the equality where δ 0 is the Dirac function centred at 0. Setting C µ,ν : where ν is the Poisson ratio which is related to λ and µ by the identity where δ ij is the Kronecker symbol and Γ ij stands for the i-th component of the displacement when a force is applied in the j-th direction at the point 0. For the reader's convenience, we write also the gradient of Γ to highlight its behaviour at infinity Therefore, from (3.4) and (3.5) it is straightforward to see that With the Kelvin matrix Γ at hand, we recall the definition of the single and double layer potentials corresponding to the operator L. Given ϕ ∈ L 2 (∂C), we set (see [4,5,27]) where ∂Γ/∂ν denotes the conormal derivative applied to each column of the matrix Γ.
In the sequel the subscripts + and − indicate the limits from outside and inside C, respectively. The double layer potential and the conormal derivative of the single layer potential satisfy the jump relations where K and K * are the L 2 -adjoint Neumann-Poincaré boundary integral operators defined, in the sense of Cauchy principal value, by It is worth noticing that these two operators are not compact even on smooth domains, in contrast with the analogous operators for the Laplace equation (see [5]), due to the presence in their kernels of the terms which make the kernel not integrable. Indeed, even in the case of smooth domains, we cannot approximate locally the terms n × (x − y) with a smooth function, that is, by power of |x − y| via Taylor expansion, in order to obtain an integrable kernel on ∂C. Therefore, the analysis to prove the invertibility of the operators in (3.8) is intricate and usually based on regularizing operators (see [27]) in the case of smooth domains. For the Lipschitz domains the analysis is much more involved and based on Rellich formulas (see [15] and its companion article [19]).

3.2.
Fundamental solution of the half-space.
-In this subsection we show the explicit expression of the Neumann function of the half-space presented for the first time in [31] by means of Galerkin vector and nuclei of strain of the theory of linear elasticity, and secondly in [32] using the Papkovich-Neuber representation of the displacement vector and the potential theory. We consider the boundary value problem The Neumann function of (3.9) is the kernel N of the integral operator giving the solution to the problem. Given y = (y 1 , y 2 , y 3 ), we set y := (y 1 , y 2 , −y 3 ).
Theorem 3.1. -The Neumann function N of problem (3.9) can be decomposed as For the proof of Theorem 3.1 see the appendix. Uniqueness of the solution to (3.9) is similar to the one for problem (3.1) which we present in the following section.
The matrix R, defined by gives the regular part of the Neumann function since the singular point η = 0 corresponds to y = (y 1 , y 2 , −y 3 ) with y 3 < 0, which belongs to R 3 + . To convert the problem (3.1) into an integral form, bounds on the decay at infinity of the Neumann function and its derivative at infinity are needed. for any x, y ∈ R 3 − with |x| M x and |y| M y .
The proof of this proposition is a consequence of the properties on homogeneous functions and their derivatives, noticing that R k are homogeneous of degree −k, for k = 1, 2, 3.

3.3.
Representation formula. -Next, we derive an integral representation formula for u solution to problem (3.1). For, we make use of single and double layer potentials defined in (3.7) and integral contributions relative to the regular part R of the Neumann function N, defined by (3.13) where ϕ ∈ L 2 (∂C). (3.13), pn is the boundary condition in (3.1) and f is the trace of u on ∂C.
Before proving this theorem, we observe that if a solution to (3.1) exists, then f solves the integral equation obtained by the application of the trace properties of the double layer potential (3.8) in formula (3.14).
with r sufficiently large such that to contain the cavity C; additionally, we define ∂B h r (0) as the intersection of the hemisphere with the boundary of the half-space, and with ∂B b r (0) the spherical cap (see Figure 3.1). Now, we apply Betti's formula (3.3) to u and the k-th column vector of N, indicated by N (k) , for k = 1, 2, 3, in Ω r,ε , hence since, from (3.1) and the boundary condition in (3.9), We show that the term I 1 goes to zero by using the behaviour at infinity of u given in (3.1) and of the Neumann function given in (3.12). Indeed, we have This last integral can be estimated by means of the spherical coordinates Again, passing through spherical coordinates, we get as r → +∞, since |∇u| = o(r −1 ). The integral I 2 gives the value of the function u in y as ε goes to zero. Indeed, we have ∂Bε(y) since the second integral has a continuous kernel. On the other hand − ∂Bε(y) The latter integral tends to zero when ε goes to zero because and this last integral is bounded when ε goes to zero. Let finally observe that −u(y)· ∂Bε(y) where the latter integral tends to zero as ε → 0, since R (k) represents the regular part of the Neumann function. To deal with the first integral, we preliminarily observe that direct differentiation gives where c ν := (1−2ν)/(8π (1−ν)). We substitute this expression into the integral (3.17) and we take into account that To compute this last integral we use spherical coordinates, that is, where ϕ ∈ [0, π] and θ ∈ [0, 2π). From a simple calculation it follows for any h and k. Hence, (3.17) becomes − u(y) · ∂Bε(y) Employing again the spherical coordinates and the definition of c ν , we find that Putting together all the results in (3.21) and (3.22), we find that Using the definition of single and double layer potentials (3.7), (3.13) and splitting N as Γ + R, Formula (3.14) holds.
From the behaviour of the Neumann function given in (3.12) and the representation formula in (3.14), we immediately get . In particular, in order to prove the injectivity of the operator (3.24) we show the uniqueness of u following the classical approach based on the application of the Betti's formula (3.2) and the energy method, see [21,27]. From the injectivity, it follows the existence of u proving the surjectivity of (3.24) which is obtained by the application of the index theory regarding bounded and linear operators.
First of all, let us recall the closed range theorem due to Banach (see [25,41]). Proof. -The assertion of this lemma is based on the invertibility of the operator 1 2 I + K * studied in [15]; it is known that is a bounded linear operator, injective and with dense and closed range. Therefore, from Theorem 3.5 we have and Im(1/2I + K) is closed. Then, it follows that the operator 1 2 I + K : L 2 (∂C) → L 2 (∂C) is bijective and the assertion follows exploiting the bounded inverse theorem.
Thus D R (S) ⊂ C(∂C), where C(∂C) indicates the space of continuous function on ∂C. Since each component of the matrix H is uniformly continuous on the compact set ∂C × ∂C, for every ε > 0 there exists δ > 0 such that for all x, y, z ∈ ∂C with |x − y| < δ. Since for all x, y ∈ ∂C and ϕ ∈ S, that is, D R (S) is equicontinuous. The assertion follows from Ascoli-Arzelà Theorem and noticing that C(∂C) is dense in L 2 (∂C).
We now prove Proof. -Let u 1 and u 2 be solutions to (3.1). Then the difference v := u 1 − u 2 solves the homogeneous version of (3.1), that is, where we make use of the decay condition at infinity comes from Corollary 3.4. Applying Betti's formula (3.2) to v in Ω r = (R 3 − ∩ B r (0)) C, we find where Q is the quadratic form Q(v, v) = λ(div v) 2 + 2µ| ∇v| 2 . From the behaviour of v and the boundary conditions (3.26), we estimate the previous integral defined on the surface of Ω r ; contributions over the surface of the cavity and the intersection of the hemisphere with the half-space are null by means of (3.26), whereas on the spherical cap As already done in (3.16) to obtain the representation formula, this integral can be evaluated by spherical coordinates; in particular, it tends to zero when r → +∞. Therefore Since the quadratic form is positive definite for the parameters range 3λ + 2µ > 0 and µ > 0, we have that It follows that the rigid displacements v = a + Ax, with a ∈ R 3 and A belonging to the space of the anti-symmetric matrices (see [4,12]), could be the only nonzero solutions which satisfy (3.25), the boundary conditions in (3.26) and (3.27). However, in this case they are excluded thanks to the behaviour of the function v at infinity. Hence, we obtain that v = 0, that is, The uniqueness result for problem (3.1) ensures the injectivity of the operator (3.24). In order to prove the surjectivity of the operator (3.24) we recall, for the reader's convenience, the definition of the index of an operator (see [1,25]).
Definition 3.1. -Given a bounded operator T : X → Y between two Banach spaces, the index of the operator T is the extended real number defined as where dim(Ker(T )) is called the nullity and dim(Y / Im(T )) the defect of T . In particular, when the nullity and the defect are both finite the operator T is said to be Fredholm.
We recall also an important theorem regarding the index of a bounded linear operator perturbed with a compact operator (see [1]). Proof. -Uniqueness follows from Theorem 3.8 and the existence from Theorem 3.10.

Rigorous derivation of the asymptotic expansion
In this section, with the integral representation formula (3.14) at hand, we consider the hypothesis that the cavity C is small compared to the distance from the boundary of the half-space. The aim is to derive an asymptotic expansion of the solution u.
In particular, let us take the cavity, that from now on we denote by C ε to highlight the dependence on ε, as where Ω is a bounded Lipschitz domain containing the origin. At the same time, we write the solution of the boundary value problem (3.1) as u ε . From (3.14), recalling that N = Γ + R, we have for k = 1, 2, 3, where u k ε indicates the k-th component of the displacement vector and f = f ε is the solution of (3.15), that is, where we add the dependence on ε to all the layer potentials to distinguish them, in the sequel, from the layer potential defined over a domain independent of ε. In what follows, with I we indicate the fourth-order symmetric tensor such that IA = A and for any fixed value of ε > 0, given h : ∂C ε → R 3 , we introduce the function h : ∂Ω → R 3 defined by h(ζ) := h(z + εζ), ζ ∈ ∂Ω.
Moreover, we consider the functions θ qr , for q, r = 1, 2, 3, solutions to Cn on ∂Ω, with the decay conditions at infinity where the condition ∂θ qr /∂ν has to be read as We now state our main result.
Theorem 4.1 (asymptotic expansion). -There exists ε 0 > 0 such that for all ε ∈ (0, ε 0 ) the following expansion holds at any y ∈ R 2 : Before proving the theorem on the asymptotic expansion of u ε , we need to present some results. , when x = z + εζ, with ζ ∈ ∂Ω, is such that is uniformly bounded in ε. Moreover, when ε is sufficiently small, we have Proof. -At the point z + εζ, where ζ ∈ ∂Ω, we obtain Therefore, recalling that the kernel ∂R/∂ν(η) is continuous, we get where Λ Ω,ε C , with C independent of ε. For the integral we use the explicit expression of the conormal derivative of the fundamental solution of the Lamé operator given in (3.18). In particular, since (3.18) is a homogeneous function of degree −2, with the substitution t = z + εη, we find for h, k = 1, 2, 3. Therefore, we immediately obtain that Evaluating the other integrals in (4.2) we obtain hence, choosing t = z + εη, with η ∈ ∂Ω, we find where the last equality follows noticing that the fundamental solution is homogeneous of degree −1. In a similar way hence, taking again t = z + εη, we find and since R is regular it follows that . Relation (4.7) follows putting together the result in (4.8), (4.9), (4.10) and (4.11).
For ease of reading, we define the function w : ∂Ω → ∂Ω as Taking the problem with decay conditions at infinity we show that w(x), for x ∈ ∂Ω, is the trace of v on the boundary of Ω. The wellposedness of this problem is a classical result in the theory of linear elasticity so we remind the reader, for example, to [21,23,27].
where v is the solution to (4.13) and (4.14).
Proof. -Applying the second Betti formula to the fundamental solution Γ and the function v into the domain B r (0) (Ω ∪ B ε (x)), with ε > 0 and r > 0 sufficiently large so that the domain contains the cavity Ω, we obtain, as done in a similar way in the proof of the Theorem 3.3, Therefore, from the single and double layer potential properties for the elastostatic equations, we find that is, the assertion.
We note that the function v, as well as its trace w on ∂Ω, can be written in terms of the functions θ qr . Indeed, taking v = θ qr δ qr , q, r = 1, 2, 3, defined in (4.1). Since y ∈ R 2 and x ∈ ∂C ε = z + εζ, with ζ ∈ ∂Ω, we consider the Taylor expansion for the Neumann function, that is, for k = 1, 2, 3. By the change of variable x = z + εζ and substituting (4.15) in I (k) 1 , we find The integral I For the integral I (k) 12 , we use the equality n · ∇N (k) ζ = ∇N (k) : (n(ζ) ⊗ ζ), therefore For the term I (k) 2 we use the result in Lemma 4.2 and the Taylor expansion of the conormal derivative of N (k) (x, y), for k = 1, 2, 3. In particular, for x = z + εζ, when ζ ∈ ∂Ω and y ∈ R 2 , we consider only the first term of the asymptotic expansion, that is, for any k, where w is defined in (4.12). Since ∂N (k) /∂ν(ζ) = C ∇N (k) n(ζ), we have Collecting the result in (4.16) and (4.17), Equation (4.1) becomes Now, handling this expression, we highlight the moment elastic tensor. We have indeed, for any i, j = 1, 2, 3, it follows where e j is the j-th unit vector of R 3 . Hence, by (4.18) and taking the symmetric part of ∇N (k) , for any k, we find Using the symmetries of C, we have for k = 1, 2, 3. Now, taking into account that I = II and using the equality w = θ qr δ qr , with q, r = 1, 2, 3, we have the assertion.
The Mogi model. -In this subsection, starting from the asymptotic expansion (4.5), that is, u k ε (y) = ε 3 |Ω| p ∇ z N (k) (z, y) : MI + O(ε 4 ), k = 1, 2, 3, where M is the tensor given in (4.6), we recover the Mogi model, presented within the Section 2, related to a spherical cavity. We first recall that where, in the last equality, we use the link between the functions w and θ qr , that is, w = θ qr δ qr , q, r = 1, 2, 3. Therefore, to get the Mogi's formula, we first find the explicit expression of w when the cavity Ω is the unit sphere and then we calculate the gradient of the Neumann function N.
We recall that w is the trace on the boundary of the cavity of the solution to the external problem We look for a solution of the form v(x) = φ(r) x with r := |x|, so that By direct substitution, since n = x on ∂B, we get Thus, we need to find a function φ : The condition at infinity implies that B = 0 and A = 1/4µ. Therefore, the solution is v(x) = x/4µ|x| 3 , which implies that With the function w at hand, we have that Through the use of spherical coordinates and orthogonality relations for the circular functions, it holds hence the second-order tensor MI is given by For the Neumann's function N (see the appendix for its explicit expression), we are interested only in the trace of ∇ z N(z, y) computed at y 3 = 0. Evaluating N = N(z, y) at y 3 = 0, we get 3 f 3 where α, β = 1, 2 and κ µ = 1/(4πµ), with f = 1/|z − y| and g = 1/ |z − y| − z 3 . Let ρ 2 := (z 1 − y 1 ) 2 + (z 2 − y 2 ) 2 . Using the identities we deduce the following formulas for some of the derivatives of κ −1 As a consequence, we obtain (4.20) Tr Combining (4.19), (4.20) and using the explicit expression for f , we find that are the components given in (2.1).

Appendix. Neumann function for the half-space with zero traction
For the reader's convenience, we provide here the complete derivation of the explicit formula for the Neumann function N of the problem as stated in Theorem 3.1. By definition, N is the kernel of the integral operator which associates to the forcing term b the solution v to the boundary value problem, viz. v The explicit expression for the Neumann function has been determined by Mindlin in [31,32] using different approaches. Here, we follow the second one which is based on the Papkovich-Neuber representation of the displacement v and potential theory.
There exist other mathematical approaches to obtain the Neumann function. See for example [38] for an application of a Fourier method. Let us introduce the functions observing that, apart from ∂ i φ 0 = −x i φ 3 0 , i = 1, 2, 3, the following identities hold true for α = 1, 2, We denote by φ and φ the values φ 0 (x + e 3 ) and φ 0 (x − e 3 ), respectively, with analogous notation for ψ 0 , that is, ψ and ψ are the values ψ 0 (x + e 3 ) and ψ 0 (x − e 3 ) respectively.
To establish (A.1), we observe that the columns N (i) of N are determined by solving the equation Lv = e i δ for i = 1, 2, 3 and using the Papkovich-Neuber representation where δ ij is the Kronecker symbol. The coupling between h and β is determined by the boundary conditions on the plane {x 3 = 0}, which are Given y := (y 1 , y 2 , y 3 ), let y be the reflexed point. Set Denoting by f, g the action of the distribution f on the function g, we determine h and β taking advantage of the relation (which is deduced from the second Green identity) Proof of Proposition A.1.. -To determine N , we consider separately the case of horizontal and vertical forcing. By symmetry, x 1 and x 2 can be interchanged.
Horizontal force: Lv = e 1 δ. -We choose h 2 = 0, so that boundary conditions become Differentiating the first equation with respect to x 1 , the second with respect to x 2 and taking the difference, we obtain , which suggests, after integration with respect to x 2 , the choice F := ∂ 3 h 1 . Applying (A.4), and thus h 1 = −(φ + φ).
The fundamental solution N = N(x, y) in the half-space {x 3 < 0} is such that its columns v 1 , v 2 and v 3 solve Lv i = δ y e i where δ y is the Dirac delta concentrated at y = (y 1 , y 2 , y 3 ) with y 3 < 0. Thus, the Neumann function N is given by N(x, y) = 1 as a result of the homogeneity of δ and the second order degree of L.
Recalling the definitions of φ, φ, ψ and computing at (x 1 − y 1 , x 2 − y 2 , x 3 )/|y 3 |, we obtain the identities where y = (y 1 , y 2 , −y 3 ). Hence, the components of C −1 µ,ν N are given by Recollecting the expression for fundamental solution Γ in the whole space and using the relation f = g − (x 3 + y 3 ) f g, the above formulas can be rewritten as N = Γ + R, where Γ is computed at x − y and the component R ij of R are given by where η α = x α − y α for α = 1, 2 and η 3 = x 3 + y 3 , which can be recombined as for i, j = 1, 2, 3. Since x 3 = η 3 −y 3 , we obtain the decomposition R ij := R 1 ij +R 2 ij +R 3 ij , where The proof of Theorem 3.1 is complete.