Zohar Yosibash
email: zohary@bgu.ac.il

Samuel Shannon

Monique Dauge
email: monique.dauge@univ-rennes1.fr

Martin Costabel
email: martin.costabel@univ-rennes1.fr

Circular edge

INTRODUCTION.

Solutions of elliptic boundary value problems over two dimensional domains, for example those arising in heat transfer and elasticity, when posed and solved in non-smooth domains like polygons, have non-smooth parts. These are described in terms of special singular functions depending on the geometry and the differential operators on one hand, and of unknown coefficients depending on the given right hand side and boundary conditions on the other hand, see e.g. [START_REF] Williams | Stress singularities resulting from various boundary conditions in angular corners of plates in extension[END_REF].

In three dimensional domains as polyhedra both vertex and edge singularities exist, see [START_REF] Dauge | Elliptic boundary value problems in corner domains -smoothness and asymptotics of solutions[END_REF][START_REF] Von Petersdorff | Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron[END_REF]. For straight edges we have provided explicit representation of the singular solutions [START_REF] Costabel | A quasidual function method for extracting edge stress intensity functions[END_REF][START_REF] Omer | Edge flux intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF][START_REF] Yosibash | Edge stress intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF] as a series characterized by:

• an exponent α which belongs to a discrete set {α k , k ∈ N} of eigen-values depending only on the geometry and the operator, and which determines the level of non-smoothness of the singularity. Any eigen-value α k is computed by solving a 2-D problem. • eigen-functions φ k,0 (ϕ) which depends on the geometry of the domain and the operator. These eigen-functions are computed by solving a set of 2-D problems. • a function along the edge, denoted by A k (s) ( s is a coordinate along the edge) and called "Edge Flux/Stress Intensity Function" (EFIF/ESIF) which determines the "amount of energy" residing in each singularity.

Here we concentrate on circular edges (a "penny-shaped crack" being a special renown case) in 3-D domain, and derive explicitly singular series expansion in the vicinity of such an edge first for the simplest scalar elliptic operator, the Laplace operator, and then for the elasticity system. We demonstrate that our asymptotic solution for the elasticity system in the simplified case of the penny-shaped crack under axi-symmetric boundary conditions and geometry reduces to the one presented in [START_REF] Leung | Eigenfunction expansion for penny-shaped and circumferential cracks[END_REF].

From the engineering perspective the edge flux/stress intensity functions A k (s) (EFIFs/ESIFs) when α k < 1 ( α 1 = 1 2 for the penny-shaped crack) are of major importance. These are used to predict failure initiation and propagation, and are an important ingredient in any failure law for cracked and V-notched structures. To efficiently and accurately compute them, the asymptotic solution has first to be explicitly derived.

This work is motivated by the need to compute edge stress intensity functions (ESIFs) for elasticity problems in 3-D domains. These are of significant engineering importance in cracked and V-notched structures, in which the ESIFs may (and often do) vary along the crack front.

In order to explain the ideas of the implementation of the method and to test its efficiency, we consider the Laplace operator first. This is a simpler elliptic operator that allows more transparent analytic computations and invokes all necessary characteristics of the elasticity system. Thus, the characteristics of the solution can be more easily addressed.

The first three singular terms for the solution of the Laplace equation in the vicinity of a circular edge with homogeneous Dirichlet boundary conditions were analyzed from a theoretical viewpoint in [START_REF] Von Petersdorff | Singularities of the solution of the Laplacian in domains with circular edges[END_REF]. The first two terms in the Neumann case are provided in [START_REF] Bach | Propagation of a penny shaped crack under the Irwin criterion[END_REF] when the edge is the boundary of a smooth plane crack surface. For the elasticity system, Leblond&Torlai [START_REF] Leblond | The stress-field near the front of an arbitrarily shaped crack in a 3-dimensional elastic body[END_REF] provided the machinery for the pointwise derivation of the solution up to second order for a general curved crack, whereas Leung&Su derived the asymptotic series for the axi-symmetric case in [START_REF] Leung | Eigenfunction expansion for penny-shaped and circumferential cracks[END_REF]. Herein, we present a different approach enabling the computation of the entire series solution up to an arbitrary order for any circular edge be it in an axi-symmetric or non-axisymmetric setting. This explicit representation illustrates the distinct two levels of complexity of shadow terms associated with the curved singular edge.

The Laplace equation and the notation are introduced in section 2. The systematic derivation of the singular series expansion is presented for homogeneous Dirichlet and Neumann boundary conditions. Both axi-symmetric and non-axi-symmetric configurations are addressed, and for a penny-shaped crack, a circumferential crack, an external crack and a 3π/2 reentrant corner we also present closed form explicit singular series expansion. The asymptotic expansion is provided in terms of eigen-functions, their shadows, the EFIFs and their derivatives.

We then use the machinery developed in the framework of the Laplace equation to provide explicit representation of similar cases associated with the elasticity system (with clamped or traction free boundary conditions) in section 3. The elasticity explicit asymptotic solution is mandatory for the computation of edge stress intensity functions for cracks occurring usually in pipes and pressure vessels.

ASYMPTOTIC SOLUTION FOR THE LAPLACE EQUATION

As a model, we choose a domain generated by rotating the 2-D plane Ω having a reentrant corner with an opening ω ∈ (0, 2π] (the case of a crack, ω = 2π , is included) along the axis x 3 , as shown in Figure 1. The cylindrical coordinate system r, θ, x 3 and the coordinate system attached to the circular edge ρ, ϕ, θ are shown in Figure 1. It is important to emphasize that domain's geometry does not need

H N ! ! " 4 # $ % 1 % 2 % 3 & # 1 # 2 FIGURE 1.
Model domain of interest Ω and the coordinate systems. to be axi-symmetric, but only the generated circular singular edge. An example of several different circular singular edges to which the analysis in this manuscript is applicable are shown in Figure 2. For example, the lower singular edge in Figure 2 (a) is determined by ϕ ∈ (-π, π/2) , the outer circular crack in Figure 2 (b) is determined by ϕ ∈ (0, 2π) whereas the penny-shaped crack in (c) is determined by ϕ ∈ (-π, π) . Finally the re-entrant corner with the solid angle ω in Figure 2 (d) is determined by ϕ ∈ ((π -ω)/2, (π + ω)/2) .

= > ? @

FIGURE 2. Different types of singular circular edges (only a sector is plotted so the circular edge is clearly visible. The domain in (c) includes the renown "penny-shaped" crack.

For the Laplace operator, we are interested in solutions τ (x) of the equation:

3D τ def = ∂ rr + 1 r ∂ r + 1 r 2 ∂ θθ + ∂ 33 τ = 0, (1)
where

∂ r def = ∂ ∂r , ∂ rr def = ∂ 2 ∂r 2 ∂ θθ def = ∂ 2 ∂θ 2 and ∂ 33 def = ∂ 2 ∂x 2 3 . Homogeneous Dirichlet or Neumann boundary conditions are considered on Γ 1 × [0, 2π] and Γ 2 × [0, 2π] .
The solution in the vicinity of the edge is of interest so we perform a change of coordinates as follows:

r = ρ cos ϕ + R, x 3 = ρ sin ϕ. (2)
The Laplace operator in the new coordinates is given by:

3D = ∂ ρρ + 1 ρ ∂ ρ + 1 ρ 2 ∂ ϕϕ + 1 r cos ϕ∂ ρ - 1 ρ sin ϕ∂ ϕ + 1 r 2 ∂ θθ (3)
2.1. Axi-symmetric case. For an axi-symmetric domain and boundary conditions the solution is independent of θ . Then the last term in (3) vanishes and the Laplace operator for ρ/R 1 reads:

Axi = ∂ ρρ + 1 ρ ∂ ρ + 1 ρ 2 ∂ ϕϕ + 1 r cos ϕ∂ ρ - 1 ρ sin ϕ∂ ϕ . ( 4
) Remark 1. Since r → ∞ as R → ∞ , one may observe that Axi R→∞ -→ 2D .
Axisymmetric solutions τ of (1) are equivalently the solutions of r R Axi τ = 0 , i.e.

(

+ ρ R cos ϕ) ∂ ρρ + 1 ρ ∂ ρ + 1 ρ 2 ∂ ϕϕ τ + 1 R cos ϕ∂ ρ - 1 ρ sin ϕ∂ ϕ τ = 0. 1
Multiplying by ρ 2 , we find another equivalent equation

(ρ∂ ρ ) 2 + ∂ ϕϕ τ + ρ R cos ϕ(ρ∂ ρ ) -sin ϕ∂ ϕ + cos ϕ (ρ∂ ρ ) 2 + ∂ ϕϕ τ = 0. (6)
The solution in the vicinity of the singular point in the 2-D cross-section Ω can be obtained in a simple form as an asymptotic series defined by eigen-pairs of a one-dimensional boundary value problem on the interval ϕ ∈ (ϕ 1 , ϕ 1 + ω) . If we denote one such eigen-pair by α and φ 0 (ϕ) , then it is conceivable (as to be shown in the sequel) that for the axi-symmetric case a solution is formed as an asymptotic series of the form:

τ = Aρ α ∞ i=0 ρ R i φ i (ϕ) (7)
Boundary Conditions: To satisfy the homogeneous boundary conditions, the series representation has to satisfy the following constraints on ϕ = ϕ 1 and ϕ = ϕ 2 = ϕ 1 + ω :

φ i (ϕ = ϕ 1 , ϕ 2 ) = 0 in Dirichlet case (8) φ i (ϕ = ϕ 1 , ϕ 2 ) = 0 in Neumann case (9)
Substitute ( 7) in ( 6) to obtain:

A α 2 φ 0 + φ 0 (10)
+ ρ R (α + 1) 2 φ 1 + φ 1 + α cos ϕφ 0 -sin ϕφ 0 + cos ϕ α 2 φ 0 + φ 0 + ρ 2 R 2 (α + 2) 2 φ 2 + φ 2 + (α + 1) cos ϕφ 1 -sin ϕφ 1 + cos ϕ (α + 1) 2 φ 1 + φ 1 + ρ 3 R 3 (α + 3) 2 φ 3 + φ 3 + (α + 2) cos ϕφ 2 -sin ϕφ 2 + cos ϕ (α + 2) 2 φ 2 + φ 2 + • • • = 0
To satisfy the above equation for any A and ρ , the following relationships must hold:

α 2 φ 0 + φ 0 = 0 (11) (α + 1) 2 φ 1 + φ 1 = -α cos ϕφ 0 -sin ϕφ 0 -cos ϕ α 2 φ 0 + φ 0 (12) (α + 2) 2 φ 2 + φ 2 = -(α + 1) cos ϕφ 1 -sin ϕφ 1 -cos ϕ (α + 1) 2 φ 1 + φ 1 (13) (α + 3) 2 φ 3 + φ 3 = -(α + 2) cos ϕφ 2 -sin ϕφ 2 -cos ϕ (α + 2) 2 φ 2 + φ 2 (14) • • •
Substituting the RHS of equation [START_REF] Sokolnikoff | Mathematical Theory of Elasticity[END_REF] in [START_REF] Williams | Stress singularities resulting from various boundary conditions in angular corners of plates in extension[END_REF] one obtains:

α 2 φ 0 + φ 0 = 0, ϕ 1 < ϕ < ϕ 2 (15) (α + 1) 2 φ 1 + φ 1 = -α cos ϕφ 0 -sin ϕφ 0 , ϕ 1 < ϕ < ϕ 2 (16) (α + i) 2 φ i + φ i = -(α + i)(α + i -1) cos ϕφ i-1 -sin ϕφ i-1 + cos ϕφ i-1 (17) i ≥ 2, ϕ 1 < ϕ < ϕ 2
These equations have to be completed by the boundary conditions [START_REF] Omer | Edge flux intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF] or [START_REF] Von Petersdorff | Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron[END_REF]. Note the following:

• The equation (15) with BCs (8) or ( 9) is the one dimensional eigenvalue problem corresponding to the 2-D problem over Ω , with eigen-value α and eigen-function φ 0 . Traditionally φ 0 is called primal eigen-function. • A recursive system of ordinary differential equations is obtained -once φ 0 is computed from (15) it can be inserted in (16) to obtain φ 1 then these both can be inserted in (17) to obtain φ 2 , etc. • Only particular solutions in ( 16) and (17) are required. Because [START_REF] Leung | Eigenfunction expansion for penny-shaped and circumferential cracks[END_REF] corresponds only to one representative eigen-pair, the complete solution should be a sum over all eigen-pairs α k , φ k,i , thus is a double sum series:

τ = k A k ρ α k ∞ i=0 ρ R i φ k,i (ϕ) (18)
Remark 2. For each eigen-function and shadow φ k,i (ϕ) the first index k represents the eigen-value α k to which this eigen-function is associated, whereas the second index i ≥ 1 represents the rank of the shadow terms. Here α k = kπ ω , where k = 0, 1, 2, . . . for homogeneous Neumann BCs, and k = 1, 2, 3, . . . for homogeneous Dirichlet BCs. Remark 3. If α + 1 is not an eigenvalue of equation (15) with BCs (8) or [START_REF] Von Petersdorff | Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron[END_REF], there exists a unique solution Φ 1 to equation (16). On the other hand, if α + 1 is itself an eigenvalue, then it can either happen that (16) has no solution (then the ansatz [START_REF] Leung | Eigenfunction expansion for penny-shaped and circumferential cracks[END_REF] has to be completed with logarithmic terms), or (16) has infinitely many solutions. The same situation holds for equation (17), depending on whether α + i is an eigenvalue or not.

In the special case of a crack, we have α k = k 2 , therefore resonances (i.e. α + i is an eigenvalue) always occur. Nevertheless, as proved in [START_REF] Costabel | Asymptotics without logarithmic terms for crack problems[END_REF], logarithmic terms never appear: Equations ( 16) and (17) with Dirichlet or Neumann BCs are always solvable. An orthogonality condition against the eigenvector makes the solution unique, see (20).

k α k φ k,0 0 0 1 1 1 2 sin ϕ 2 2 1 cos ϕ 3 3 2 sin 3ϕ 2 4 2 cos 2ϕ TABLE 1.
First four eigen-pairs for a crack with homogeneous Neumann BCs.

Note that the angular part of the first singular function is sin ϕ 2 instead of the "usual" formula cos ϕ 2 : This is due to the choice of the angular coordinate ϕ ∈ (-π, π) instead of ϕ ∈ (0, 2π) .

Equations ( 15)-( 17) can be solved (cf. Remark 3) for α k = 0, 1/2, 1, 3/2 , obtaining φ 0,i , φ 1,i , φ 2,i , φ 3,i . They yield the following series solution for a penny-shaped crack with homogeneous Neumann BCs:

τ = A 0 (19) + A 1 ρ 1 2 sin ϕ 2 + ρ R 1 4 sin ϕ 2 + ρ R 2 1 12 sin ϕ 2 - 3 32 sin 3ϕ 2 + ρ R 3 1 16 sin ϕ 2 - 1 30 sin 3ϕ 2 + 5 128 sin 5ϕ 2 + • • • + A 2 ρ cos ϕ - ρ R 1 4 + ρ R 2 3 16 cos ϕ - ρ R 3 9 128 + 5 64 cos 2ϕ + • • • + A 3 ρ 3 2 sin 3ϕ 2 - ρ R 1 4 sin ϕ 2 - ρ R 2 1 32 3 sin ϕ 2 - 16 5 sin 3ϕ 2 + ρ R 3 - 3 40 sin ϕ 2 + 5 128 sin 3ϕ 2 - 3 70 sin 5ϕ 2 + • • • + • • •
It is worthwhile to notice that we enforced the following orthogonality conditions on the shadow terms

ϕ 2 =π ϕ 1 =-π φ k,i (ϕ) φ k+i,0 (ϕ) dϕ = 0, k = 0, 1, 2, 3 and i = 1, 2, 3. (20)
making them unique. One may notice that for R → ∞ (the crack edge curvature tends to zero) only the first terms are non-zero so the solution (19) reduces to the 2-D solution:

τ R→∞ -→ A 0 + A 1 ρ 1/2 sin ϕ 2 + A 2 ρ cos ϕ + A 3 ρ 3/2 sin 3ϕ 2 + • • • Remark 4.
The eigen-functions and shadows associated with A 0 and A 2 above are polynomials in local Cartesian variables z 1 := ρ cos ϕ and z 2 = ρ sin ϕ . This can be predicted by the general theory [START_REF] Kondratiev | Boundary value problems for elliptic equations in domains with conical or angular points[END_REF][START_REF] Dauge | Elliptic boundary value problems in corner domains -smoothness and asymptotics of solutions[END_REF].

To verify the correctness of the solution (19) we consider a torus with an inner radius r 1 = 1.5 and an outer radius r 2 = 2.5 having a circular crack with the tip at R = 2 , see Figure 3.

Taking A 1 = 1 and A k = 0 , k = 1 , (notice that the outer boundary of the torus is ρ out = 1/2 and ρ/R = 1/4 in the considered example) we prescribed on the outer surface of the torus Dirichlet with homogeneous Neumann boundary conditions on the crack face. Because the problem is axisymmetric we construct a two-dimensional axi-symmetric finite element (FE) model and solve the Laplace equation over the axi-symmetric cross section using a high-order FE analysis. In Figure 4 left the finite element solution ( τ F E ) at polynomial level p = 8 is shown whereas in Figure 4 right the difference between the analytical and FE solution is shown τ -τ F E . As may be noticed τ -τ F E is three and a half orders of magnitude smaller compared to τ , indicating on the correctness of the derived analytical solution. If only terms up to (ρ/R) 2 are applied on the boundary of the domain, then the error τ -τ F E increases by one order of magnitude as expected.

2.1.2.

A specific example problem -penny-shaped crack with axisymmetric loading and homogeneous Dirichlet BCs. Similarly to subsection 2.1.1 we present here the first terms in the asymptotic series solution for a penny-shaped crack, ϕ 1 = -π , ω = 2π in an axisymmetric domain with homogeneous Dirichlet BCs [START_REF] Omer | Edge flux intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF]. Now the eigen-pairs are given by We obtain the following expression for the first terms in the asymptotic series solution:

k α k φ k,0 1 1 2 cos ϕ 2 2 1 sin ϕ 3 3 2 cos 3ϕ 2 TABLE 2.
τ = A 1 ρ 1 2 cos ϕ 2 - ρ R 1 4 cos ϕ 2 + ρ R 2 1 12 cos ϕ 2 + 3 32 cos 3ϕ 2 + ρ R 3 - 1 16 cos ϕ 2 - 1 30 cos 3ϕ 2 - 5 128 cos 5ϕ 2 + • • • + A 2 ρ sin ϕ + A 3 ρ 3 2 cos 3ϕ 2 - ρ R 1 4 cos ϕ 2 + ρ R 2 3 32 cos ϕ 2 + 1 10 cos 3ϕ 2 + ρ R 3 - 3 40 cos ϕ 2 - 5 128 cos 3ϕ 2 - 3 70 cos 5ϕ 2 + • • • + • • •
Here we still enforce the orthogonality conditions (20) in order to have uniqueness.

The first terms in the asymptotic solution for the same specific problem are provided in [10, p. 293], where the angular coordinate is measured from the crack face, denoted by ϕ ∈ (0, 2π) , and R = 1 :

τ = c 1 ρ 1 2 sin ϕ 2 - ρ 4 sin ϕ 2 - ρ 2 32 sin 3 ϕ 2 + • • • . ( 21
)
If we replace ϕ by ϕ + π and c 1 by -A 1 , one obtains:

τ = A 1 ρ 1 2 cos ϕ 2 - ρ 4 cos ϕ 2 + ρ 2 32 cos 3ϕ 2 + • • • (22)
whereas our formula ( 46) is (with R = 1 )

τ = A 1 ρ 1 2 cos ϕ 2 - ρ 4 cos ϕ 2 + ρ 2 1 12 cos ϕ 2 + 3 32 cos 3ϕ 2 + • • • (23)
So we note a discrepancy between ( 22) and ( 23) at the level of the third term in the asymptotics. The source of this error in [START_REF] Von Petersdorff | Singularities of the solution of the Laplacian in domains with circular edges[END_REF] is a typo in Proposition 1 where instead of 3 32 the term 1 32 appears. In the proof of this proposition, the right expression is given in equation ( 32) in [START_REF] Von Petersdorff | Singularities of the solution of the Laplacian in domains with circular edges[END_REF], however again in the example at the end of that paper the term 1 12 sin ϕ 2 was forgotten. To verify numerically these formulas, we considered a circular crack in the same torus as in example problem in subsection 2.1.1 with R = 1, ρ out = 1/2 and constructed an axisymmetric FE model. Taking A 1 = 1 and A i = 0 , i = 2, 3, • • • , we prescribed on the outer domain boundary condition either according to ( 22) or (23) with one term (up to ρ 1 2 ), two terms (up two ρ 3/2 ) or three terms (up two ρ 5/2 ). We then computed the discrete L 2 norm of the relative difference between the FE solution and the anticipated "exact solution". Of course that both ( 22) and ( 23) provide the same results if up to two terms in the expansion are considered (which are the same), and the relative difference with three terms is different.

The relative difference is defined as:

||e|| 2 L 2 = 2π π -π ρout 0 |τ -τ F E | 2 × ρ(R + ρ cos ϕ) dρdϕ ||τ || 2 L 2 with ||τ || 2 L 2 = 2π π -π ρout 0 |τ | 2 × ρ(R + ρ cos ϕ) dρdϕ
The FE solution converged to an estimated 0.2% relative error in energy norm, the and integration was performed numerically using 90 Gauss point. The obtained relative errors are shown in Table 3.

number of terms ||e|| 2 L 2 ||τ || 2 L 2 || e || 2 L 2 || τ || 2 L 2 1 1.2 × 10 -3 1.554 × 10 -1 1.2 × 10 -3 1.554 × 10 -1 2 9.5 × 10 -5 1.274 × 10 -1 9.5 × 10 -5 1.274 × 10 -1 3 1.77 × 10 -5 1.317 × 10 -1 7.74 × 10 -5 1.277 × 10 -1 TABLE 3. L 2
relative error -verification of our solution compared to [START_REF] Von Petersdorff | Singularities of the solution of the Laplacian in domains with circular edges[END_REF].

Inspecting the values in Table 3 one may notice that the relative error of our solution decreases by a factor of ≈ 4 when adding the third term, compared to a factor of less than 1.2 for the solution in [START_REF] Von Petersdorff | Singularities of the solution of the Laplacian in domains with circular edges[END_REF], indicating that (23) is the right solution.

A specific example problem -circumferential crack with axisymmetric loading and homogeneous

Neumann BCs. Similarly to subsection 2.1.1 we present herein the first terms in the asymptotic series solution for a circumferential crack, (see Figure 2 (d)) ϕ 1 = -π 2 , ω = 2π in an axisymmetric domain with homogeneous Neumann BCs [START_REF] Von Petersdorff | Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron[END_REF], still taking the orthogonality condition (20) into account.

τ = A 0 + A 1 ρ 1 2 sin ϕ 2 -cos ϕ 2 + ρ R 1 4 sin ϕ 2 + cos ϕ 2 + 1 12 sin 3ϕ 2 -cos 3ϕ 2 + • • • (24)
Note that using the modified angular variable ϕ := ϕ -π 2 ∈ (-π, π) we obtain the expression

τ = A 0 + A 1 ρ 1 2 sin ϕ 2 + ρ R 1 4 cos ϕ 2 + 1 12 cos 3 ϕ 2 + • • • (25)
which can be compared with (19). More generally, if ω = 2π , using the angular variable ϕ := ϕ -ϕ 1 -π ∈ (-π, π) we obtain

τ = A 0 + A 1 ρ 1 2 sin ϕ 2 - ρ R 1 4 cos ϕ 1 sin ϕ 2 + 1 4 sin ϕ 1 cos ϕ 2 + 1 12 sin ϕ 1 cos 3 ϕ 2 + • • • (26)
As particular cases ( ϕ 1 = -π and ϕ 1 = 0 ) we find the penny-shaped and the outer circular crack, for which formula (26) coincides with the formula given in [START_REF] Bach | Propagation of a penny shaped crack under the Irwin criterion[END_REF]:

K 1 (s)ρ 1 2 sin ϕ 2 + 1 4 κ(s)ρ sin ϕ 2 + • • • (27)
since in the first case κ(s) = 1 R and in the second case κ(s) = -1 R .

2.2. General case. If no axi-symmetric assumption is imposed on the data (only the edge is circular) then the full Laplace operator 3D in (3) has to be considered. Like for ( 5)-( 6), we find that solutions τ of (1) are equivalently the solutions of ( r R ) 2 ρ 2 3D τ = 0 , i.e.

(

+ ρ R cos ϕ) 2 (ρ∂ ρ ) 2 + ∂ ϕϕ τ + ρ R (1 + ρ R cos ϕ) cos ϕ(ρ∂ ρ ) -sin ϕ∂ ϕ τ + ρ R 2 ∂ θθ τ = 0. ( 1
) 28
To condense formulas, let us introduce the operators

m 0 (ρ∂ ρ ; ∂ ϕ ) = (ρ∂ ρ ) 2 + ∂ ϕϕ , m 01 (ρ∂ ρ ; ∂ ϕ ) = cos ϕ(ρ∂ ρ ) -sin ϕ∂ ϕ . ( 29
)
Then equation ( 28) is equivalent to

m 0 τ + ρ R 2 cos ϕ m 0 + m 01 τ + ρ R 2 cos 2 ϕ m 0 + cos ϕ m 01 + ∂ θθ τ = 0. (30)
In the general case for a circular edge the following form of expansion series is appropriate:

τ = =0,2,4,... k=0 ∂ θ A k (θ) ρ α k ρ R ∞ i=0 ρ R i φ ,k,i (ϕ) (31)
Remark 5. Notice that φ 0,k,i = φ k,i (associated with the curvature for an axisymmetric case), so these are known for the axi-symmetric analysis.

Comparing this asymptotic expansion to the case of a straight edge [START_REF] Costabel | A quasidual function method for extracting edge stress intensity functions[END_REF], one notices one extra sum, implying that for each primal eigen-function there are two levels of shadow-functions -one set is associated with the derivatives of A k (the index ), and the other set associated with the "curvature terms", i.e. the powers ρ/R (index i ).

The splitting in (30) provides an elegant and convenient way to the formulation of the series expansion of the solution. Introducing the definition for a general term in the expansion (31):

Φ ,k,i def = ρ α k ρ R +i φ ,k,i (ϕ) (32)
we observe that

m 0 (ρ∂ ρ ; ∂ ϕ )Φ ,k,i (ρ, ϕ) = ρ α k ρ R +i m 0 (α k + + i; ∂ ϕ )φ ,k,i (ϕ), (33)
m 01 (ρ∂ ρ ; ∂ ϕ )Φ ,k,i (ρ, ϕ) = ρ α k ρ R +i m 01 (α k + + i; ∂ ϕ )φ ,k,i (ϕ),
Substituting (31) into (30) one deduces:

0 = A k (θ) × m 0 (α k )φ 0,k,0 (34)
+ ρ R m 0 (α k + 1)φ 0,k,1 + 2 cos ϕ m 0 (α k ) + m 01 (α k ) φ 0,k,0 + ρ R 2 m 0 (α k + 2)φ 0,k,2 + 2 cos ϕ m 0 (α k + 1) + m 01 (α k + 1) φ 0,k,1 + cos 2 ϕ m 0 (α k ) + cos ϕ m 01 (α k ) φ 0,k,0 + • • • + A k (θ) × ρ R 2 m 0 (α k + 2)φ 2,k,0 + φ 0,k,0 + ρ R 3 m 0 (α k + 3)φ 2,k,1 + 2 cos ϕ m 0 (α k + 2) + m 01 (α k + 2) φ 2,k,0 + φ 0,k,1 + ρ R 4 m 0 (α k + 4)φ 2,k,2 + 2 cos ϕ m 0 (α k + 3) + m 01 (α k + 3) φ 2,k,1 + cos 2 ϕ m 0 (α k + 2) + cos ϕ m 01 (α k + 2) φ 2,k,0 + φ 0,k,2 + • • •
Equation (34) has to hold true for any (ρ/R) i , and for any ∂ θ A k , resulting in the following recursive set of ordinary differential equations for the determination of the eigen-functions and shadows φ ,k,i (ϕ) :

m 0 (α k + + i)φ ,k,i = (35) -2 cos ϕ m 0 (α k + + i -1) + m 01 (α k + + i -1) φ ,k,i-1 -cos 2 ϕ m 0 (α k + + i -2) + cos ϕ m 01 (α k + + i -2) φ ,k,i-2 -φ -2,k,i , for = 0, 2, 4, 6, • • • , and i ≥ 0.
Here, by convention, φ 's with negative indices are zero. Equations (35) for = 0 are equivalent to equations (15-17) associated with the axi-symmetric case, and for = 2, 4, 6, • • • results in (36) associated with the non-axi-symmetric case.

= 0

Equations (15)-( 17) for the axi-symmetric case hold.

= 2, 4, 6 • • • , i ≥ 0 (36) (α k + i + ) 2 φ ,k,i + φ ,k,i = -( + i + α k -1) [2( + i + α k ) -1] cos ϕ φ ,k,(i-1) + sin ϕ φ ,k,(i-1) -2 cos ϕ φ ,k,(i-1) -( + α k + i -2)( + α k + i -1) cos 2 ϕ φ ,k,(i-2) + cos ϕ sin ϕ φ ,k,(i-2) -cos 2 ϕ φ ,k,(i-2) -φ ( -2),k,i
Equations (36) are complemented by the homogeneous Dirichlet or Neumann boundary conditions:

φ ,k,i (ϕ) = 0, (ϕ = ϕ 1 , ϕ 1 + ω) in case of Drichlet BCs (37) (φ ,k,i ) (ϕ) = 0, (ϕ = ϕ 1 , ϕ 1 + ω) in case of Neumann BCs (38)
2.2.1. A specific example problem -penny-shaped crack for a non-axisymmetric loading and homogeneous Neumann BCs. Again we consider as an example problem a penny-shaped crack ϕ 1 = -π , ω = 2π , however, the loading may by non-axisymmetric, as well as the outer boundary of the 3-D domain of interest. The eigen-functions and a part of the shadow-functions, φ 0,k,i (ϕ) have been provided by (19).As an example, the solution of φ 2,1,0 (ϕ) ( = 2, k = 1, i = 0) may be obtained from (36), for k = 1 α 1 = 1/2 and i = 0, = 2 . All φ s with negative indices in the RHS vanish except one:

1 2 + 0 + 2 2 φ 2,1,0 + φ 2,1,0 = -φ 0,1,0 , (39)
and the homogeneous Neumann BCs read:

φ 2,1,0 (ϕ = ±π) = 0.
From (19), φ 0,1,0 = sin ϕ 2 , thus the solution of (39) can be taken as the particular solution alone:

φ 2,1,0 = - 1 6 sin ϕ 2 (40)
Once φ 2,1,0 is available one may proceed to the computation of φ 2,1,1 (ϕ) ( = 2, k = 1, i = 1) obtained from (36), for k = 1 α 1 = 1/2 and i = 1, = 2 :

1 2 + 1 + 2 2 φ 2,1,1 + φ 2,1,1 = -2 + 1 + 1 2 -1 2 2 + 1 + 1 2 -1 cos ϕφ 2,1,0 + sin ϕφ 2,1,0 -2 cos ϕφ 0,1,0 -φ 0,1,1 (41)
Substituting φ 2,1,0 from (40) and φ 0,1,1 = 1 4 sin ϕ 2 from (19), the particular solution to (41) that satisfies the homogeneous Neumann BCs is:

φ 2,1,1 = - 1 8 sin ϕ 2 + 7 60 sin 3ϕ 2 (42)
This procedure may be continued, to finally obtain the terms in the series expansion:

τ = A 0 (θ) (43)
+ A 0 (θ) ρ R 2 - 1 4 + ρ R 5 16 cos ϕ - ρ R 2 19 128 + 11 64 cos 2ϕ + • • • + • • • + A 1 (θ)ρ 1 2 sin ϕ 2 + ρ R 1 4 sin ϕ 2 + ρ R 2 1 12 sin ϕ 2 - 3 32 sin 3ϕ 2 + + ρ R 3 1 16 sin ϕ 2 - 1 30 sin 3ϕ 2 + 5 128 sin 5ϕ 2 + • • • + A 1 (θ)ρ 1 2 ρ R 2 - 1 6 sin ϕ 2 + - 1 8 sin ϕ 2 + 7 60 sin 3ϕ 2 ρ R + • • • + • • •
Again, the factors corresponding to A 0 (and all terms of even order) and their derivatives are polynomial in (z 1 , z 2 ) . Remark 6. In the vicinity of a crack with a straight edge along the axis z 3 the solution admits the expansion:

τ = A 0 (z 3 ) + A 0 (z 3 )r 2 - 1 4 + • • • (44) +A 1 (z 3 )r 1 2 sin ϕ 2 + A 1 (z 3 )r 5 2 - 1 6 sin ϕ 2 + • • •
One may notice that (44) is composed of the same leading terms associated with i = 0 as in the expansion (43), as expected.

To assess the correctness of the solution (43) we consider a torus with an inner radius r 1 = 9 and an outer radius r 2 = 11 having a circular crack with the tip at R = 10 shown in Figure 3. This time we apply non-axisymmetric boundary conditions, so a fully 3-D FE model is constructed as presented in Figure 5. Taking A 1 = 10 cos θ and A k = 0 , k = 1 , (notice that ρ/R = 1/10 in the considered with homogeneous Neumann boundary conditions on the crack face. In Figure 6 left the finite element solution ( τ F E ) at polynomial level p = 8 is shown whereas in the right the difference between the analytical and FE solution is shown. As may be noticed τ -τ F E is three orders of magnitude smaller compared to τ , assuring the correctness of the derived analytical solution.

2.2.2.

A specific example problem -penny-shaped crack for a non-axisymmetric loading and homogeneous Dirichlet BCs. As an example problem the first terms of the asymptotic solution for a pennyshaped crack ϕ 1 = -π , ω = 2π with homogeneous Dirichlet boundary conditions is provided:

τ = A 1 (θ)ρ 1 2 cos ϕ 2 - 1 4 cos ϕ 2 ρ R + 1 12 cos ϕ 2 + 3 32 cos 3ϕ 2 ρ R 2 + (46) - 1 16 cos ϕ 2 - 1 30 cos 3ϕ 2 - 5 128 cos 5ϕ 2 ρ R 3 + • • • + A 1 (θ)ρ 1 2 ρ R 2 - 1 6 cos ϕ 2 + 1 8 cos ϕ 2 + 7 60 cos 3ϕ 2 ρ R + • • • + • • • 2.2.
τ = A 0 + A 0 ρ R 2 - 1 4 (47) +A 1 ρ 2/3 sin 2ϕ 3 - 1 √ 3 cos 2ϕ 3 + ρ R 1 60 5 √ 3 cos ϕ 3 - √ 3 cos 5ϕ 3 + 15 sin ϕ 3 + 3 sin 5ϕ 3 + ρ R 2 1 160 12 sin 2ϕ 3 -4 √ 3 cos 2ϕ 3 -15 sin 4ϕ 3 -5 √ 3 cos 4ϕ 3 + +A 1 ρ 2/3 ρ R 2 1 20 √ 3 cos 2ϕ 3 -3 sin 2ϕ 3 + • • • 2.2.4.
A specific example problem -exterior circular crack with non-axisymmetric loading and homogeneous Neumann BCs. Consider the circular external crack as in Figure 2 (b) with homogeneous Neumann BCs. In this case ϕ ∈ (0, 2π) and the first few terms in the asymptotic solution are given by:

τ = A 0 + A 0 ρ R 2 - 1 4 (48) +A 1 ρ 1/2 cos ϕ 2 - ρ R 1 4 cos ϕ 2 + ρ R 2 1 12 cos ϕ 2 + 3 32 cos 3ϕ 2 +A 1 ρ 1/2 ρ R 2 - 1 6 cos ϕ 2 + • • •

ASYMPTOTIC SOLUTION FOR THE ELASTICITY SYSTEM

The point of departure for the elasticity system is the three equilibrium equations given in a Cartesian system by:

3 i=1 ∂σ ij ∂x i + F j = 0, j = 1, 2, 3,
where σ ij are the components of the Cartesian stress tensor, and F j are the body forces. For vanishing body forces in the vicinity of the singular edge, i.e. F 1 = F 2 = F 3 = 0 we reformulate the equilibrium equations to be expressed in the coordinate system ρ, ϕ, θ shown in Figure 1:

x 1 = (R + ρ cos ϕ) cos θ, x 2 = (R + ρ cos ϕ) sin θ, x 3 = ρ sin ϕ (49)
The Cartesian components of the displacements (u 1 , u 2 , u 3 ) are connected to the local coordinate system displacements (u ρ , u ϕ , u θ ) by:

u 1 = (u ρ cos ϕ -u ϕ sin ϕ) cos θ -u θ sin θ (50) u 2 = (u ρ cos ϕ -u ϕ sin ϕ) sin θ + u θ cos θ u 3 = u ρ sin ϕ + u ϕ cos ϕ
The kinematic connection between displacements and strains in the coordinate system ρ, ϕ, θ can be derived using the metric tensor [START_REF] Sokolnikoff | Mathematical Theory of Elasticity[END_REF]:

ε ρρ = ∂u ρ ∂ρ (51)
ε ϕϕ = 1 ρ ∂u ϕ ∂ϕ + u ρ ρ ε θθ = 1 r ∂u θ ∂θ + u ρ cos ϕ -u ϕ sin ϕ ε ρϕ = 1 2 1 ρ ∂u ρ ∂ϕ + u ϕ ρ - u ϕ ρ ε ϕθ = 1 2 1 r ∂u ϕ ∂θ + 1 ρ ∂u θ ϕ + 1 r u θ sin ϕ ε ρθ = 1 2 1 r ∂u ρ ∂θ + ∂u θ ρ - 1 r u θ cos ϕ
These connections are identical to [7, Eq. ( 3)]. In the new curvilinear coordinate system the equilibrium equations in ρ, ϕ, θ directions read [START_REF] Sokolnikoff | Mathematical Theory of Elasticity[END_REF]:

0 = ∂σ ρρ ∂ρ + 1 ρ ∂σ ρϕ ∂ϕ + σ ρρ -σ ϕϕ ρ + 1 r ∂σ ρθ ∂θ + (σ ρρ -σ θθ ) cos ϕ -σ ρϕ sin ϕ (52) 0 = 1 ρ ∂σ ϕϕ ∂ϕ + ∂σ ρϕ ∂ρ + 2 ρ σ ρϕ + 1 r ∂σ ϕθ ∂θ + (σ θθ -σ ϕϕ ) sin ϕ + σ ρϕ cos ϕ (53) 0 = ∂σ ρθ ∂ρ + 1 ρ σ ρθ + 1 ρ ∂σ ϕθ ∂ϕ + 1 r ∂σ θθ ∂θ + 2σ ρθ cos ϕ -2σ ϕθ sin ϕ . (54)
We consider an isotropic elastic homogeneous material with Lamé constants λ and µ (kinematic equations and Hooke's law given by:

σ ij = λδ ij 3 k=1 ε kk + 2µε ij , i, j = ρ, ϕ, θ (55)
The stresses in terms of displacements are given by: 56) into ( 52)-( 54) and multiplying by r 2 (ρ/R) 2 = ρ 2 1 + ρ R cos ϕ 2 one obtains the Navier-Lamé system:

        σ ρρ σ θθ σ ϕϕ σ ρθ σ ρϕ σ θϕ         =           λ 1 ρ + (λ + 2µ)∂ ρ λ 1 ρ ∂ ϕ 0 λ 1 ρ + ∂ ρ ) λ 1 ρ ∂ ϕ 0 (λ + 2µ) 1 ρ + λ∂ ρ (λ + 2µ) 1 ρ ∂ ϕ 0 0 0 µ∂ ρ µ 1 ρ ∂ ϕ µ -1 ρ + ∂ ρ 0 0 0 µ 1 ρ ∂ ϕ             u ρ u ϕ u θ   + 1 r         λ cos ϕ -λ sin ϕ 0 (λ + 2µ) cos ϕ -(λ + 2µ) sin ϕ 0 λ cos ϕ -λ sin ϕ 0 0 0 -µ cos ϕ 0 0 0 0 0 µ sin ϕ           u ρ u ϕ u θ   + 1 r         0 0 λ 0 0 λ + 2µ 0 0 λ µ 0 0 0 0 0 0 µ 0         ∂ θ   u ρ u ϕ u θ   (56) Inserting (
0 = 1 + ρ R cos ϕ 2 (λ + 2µ) (ρ∂ ρ ) 2 -1 u ρ + µ ∂ ϕϕ u ρ -(λ + 3µ)∂ ϕ u ϕ + (λ + µ)∂ ρϕ u ϕ (57) + 1 + ρ R cos ϕ ρ R [(λ + 2µ) cos ϕ ρ ∂ ρ u ρ -(λ + µ) sin ϕ ρ ∂ ρ u ϕ + µ sin ϕ (u ϕ -∂ ϕ u ρ ) + (λ + µ)ρ ∂ ρθ u θ ] + ρ R 2 [(λ + 2µ) cos ϕ(u ϕ sin ϕ -u ρ cos ϕ) + µ∂ θθ u ρ -(λ + 3µ) cos ϕ ∂ θ u θ ] 0 = 1 + ρ R cos ϕ 2 µ (ρ∂ ρ ) 2 -1 u ϕ + (λ + 2µ)∂ ϕϕ u ϕ + (λ + 3µ)∂ ϕ u ρ + (λ + µ)ρ ∂ ρϕ u ρ (58) + 1 + ρ R cos ϕ ρ R [(λ + µ) cos ϕ (∂ ϕ u ρ -u ϕ ) + µ cos ϕ ρ∂ ρ u ϕ -(λ + 2µ) sin ϕ (∂ ϕ u ϕ + u ρ ) + (λ + µ)∂ ϕθ u θ ] + ρ R 2 [(λ + 2µ) sin ϕ(u ρ cos ϕ -u ϕ sin ϕ) + µ∂ θθ u ϕ + (λ + 3µ) sin ϕ∂ θ u θ ] 0 = 1 + ρ R cos ϕ 2 µ (ρ∂ ρ ) 2 + ∂ ϕϕ u θ (59) + 1 + ρ R cos ϕ ρ R [µ (cos ϕ ρ∂ ρ u θ -sin ϕ ∂ ϕ u θ ) + (λ + µ) (∂ θ u ρ + ∂ ϕθ u ϕ + ρ∂ ρθ u ρ )] + ρ R 2 [-µu θ + (λ + 2µ)∂ θθ u θ + (λ + 3µ) (cos ϕ ∂ θ u ρ -sin ϕ ∂ θ u ϕ )]
The Navier-Lamé equations are complemented by homogeneous boundary conditions on the faces intersecting at the singular edge:

u ρ = u ϕ = u θ = 0 on Γ 1 ∪ Γ 2 Clamped BCs (60) t ϕ = t ρ = t θ = 0 on Γ 1 ∪ Γ 2 Traction Free BCs, ( 61
)
where t is the traction vector on the boundary. On the boundaries Γ 1 , Γ 2 , i.e. for ϕ = ϕ 1 , ϕ 2 , the traction free BCs (61) are expressed in terms of the stresses using Cauchy's law:

σ ϕϕ = σ ρϕ = σ θϕ = 0, ϕ = ϕ 1 , ϕ 2
Denoting the displacement vector by u = (u ρ , u ϕ , u θ ) , the Navier-Lamé (57-59) system is split as follows:

1 + ρ R cos ϕ 2 [M 0 ] + 1 + ρ R cos ϕ ρ R [M 01 ] + ρ R 2 [M 02 ] (62) + 1 + ρ R cos ϕ ρ R [M 10 ]∂ θ + ρ R 2 [M 11 ]∂ θ + ρ R 2 [M 2 ]∂ θθ u = 0, with [M 0 ] =   (λ + 2µ) (ρ∂ ρ ) 2 -1 + µ∂ ϕϕ -(λ + 3µ)∂ ϕ + (λ + µ)ρ∂ ρϕ 0 (λ + µ)ρ∂ ρϕ + (λ + 3µ)∂ ϕ (λ + 2µ)∂ ϕϕ + µ (ρ∂ ρ ) 2 -1 0 0 0 µ (ρ∂ ρ ) 2 + ∂ ϕϕ   , (63)
[M 01 ] =   (λ + 2µ) cos ϕ ρ∂ ρ -µ sin ϕ∂ ϕ sin ϕ [-(λ + µ)ρ∂ ρ + µ] 0 (λ + µ) cos ϕ ∂ ϕ -(λ + 2µ) sin ϕ cos ϕ [-(λ + µ) + µρ∂ ρ ] -(λ + 2µ) sin ϕ∂ ϕ 0 0 0 µ (cos ϕ ρ∂ ρ -sin ϕ∂ ϕ )   , (64)
[M 02 ] =   -(λ + 2µ) cos 2 ϕ (λ + 2µ) sin ϕ cos ϕ 0 (λ + 2µ) sin ϕ cos ϕ -(λ + 2µ) sin 2 ϕ 0 0 0 -µ   , (65)
[M 10 ] =   0 0 (λ + µ)ρ∂ ρ 0 0 (λ + µ)∂ ϕ (λ + µ) (ρ∂ ρ + 1) (λ + µ)∂ ϕ 0   , (66)
[M 11 ] =   0 0 -(λ + 3µ) cos ϕ 0 0 (λ + 3µ) sin ϕ (λ + 3µ) cos ϕ -(λ + 3µ) sin ϕ 0   , ( 67
)
[M 2 ] =   µ 0 0 0 µ 0 0 0 (λ + 2µ)   . ( 68
)
Following same asymptotic series expansion shown to be appropriate for the Laplace operator, we assume herein an expansion for the displacements of the form:

u = =0 k=0 ∂ θ A k (θ) ρ α k ∞ i=0 ρ R i+    φ ρ (ϕ) φ ϕ (ϕ) φ θ (ϕ)    ,k,i = =0 k=0 ∂ θ A k (θ) ρ α k ∞ i=0 ρ R i+ φ ,k,i (69)
Comparing this asymptotic expansion to the case of a straight edge, one notices one extra sum, implying that for each primal eigen-function there are two sets of shadow-functions -one set is associated with the derivatives of A k , and the other set associated with the "curvature terms", i.e. the powers ρ/R . Inserting (69) in (62) so to gather terms of same order of derivatives of A k , and same order of powers ρ/R , we obtain the following recursive formula for the computation of the primal and shadow functions:

[m 0 ]φ ,k,i = -(2 cos ϕ[m 0 ] + [m 01 ]) φ ,k,i-1 -cos 2 ϕ[m 0 ] + cos ϕ[m 01 ] + [m 02 ] φ ,k,i-2 -[m 10 ]φ -1,k,i -(cos ϕ[m 10 ] + [m 11 ]) φ -1,k,i-1 -[m 2 ]φ -2,k,i , ≥ 0, i ≥ 0 (70)
where φ 's with negative indices are set to zero, and

[m 0 ]φ ,k,i =   (λ + 2µ) β 2 -1 + µ∂ ϕϕ ((λ + µ)β -(λ + 3µ)) ∂ ϕ 0 ((λ + µ)β + (λ + 3µ)) ∂ ϕ µ β 2 -1 + (λ + 2µ)∂ ϕϕ 0 0 0 µ β 2 + ∂ ϕϕ   φ ,k,i (71)
[m 01 ]φ ,k,i =   (λ + 2µ) cos ϕβ -µ sin ϕ∂ ϕ sin ϕ (µ -(λ + µ)β) 0 -(λ + 2µ) sin ϕ + (λ + µ) cos ϕ∂ ϕ cos ϕ (µ (β -1) -λ) -(λ + 2µ) sin ϕ∂ ϕ 0 0 0 µ (β cos ϕ -sin ϕ∂ ϕ )   φ ,k,i (72)
[m 02 ]φ ,k,i =   -(λ + 2µ) cos 2 ϕ (λ + 2µ) cos ϕ sin ϕ 0 (λ + 2µ) sin ϕ cos ϕ -(λ + 2µ) sin 2 ϕ 0 0 0 -µ   φ ,k,i (73)
[m 10 ]φ ,k,i =   0 0 (λ + µ)β 0 0 (λ + µ)∂ ϕ (λ + µ)β (λ + µ)∂ ϕ 0   φ ,k,i (74)
[m 11 ]φ ,k,i =   0 0 -(λ + 3µ) cos ϕ 0 0 (λ + 3µ) sin ϕ (λ + 3µ) cos ϕ -(λ + 3µ) sin ϕ 0   φ ,k,i (75)
[m 2 ]φ ,k,i =   µ 0 0 0 µ 0 0 0 (λ + 2µ)   φ ,k,i (76)
where β = (α k + + i).

φ ,k,i (ϕ 1 ) = φ ,k,i (ϕ 2 ) = 0, ∀ , α k , i (77)
The traction free boundary conditions imply:

(σ ϕϕ , σ ρϕ , σ θϕ ) = 0, ϕ = ϕ 1 , ϕ 2 (78)
Multiplying (78) by r ρ R yields:

1 + ρ R cos ϕ [T 0 ] + ρ R [T 01 ] + ρ R [T 1 ]∂ θ u = 0, ϕ = ϕ 1 , ϕ 2 (79)
where,

[T 0 ] =   (λ + 2µ) + λρ∂ ρ (λ + 2µ)∂ ϕ 0 µ∂ ϕ -µ(1 -ρ∂ ρ ) 0 0 0 µ∂ ϕ   , (80)
[T 01 ] =   λ cos ϕ -λ sin ϕ 0 0 0 0 0 0 µ sin ϕ   , [T 1 ] =   0 0 λ 0 0 0 0 µ 0   (81)
Inserting ( 69) in ( 79) we obtain the following boundary conditions on the primal and shadow functions ∀ , α k , i :

[t 0 ]φ ,k,i = -(cos ϕ[t 0 ] + [t 01 ]) φ ,k,i-1 -[t 1 ]φ -1,k,i , ϕ = ϕ 1 , ϕ 2 (82)
again φ 's with negative indices are set to zero, and

[t 0 ]φ ,k,i =   2µ + λ (β + 1) (λ + 2µ)∂ ϕ 0 µ∂ ϕ µ (β -1) 0 0 0 µ∂ ϕ   φ ,k,i (83)
[t 01 ]φ ,k,i =   λ cos ϕ -λ sin ϕ 0 0 0 0 0 0 µ sin ϕ   φ ,k,i , [t 1 ]φ ,k,i =   0 0 λ 0 0 0 0 µ 0   φ ,k,i (84)
3.2. Axi-symmetric case. For an axi-symmetric domain and boundary conditions all derivatives with respect to θ vanish so (62) is simplified to:

1 + ρ R cos ϕ 2 [M 0 ] + 1 + ρ R cos ϕ ρ R [M 01 ] + ρ R 2 [M 02 ] u = 0,
Notice that third equation in the set (85) is in terms of u θ alone and decoupled from the first two equations, thus the solution to the displacement u θ is decoupled from u ρ and u ϕ for the axi-symmetric case.

Remark 7. The equation for determining u θ is very similar (but not identical) to the Laplace equation (compare to (4)). The difference is the last term 1 r 2 u θ . This term results in a slightly different solution of u θ compared to the Laplace solution only when inspecting higher shadows ( i ≥ 2 ).

The asymptotic series expansion in case of an axi-symmetric solution is obtained by taking = 0 in (69):

u = ∞ k=0 A k ρ α k ∞ i=0 ρ R i    φ ρ (ϕ) φ ϕ (ϕ) φ θ (ϕ)    0,k,i (85)
In the axi-symmetric case, ≡ 0 , therefore the recursive formulas (70) and (82) are simplified to:

[m 0 ]φ 0,k,i = -(2 cos ϕ[m 0 ] + [m 01 ]) φ 0,k,i-1 -cos 2 ϕ[m 0 ] + cos ϕ[m 01 ] + [m 02 ] φ 0,k,i-2 i ≥ 0 (86)
and

[t 0 ]φ 0,k,i = -(cos ϕ[t 0 ] + [t 01 ]) φ 0,k,i-1 , ϕ = ϕ 1 , ϕ 2 (87)
3.2.1. A specific example problem -penny-shaped crack with axisymmetric loading and traction free BCs. In the case of a penny shaped crack ω = 2π , α 1 = α 2 = α 3 = 0 and α 4 = α 5 = α 6 = 1 2 , the expressions for the traction free boundary condition (87) are further simplified:

[t 0 ]Φ ,k,i =   2µ + λ (α k + + i + 1) (λ + 2µ)∂ ϕ 0 µ∂ ϕ µ (α k + + i -1) 0 0 0 µ∂ ϕ   φ ,k,i , (88)
[t 01 ]Φ ,k,i =   -λ 0 0 0 0 0 0 0 0   φ ,k,i , [t 1 ]Φ ,k,i =   0 0 λ 0 0 0 0 µ 0   φ ,k,i (89)
The solution to the system (86) with traction free BCs (87) for a penny shaped crack ( ω = 2π ) is:

   u ρ u ϕ u θ    = A 1   cos ϕ -sin ϕ 0   1 - ρ R cos ϕ + ρ R 2 cos 2 ϕ + • • • + A 2   sin ϕ cos ϕ 0   + A 3   0 0 1   1 + ρ R cos ϕ + A 4 ρ 1 2        -λ+5µ λ+µ cos ϕ 2 + cos 3ϕ 2 3λ+7µ λ+µ sin ϕ 2 -sin 3ϕ 2 0    + ρ R     3λ 2 +22λµ-13µ 2 12(λ+µ) 2 cos ϕ 2 + 3 4 -λ λ+µ cos 3ϕ 2 9λ 2 -14λµ-55µ 2 12(λ+µ) 2 sin ϕ 2 + -5 4 + λ λ+µ sin 3ϕ 2 0     + ρ R 2    -3λ 2 +134λµ+67µ 2 96(λ+µ) 2 cos ϕ 2 -45λ 3 +255λ 2 µ+199λµ 2 +53µ 3 180(λ+µ) 3 cos 3ϕ 2 + 3(3λ-µ) 32(λ+µ) cos 5ϕ 2 69λ 2 +170λµ+37µ 2 96(λ+µ) 2 sin ϕ 2 + -15λ 3 +155λ 2 µ+371λµ 2 +137µ 3 180(λ+µ) 3 sin 3ϕ 2 -3(λ-3µ) 32(λ+µ) sin 5ϕ 2 0    + • • •    + A 5 ρ 1 2       -λ+5µ 3(λ+µ) sin ϕ 2 + sin 3ϕ 2 -3λ+7µ 3(λ+µ) cos ϕ 2 + cos 3ϕ 2 0    + ρ R    -123λ 2 +262λµ+107µ 2 180(λ+µ) 2 sin ϕ 2 -λ-3µ 12(λ+µ) sin 3ϕ 2 -15λ 2 +62λµ+79µ 2 180(λ+µ) 2 cos ϕ 2 + (λ+5µ) 12(λ+µ) cos 3ϕ 2 0    + ρ R 2 1 1440    -483λ 2 +1350λµ+803µ 2 (λ+µ) 2 sin ϕ 2 + 16 2625λ 3 +8779λ 2 µ+9251λµ 2 +3161µ 3 35(λ+µ) 3 sin 3ϕ 2 + 45 3λ-µ (λ+µ) sin 5ϕ 2 -16 357λ 2 +1002λµ+581µ 2 (λ+µ) 2 cos ϕ 2 + 1365λ 3 +5271λ 2 µ+6719λµ 2 +2749µ 3 35(λ+µ) 3 cos 3ϕ 2 + 45 λ-3µ (λ+µ) cos 5ϕ 2 0    + • • • ] + A 6 ρ 1 2     0 0 sin ϕ 2   + ρ R   0 0 1 4 sin ϕ 2   + ρ R 2   0 0 1 4 sin ϕ 2 -3 32 sin 3ϕ 2   + • • •   + • • •
One may notice that the solutions of φ ρ and φ ϕ are decoupled from the solution of φ θ . The eigenfunctions φ θ are not identical to the ones obtained by the Laplace equation -if one observes the last line in (90) in comparison with (19), then the coefficient of the term ρ 1/2 ρ R 2 sin ϕ 2 is 1 4 instead of 1 12 in (19). More differences would appear if inspecting terms of higher orders of ρ R . Remark 8. The second and third terms corresponding to zero eigen-values are associated with rigid body motions: The term A 2 (sin ϕ, cos ϕ, 0) expressed in Cartesian components becomes (u 1 , u 2 , u 3 ) 2 = A 2 (0, 0, 1) , which is clearly a translation along z 3 . Likewise, the factor term of A 3 becomes (u 1 , u 2 , u 3 ) 3 = A 3 r R (-sin θ, cos θ, 0) , which is proportional to the rotation (-x 2 , x 1 , 0) around the axis of symmetry. Note finally that for ρ/R 1

1 1 + ρ R cos ϕ = 1 - ρ R cos ϕ + ρ R cos ϕ 2 - ρ R cos ϕ 3 + • • • therefore the factor of A 1 is not linear in Cartesian coordinates: (u 1 , u 2 , u 3 ) 1 = A 1 R( x 1 r 2 , x 2
r 2 , 0) . Details on the derivation of the series solution associated with A 1 , A 2 , A 3 are provided in Appendix A.

Remark 9. The solution for the circular singular crack with traction free boundary condition was presented also by Leung and Su in [START_REF] Leung | Eigenfunction expansion for penny-shaped and circumferential cracks[END_REF], and is almost identical to (90). One difference is the factor 1/1440 in the term associated with A 2 ρ 1 2 ρ R 2 whereas in [START_REF] Leung | Eigenfunction expansion for penny-shaped and circumferential cracks[END_REF] the factor is 1/1400 . The solution herein is the correct one because it satisfies (86), whereas Leung and Su's solution does not. The second difference is in the shadow terms associated with α 1 = 0 and α 3 = 0 for A 1 and A 3 which are missing in Leung and Su's solution.

Defining:

K 1 def = A 1 µ, -1 4
K I √ 2π def = A 4 µ, 3 4
K II √ 2π def = A 5 µ, 2 K III √ 2π def = A 6 µ (90)
with K I , K II , K III the renown stress intensity factors in fracture mechanics, and using (51)-( 55), the stress tensor associated with an axisymmetric traction free problem is:

               σ ρρ σ θθ σ ϕϕ σ ρθ σ ρϕ σ θϕ                = K 1         -1 -cos 2ϕ 2 -1 + cos 2ϕ 0 sin 2ϕ 0         ρ R 1 -2 cos ϕ ρ R + 3 cos 2 ϕ ρ R 2 -4 cos 3 ϕ ρ R 3 + • • • (91) + -1 4
K I √ 2πρ                   -5 cos ϕ 2 + cos 3ϕ 2 -4λ λ+µ cos ϕ 2 -3 cos ϕ 2 -cos 3ϕ 2 0 -sin ϕ 2 -sin 3ϕ 2 0         + ρ R           -5λ+13µ 4(λ+µ) cos ϕ 2 + λ+9µ 4(λ+µ) cos 3ϕ 2 -2(2λ+µ)(λ+5µ) (λ+µ) 2 cos ϕ 2 + 3λ+2µ λ+µ cos 3ϕ 2 -3(λ+9µ) 4(λ+µ) cos ϕ 2 -λ+9µ 4(λ+µ) cos 3ϕ 2 0 λ-7µ 4(λ+µ) sin ϕ 2 + λ-7µ 4(λ+µ) sin 3ϕ 2 0           + • • •           + 3 4
K II √ 2πρ                    -5 3 sin ϕ 2 + sin 3ϕ 2 -4λ 3(λ+µ) sin ϕ 2 -sin ϕ 2 -sin 3ϕ 2 0 1 3 cos ϕ 2 + 3 cos 3ϕ 2 0          + ρ R           -51λ+107µ 60
0           + • • •           + 2 K III √ 2πρ                 0 0 0 1 2 sin ϕ 2 0 1 2 cos ϕ 2         + ρ R         0 0 0 7 8 sin ϕ 2 -1 2 sin 3ϕ 2 0 5 8 cos ϕ 2 -1 2 cos 3ϕ 2         + • • •        
Remark 10. For R → ∞ the stresses state should tend to a plane-strain state. Indeed, by computing ν(σ ρρ + σ ϕϕ ) from (91), one obtains ν

2 K I √ 2πρ cos ϕ 2 -2 K II √
2πρ sin ϕ 2 that equals σ θθ . This is exactly the connection: σ 33 = ν(σ 11 + σ 22 ) according to a plane-strain situation.

Remark 11. The primal eigen-stresses σ ρρ and σ ϕϕ do not depend on the material properties for traction free boundary conditions on crack faces. However, their shadows do depend on the material properties.

3.2.2.

A specific example problem -penny-shaped crack with axisymmetric loading and clamped BCs. For the clamped BCs (Dirichlet BCs) and a penny shaped crack ( ω = 2π ) the series expansion of the solution is:

   u ρ u ϕ u θ    = A 1 ρ 1 2      λ+5µ 3λ+7µ cos ϕ 2 + cos 3ϕ 2 -sin ϕ 2 -sin 3ϕ 2 0   + ρ R    -41λ 2 +154λµ+129µ 2 4(3λ+7µ)(5λ+9µ) cos ϕ 2 + λ-3µ 4(3λ+7µ) cos 3ϕ 2 λ+5µ 4(3λ+7µ) sin ϕ 2 + λ+5µ 4(3λ+7µ) sin 3ϕ 2 0    (92) + ρ R 2    (7λ+11µ)(23λ+67µ) 32(3λ+7µ)(5λ+9µ) cos ϕ 2 + 25λ 2 +98λµ+89µ 2 4(3λ+7µ)(5λ+9µ) cos 3ϕ 2 + 3(-3λ+µ) 32(3λ+7µ) cos 5ϕ 2 -(7λ+11µ)(17λ+37µ) 32(3λ+7µ)(5λ+9µ) sin ϕ 2 -13λ 2 +58λµ+61µ 2 4(3λ+7µ)(5λ+9µ) sin 3ϕ 2 + 3(λ-3µ) 32(3λ+7µ) sin 5ϕ 2 0    + • • •    + A 2 ρ 1 2      sin ϕ 2 + sin 3ϕ 2 3λ+7µ λ+5µ cos ϕ 2 + cos 3ϕ 2 0   + ρ R    λ-3µ 4(λ+5µ) sin ϕ 2 + λ-3µ 4(λ+5µ) sin 3ϕ 2 -3λ+7µ 4(λ+5µ) cos ϕ 2 -1 4 cos 3ϕ 2 0    + ρ R 2    3λ-µ 96(λ+5µ) sin ϕ 2 + -3λ+µ 12(λ+5µ) sin 3ϕ 2 + 3(-3λ+µ) 32(λ+5µ) sin 5ϕ 2 69λ+137µ 96(λ+5µ) cos ϕ 2 + 1 12 cos 3ϕ 2 + -3(λ-3µ) 32(λ+5µ) cos 5ϕ 2 0    + • • •    + A 3 ρ 1 2     0 0 cos ϕ 2   + ρ R   0 0 -1 4 cos ϕ 2   + ρ R 2   0 0 1 4 cos ϕ 2 + 3 32 cos 3ϕ 2   + • • •  
Here we define:

A 1 def = 3λ + 7µ -4µ K I √ 2π , A 2 def = λ + 5µ 4µ K II √ 2π , A 3 def = 2 K III √ 2π (93)
then the corresponding stress tensor is:

               σ ρρ σ θθ σ ϕϕ σ ρθ σ ρϕ σ θϕ                = 3λ + 7µ -4µ K I √ 2πρ                       5(λ+µ) 3λ+7µ cos ϕ 2 + cos 3ϕ 2 4λ 3λ+7µ cos ϕ 2 3(λ+µ) 3λ+7µ cos ϕ 2 -cos 3ϕ 2 0 λ+µ 3λ+7µ sin ϕ 2 -sin 3ϕ 2 0           + ρ R             - 3(17λ 2 +98λµ+129µ 2 ) 4(3λ+7µ)(5λ+9µ) cos ϕ 2 -λ+9µ 4(3λ+7µ) cos 3ϕ 2 2(34λ 2 +131λµ+117µ 2 ) (3λ+7µ)(5λ+9µ) cos ϕ 2 -3λ+2µ 3λ+7µ cos 3ϕ 2 - (5λ 2 +106λµ+213µ 2 ) 4(15λ 2 +62λµ+63µ 2 ) cos ϕ 2 + λ+9µ 4(3λ+7µ) cos 3ϕ 2 0 5λ 2 +106λµ+213µ 2 4(3λ+7µ)(5λ+9µ) sin ϕ 2 -λ-7µ 4(3λ+7µ) sin 3ϕ 2 0             + • • •             + λ + 5µ 4µ K II √ 2πρ                     5(λ+µ) λ+5µ sin ϕ 2 + sin 3ϕ 2 4λ λ+5µ sin ϕ 2 3(λ+µ) λ+5µ sin ϕ 2 -sin 3ϕ 2 0 -λ+µ λ+5µ cos ϕ 2 + cos 3ϕ 2 0           + ρ R           -5λ+9µ 4(λ+5µ) sin ϕ 2 -λ+9µ 4(λ+5µ) sin 3ϕ 2 -2(2λ+µ) λ+5µ sin ϕ 2 -3λ+2µ λ+5µ sin 3ϕ 2 -3λ+µ 4(λ+5µ) sin ϕ 2 + λ+9µ 4(λ+5µ) sin 3ϕ 2 0 -λ+5µ 4(λ+5µ) cos ϕ 2 + λ-7µ 4(λ+5µ) cos 3ϕ 2 0           + • • •           + 2 K III √ 2πρ                 0 0 0 1 2 cos ϕ 2 0 -1 2 sin ϕ 2         + ρ R         0 0 0 -7 8 cos ϕ 2 -1 2 cos 3ϕ 2 0 5 8 sin ϕ 2 + 1 2 sin 3ϕ 2         + • • •         (94)
3.3. Non axi-symmetric case. The set of equations corresponding to = 0 is exactly the axisymmertric set of equations presented in section 3.2. Herein we present the entire solution including the terms for ≥ 1 for a non axi-symmetric case.

3.3.1.

A specific example problem -penny-shaped crack with non-axisymmetric loading and traction free BCs. For the traction free BCs (82) and a penny shaped crack ( ω = 2π ) we solved the system (70), obtaining: The stresses in this case with the definition of SIFs as in (90) are: Comparing the terms associated with the first derivatives of the SIFs ( K I (θ), K II (θ), K III (θ) ) in [START_REF] Leblond | The stress-field near the front of an arbitrarily shaped crack in a 3-dimensional elastic body[END_REF] with (96), one notices that these are identical for K I (θ) and K III (θ) . The term that multiplies K II (θ) in [START_REF] Leblond | The stress-field near the front of an arbitrarily shaped crack in a 3-dimensional elastic body[END_REF] appears in (96), but in our expression there are another two expressions proportional to cos 3ϕ 2 and sin 3ϕ 2 that are absent in [START_REF] Leblond | The stress-field near the front of an arbitrarily shaped crack in a 3-dimensional elastic body[END_REF]. we see that the series of shadows can be stopped there. In this way we find a global rigid body motion as specified in Remark 8.

   u ρ u ϕ u θ    = A 1 (θ)   cos ϕ -sin ϕ 0   1 - ρ R cos ϕ + ρ R 2 cos 2 ϕ + • • • + A 1 (θ) ρ R   ρ R   0 0 -cos 2 ϕ   + • • •   + A 1 (θ) ρ R 2    -λ+2µ 12(λ+µ) cos ϕ 3λ+2µ 12(λ+µ) sin ϕ 0    + • • • + A 2 (θ)   sin ϕ cos ϕ 0   + A 2 (θ) ρ R     0 0 -sin ϕ   + ρ R   0 0 1 2 sin 2ϕ   + • • •   + A 2 (θ) ρ R 2    λ 4(λ+µ) sin ϕ -λ 4(λ+µ) cos ϕ 0    + • • • + A 3 (θ)   0 0 1   1 + ρ R cos ϕ + A 3 (θ) ρ R      -λ 2(λ+µ) 0 0   + ρ R    9λ 2 +16λµ+8µ 2 24(λ+µ) 2 cos ϕ -15λ 2 +24λµ+8µ 2 24(λ+µ) 2 sin ϕ 0    • • •    + A 3 (θ) ρ R 2   0 0 -1 2   + • • • + A 4 (θ)ρ 1 2        -λ+5µ λ+µ cos ϕ 2 + cos 3ϕ 2 3λ+7µ λ+µ sin ϕ 2 -sin 3ϕ 2 0    + ρ R     3λ 2 +22λµ-13µ 2 12(λ+µ) 2 cos ϕ 2 + 3 4 -λ λ+µ cos 3ϕ 2 9λ 2 -14λµ-55µ 2 12(λ+µ) 2 sin ϕ 2 + -5 4 + λ λ+µ sin 3ϕ 2 0     + ρ R 2    -3λ 2 +134λµ+67µ 2 96(λ+µ) 2 cos ϕ 2 -45λ 3 +255λ 2 µ+199λµ 2 +53µ 3 180(λ+µ) 3 cos 3ϕ 2 + 3(3λ-µ) 32(λ+µ) cos 5ϕ 2 69λ 2 +170λµ+37µ 2 96(λ+µ) 2 sin ϕ 2 + -15λ 3 +155λ 2 µ+371λµ 2 +137µ 3 180(λ+µ) 3 sin 3ϕ 2 -3(λ-3µ) 32(λ+µ) sin 5ϕ 2 0    + • • •    + A 4 (θ)ρ 1 2 ρ R     0 0 2 cos ϕ 2 -2(3λ+7µ) 3(λ+µ) cos 3ϕ 2   + ρ R   0 0 3λ+23µ 6(λ+µ) cos ϕ 2 -3 2 cos 3ϕ 2 + 15λ 2 +34λµ+3µ 2 15(λ+µ) 2 cos 5ϕ 2   + • • •   + A 4 (θ)ρ 1 2 ρ R 2    -3λ-µ 6(λ+µ) cos ϕ 2 + 45λ 2 +138λµ+61µ 2 90(λ+µ) 2 cos 3ϕ 2 -3λ+7µ 6(λ+µ) sin ϕ 2 + 15λ 2 -2λµ-49µ 2 90(λ+µ) 2 sin 3ϕ 2 0    + • • • / • • • (95)
               σ ρρ σ θθ σ ϕϕ σ ρθ σ ρϕ σ θϕ                = K 1 (θ)         -1 -cos 2ϕ 2 -1 + cos 2ϕ 0 sin 2ϕ 0         ρ R + K 1 (θ) ρ R         0 0 0 cos ϕ 0 -sin ϕ         + • • • (96) +K 2 (θ) ρ R 2         0 0 0 sin 2ϕ 0 -2 sin 2 ϕ         + K 3 (θ) ρ R         0 3λ+2µ λ+µ 0 0 0 0         + • • • + -1 4
K I (θ) √ 2πρ                   -5 cos ϕ 2 + cos 3ϕ 0           + • • •           + -1 4
K I (θ) √ 2πρ ρ R          0 0 0 2(λ-µ) λ+µ cos ϕ 2 -2(
0           + • • •           + 3 4
K II (θ) √ 2πρ ρ R          0 0 0 2(λ-µ) 3(λ+µ) sin ϕ 2 + sin 3ϕ 2 0 -2(λ+3µ) 3(λ+µ) cos ϕ 2 + cos 3ϕ 2          + • • • + 2 K III (θ) √ 2πρ                 0 0 0 1 2 sin ϕ 2 0 1 2 cos ϕ 2         + ρ R         0 0 0 7 8 sin ϕ 2 -1 2 sin 3ϕ 2 0 5 8 cos ϕ 2 -1 2 cos 3ϕ 2         + • • •         + 2 K III (θ) √ 2πρ ρ R         -2 5 sin ϕ

2. 1 . 1 .

 11 A specific example problem -penny-shaped crack with axisymmetric loading and homogeneous Neumann BCs. As an example problem, consider a penny-shaped crack (Figure 2(c)), ϕ 1 = -π , ω = 2π ( ϕ 2 = π ) in an axisymmetric domain. For the crack in a 2-D cross-section with homogeneous Neumann BCs the following 2-D eigen-pairs are known -they are obtained by solving (15) complemented by BCs (9):

  First three eigen-pairs for a crack with homogeneous Dirichlet BCs.

FIGURE 4 .

 4 FIGURE 4. Solution (left) and error (right) for the axi-symmetric Laplacian with homogeneous Neumann BCs with ρ out = 1/2 , ρ/R = 1/4 and A 1 = 1 -Series up to (ρ/R) 3 . The axis of symmetry is right to the shown domain with the crack from the center of the circle to the right.

FIGURE 5 .

 5 FIGURE 5. The 3-D FE for the torus with a circular crack.

3 .

 3 A specific example problem -hollow cylinder with non-axisymmetric loading and homogeneous Neumann BCs. Consider the circular edge having a solid angle of 3π/2 as the upper corner in Figure2

FIGURE 6 .

 6 FIGURE 6. The solution τ F E on a typical cross-section ( θ = 0 ) of the 3-D torus (Left) and the difference between the analytical and FE solution on same cross-section (Right). The axis of symmetry is right to the shown domain with the crack from the center of the circle to the right.

3. 1 .

 1 Homogeneous boundary conditions. The system of ODEs (70) is complemented by either clamped (homogenous Dirichlet) boundary conditions or traction free (homogeneous Neumann) boundary conditions on ϕ = ϕ 1 and ϕ = ϕ 2 = ϕ 1 + ω . The clamped boundary conditions u(ϕ 1 ) = 0 , u(ϕ 2 ) = 0 imply:

cos ϕ 2 + 1365λ 3 +5271λ 2 µ+6719λµ 2 +2749µ 3 6300sin ϕ 2 -45λ 2 +90λµ+77µ 2 630cos ϕ 2 + 105λ 2 +242λµ+129µ 2 525

 232222

  and the new fact is that-(2 cos ϕ[m 0 ] + [m 01 ]) φ 0,3,1 -cos 2 ϕ[m 0 ] + cos ϕ[m 01 ] + [m 02 ] φ 0,3,0 = 0 allowing the possibility to choose φ 0,3,2 = 0 . Since -cos 2 ϕ[m 0 ] + cos ϕ[m 01 ] + [m 02 ] φ 0,3,1 = 0

  • • • + A 5 (θ)ρ

	1 2	     	-λ+5µ 3(λ+µ) sin ϕ 2 + sin 3ϕ 2 -3λ+7µ 3(λ+µ) cos ϕ 2 + cos 3ϕ 2 0	   +	ρ R	  	0 sin ϕ 2 -λ-3µ 12(λ+µ) sin 3ϕ 2 cos ϕ -123λ 2 +262λµ+107µ 2 180(λ+µ) 2 -15λ 2 +62λµ+79µ 2 180(λ+µ) 2 2 + (λ+5µ) 12(λ+µ) cos 3ϕ 2	  
		+	ρ R	2	  	-483λ 2 +1350λµ+803µ 2 1440(λ+µ) 2 -357λ 2 +1002λµ+581µ 2 1440(λ+µ) 2	sin ϕ 2 + 2625λ 3 +8779λ 2 µ+9251λµ 2 +3161µ 3 6300(λ+µ) 3	sin 3ϕ 2 + 3λ-µ 32(λ+µ) sin 5ϕ 2

Acknowledgements:

The authors gratefully acknowledge an anonymous referee for valuable and constructive comments, leading to improvements in the presentation and context. The first two authors gratefully acknowledge the support of this work by the Israel Science Foundation (grant No. 750/07).

SUMMARY AND CONCLUSIONS

We have formulated the set of recursive ODEs required for the representation of the asymptotic solution in the vicinity of a circular V-notch for the Laplace equation and the elasticity system in a 3-D setting. The asymptotic solution is constructed from 2-D eigen-pairs (a V-notch in a 2-D domain) complemented by two families of shadow functions. Both homogeneous Dirichlet (clamped) or homogeneous Neumann (traction free) boundary conditions were reformulated in terms of eigen-pairs and their shadows, and explicit representation of the series was computed for a penny shaped crack, a circumferential crack and an external crack as well as a 3π/2 reentrant corner.

Having the eigen-pairs and their shadows available, one may construct bench-mark problems for the verification of various numerical schemes. In a future publication the quasidual-singular-function method [START_REF] Costabel | A quasidual function method for extracting edge stress intensity functions[END_REF] for the computation of the Edge-Flux/Stress-Intensity-Functions A k (θ) will be extended to circular singular edges using the series expansion presented herein. The series expansion (shadows) in (90) for the first three terms associated with α 1 = α 2 = α 3 = 0 is explicitly derived herein.

For α 1 = 0 :

] φ 0,2,0 = 0 , the next term φ 0,2,2 can be chosen as zero. Hence we have a solution without any shadow similar to the case of the Laplace operator.

For α 3 = 0 :

We have

Next, notice the disturbing fact is that

]) φ 0,3,0 = 0 and (cos ϕ[t 0 ] + [t 01 ]) φ 0,3,0 = 0, letting the possibility of choosing φ 0,3,1 = 0 . But, with such a choice, since

we would be obliged to calculate a non-zero φ 0,3,2 . Most presumably resulting in an infinite series of shadows.

The special circumstance is that there exists a particular solution with one shadow only if, instead of choosing φ 0,3,1 = 0 , we take


Doru Constantin

The interaction of hybrid nano-particles inserted within surfactant bilayers

Keywords: numbers: 82.70.Dd, 87.16.dt, 61.05.cf Keywords: membrane, inclusions, interaction, X-ray

We determine by small-angle X-ray scattering the structure factor of hydrophobic particles inserted within lamellar surfactant phases for various particle concentrations. The data is then analyzed by numerically solving the Ornstein-Zernicke equation, taking into account both the intra-and inter-layer interactions. We find that particles within the same layer repel each other and that the interaction potential (taken as independent of the concentration) has a contact value of 2.2 k B T and a range of about 10 Å. If the amplitude is allowed to decrease with increasing concentration, the contact value in the dilute limit is about 5 k B T , for a similar range.

I. INTRODUCTION

The elucidation of membrane-mediated interaction between inclusions in the cell membrane (such as integral proteins or membrane-active antibiotic peptides) is of paramount importance for understanding their biological activity. Since the underlying problem is the organization of the membrane -seen as a two-dimensional complex system-it has become clear over the last three decades that the concepts developed in soft matter physics for the understanding of self-assembled systems are operative in this context and that 'simple' models can yield valuable information.

It is therefore not surprising that sustained theoretical efforts attempted to provide a detailed description of these complex systems; they are either continuumelasticity theories [START_REF] Dan | Membrane-induced interactions between inclusions[END_REF][START_REF] Dan | Membraneinduced interactions between inclusions[END_REF][START_REF] Aranda-Espinoza | Interaction between inclusions embedded in membranes[END_REF][START_REF] Bohinc | Interaction between two cylindrical inclusions in a symmetric lipid bilayer[END_REF] or more detailed models taking into account the molecular structure of the lipid bilayer [START_REF] Marčelja | Lipid-mediated protein interaction in membranes[END_REF][START_REF] Sintes | Protein attraction in membranes induced by lipid fluctuations[END_REF][START_REF] Lagüe | Lipid-mediated interactions between intrinsic membrane proteins: A theoretical study based on integral equations[END_REF][START_REF] May | A molecular model for lipid-mediated interaction between proteins in membranes[END_REF] . In the few cases when the models were validated and refined using experimental data, these approaches were often successful. For instance, the lifetime of the gramicidin channel is known to depend on the thickness [START_REF] Elliott | The effects of bilayer thickness and tension on gramicidin singlechannel lifetime[END_REF] and tension [START_REF] Goulian | Gramicidin channel kinetics under tension[END_REF] of the membrane; it was shown that a continuum-elasticity model provides a satisfactory description of this phenomenon [START_REF] Huang | Deformation free energy of bilayer membrane and its effect on gramicidin channel lifetime[END_REF] .

However, this body of theoretical work has not yet been matched by the experimental results, the first of which were obtained by directly measuring the radial distribution function of membrane inclusions using freezefracture electron microscopy (FFEM) [START_REF] Lewis | Bacteriorhodopsin remains dispersed in fluid phospholipid bilayers over a wide range of bilayer thicknesses[END_REF][START_REF] Chen | Temperature-and lightdependent structural changes in rhodopsin-lipid membranes[END_REF][START_REF] James | Lipid-and temperature-dependent structural changes in Acholeplasma laidlawii cell membranes[END_REF][START_REF] Abney | Lateral interactions among membrane proteins: Implications for the organization of gap junctions[END_REF] . These data were compared to liquid state models and could be described by a hard-core model with, in some cases, an additional repulsive or attractive interaction [START_REF] Pearson | Pair distribution functions of bacteriorhodopsin and rhodopsin in model bilayers[END_REF][START_REF] Pearson | Statistical mechanics of lipid membranes, protein correlation functions and lipid ordering[END_REF][START_REF] Braun | Lateral interactions among membrane proteins: Valid estimates based on freezefracture electron microscopy[END_REF] . However, FFEM has not been extensively employed since, certainly due to the inherent experimental difficulties; furthermore, the particle distribution observed in the frozen a) Electronic mail: constantin@lps.u-psud.fr sample is not necessarily identical to that at thermal equilibrium.

A step forward was taken by Huang and collaborators, who showed that the distribution of membrane inclusions within the plane of the layers can be studied using small-angle X-ray and neutron scattering [START_REF] Harroun | Experimental evidence for hydrophobic matching and membrane-mediated interactions in lipid bilayers containing gramicidin[END_REF][START_REF] Yang | Supramolecular structures of peptide assemblies in membranes by neutron off-plane scattering: Method of analysis[END_REF] . First of all, these techniques are perfectly adapted to the typical length scales to be probed. The measurement is averaged over a large number of particles and over long times, without perturbing the system; thus, one has access to the structure factor of the interacting particles. However, in these studies, at most two peptide-to-lipid concentrations P/L were investigated for each system and no values were given for the interaction potential, the data being explained in terms of purely hard-core interactions.

Building upon this work, we recently studied systems with a varying density of inclusions. This is indispensable since, even though each structure factor (taken separately) can be described by a hard-core model, the apparent radius obtained changes with the concentration, signalling the presence of an additional interaction. For instance, we were able to measure the interaction potential of alamethicin pores in dimyristoylphosphatidylcholine (DMPC) bilayers. We showed that, aside from the expected hard-disk repulsion (with a radius corresponding to the geometrical radius of the pore) the pore interaction exhibits a repulsive contribution, with a range of about 3 nm and a contact value of 2.4 k B T 21 . On the other hand, for gramicidin pores in dilauroyl-phosphatidylcholine (DLPC) bilayers, while the interaction is still repulsive, the parameter values are quite different, with a higher contact value and much shorter range. In the latter case, the interaction was shown to decrease with the pore concentration, in agreement with the hydrophobic matching model [START_REF] Constantin | Membrane-mediated repulsion between gramicidin pores[END_REF] . These results are in qualitative agreement with recent theoretical predictions [START_REF] Lagüe | Lipid-mediated interactions between intrinsic membrane proteins: A theoretical study based on integral equations[END_REF][START_REF] Lagüe | Lipid-mediated interactions between intrinsic membrane proteins: Dependence on protein size and lipid composition[END_REF] and can be used as a test for other theoretical and numerical results.

The biological molecules cited above act specifically by insertion within cell membranes, so their affinity for lipid bilayers is assured, which is a significant advantage. However, their use as membrane probes has important shortcomings: their X-ray scattering contrast is low and the positions of the constitutive atoms are not always well-defined (membrane proteins can adopt various conformations as a function of the environment, while the pores formed by antimicrobial peptides often comprise variable numbers of monomers). These features contribute to the difficulty of obtaining high-quality scatering data on membranes with inclusions.

Fortunately, if one is interested in the general physical properties of membranes seen as two-dimensional complex fluids, rather than in the behaviour of a specific active molecule, other inclusions can be used. For instance, we have shown recently that tin oxo-clusters can be inserted within surfactant bilayers and that their interaction potential can be determined [START_REF] Constantin | Repulsion between inorganic particles inserted within surfactant bilayers[END_REF] .

These types of hybrid organic-inorganic particles have the following advantages:

• They are perfectly monodisperse (actually isomolecular) and "rigid" (their atomic configuration is well-defined).

• Their scattering contrast is high (due to the presence of metal atoms).

• Their surface properties can be tailored by changing the nature of the peripheral groups.

The purpose of this article is to develop a full analysis of the system, beyond the initial approach used in Ref. 24. Two main improvements can be noted:

• We were able to determine the complete structure factor, S(q r , q z ), while the initial study was only concerned with the equatorial slice, S(q r , q z = 0). This required new experimental data, for q z = 0.

• We employed a more elaborate statistical model (based on the numerical solution of the Ornstein-Zernicke equation with the Percus-Yevick closure) to relate the structure factor to the interaction potential; our first approach was based on the random-phase approximation (RPA).

We have also checked the reproducibility of the results (see Appendix A), by comparing the old and new data, thus verifying both the reliability of the method and the stability of the samples.

II. MATERIALS AND METHODS

We used as inclusions hybrid nanoparticles consisting of a tin oxide core decorated with butyl chains. The complete formula is:

{(BuSn) 12 O 14 (OH) 6 } 2+ (4- CH 3 C 6 H 4 SO -
3 ) 2 , shortened to BuSn12 in this paper. The synthesis and structural details are given in Ref. 25. The BuSn12 particles were dissolved in ethanol prior to use.

The membranes were composed of dimethyldodecylamine-N-oxide (DDAO), a singlechain zwitterionic surfactant. The DDAO (purchased from Sigma-Aldrich) was first dried in vacuum for 20 h to remove any residual water [START_REF] Kocherbitov | Hydration of dimethyldodecylamine-n-oxide: Enthalpy and entropy driven processes[END_REF] and then dissolved in isopropanol.

The two stock solutions were then mixed to yield the desired BuSn12/DDAO ratio and the mixtures dried in vacuum, yielding a final mass of about 200 mg for each sample. Varying amounts of water were then added so that the mixtures were in the fluid lamellar L α phase [START_REF] Kocherbitov | Hydration of dimethyldodecylamine-n-oxide: Enthalpy and entropy driven processes[END_REF] . For DDAO, the molecular weight is 229.40, the density is 0.84 g/cm [START_REF] Aranda-Espinoza | Interaction between inclusions embedded in membranes[END_REF] and the bilayer thickness is 25 ± 1 Å27,28 , resulting in an area per surfactant molecule A DDAO = 37.8

Å2 . The molecular weight of the BuSn12 is 2866.7, with a density of 1.93 g/cm 3 25 . The (two-dimensional) number density of particles in the plane of the membrane, n, is calculated using the data above and neglecting the possible increase in bilayer surface due to the inclusions.

The lamellar phases were then drawn into flat glass capillaries (VitroCom Inc., Mt. Lks, NJ), 100 µm thick and 2 mm wide by aspiration with a syringe and the capillaries were flame-sealed. Good homeotropic alignment (lamellae parallel to the flat faces of the capillary) was obtained by thermal treatment, using a Mettler FP52 heating stage. The samples were heated up to the isotropic phase (at 130 • C) and then cooled down to the lamellar phase at a rate of 1 • C/min.

We studied the samples by small-angle X-ray scattering (SAXS), on the the bending magnet beamline BM02 (D2AM) [START_REF] Simon | A new small-angle x-ray scattering instrument on the french CRG beamline at the ESRF multiwavelength anomalous scattering/diffraction beamline (D2AM)[END_REF] of the European Synchrotron Radiation Facility (ESRF, Grenoble, France). The photon energy was 11 keV and the sample-to-detector distance about 27 cm, with a scattering vector in the range 0.04 < q < 0.9 Å-1 . The detector is a charge-coupled device Peltier-cooled camera (SCX90-1300, from Princeton Instruments Inc., NJ) with a detector size of 1340 × 1300 pixels. Data preprocessing (dark current subtraction, flat field correction, radial regrouping and normalization) was performed using the bm2img software developed at the beamline. The scattering geometry is discussed in section IV.

III. MODEL AND ANALYSIS

Since the particles are identical in shape, the scattering intensity can be written as the product of a form factor (only depending on the internal constitution of the particles) and a structure factor, which describes the interaction between particles 30 : I(q) = S(q) • |F (q)| 2 , with:

S(q z , q r ) = 1 N N -1 k=1 exp (-iqr k ) 2 ( 1
)
where N is the number of objects and object "0" is taken as the origin of the coordinates.

We expect the form factor of the particles to be dominated by their inorganic core, as the electron density of the butyl chains is similar to that of the dodecyl chains within the bilayers and to that of ethanol. This assumption is confirmed by the intensity at higher scattering vectors (q > 0.5 Å-1 ), which is well described for all samples by the form factor of a sphere |F f (R, q)| 2 , with a radius (used as a free fitting parameter) R = 4.5 ± 0.2 Å, in good agreement with the average radius of the tin oxide "cage" estimated from the crystallographic data.

Dividing the measured intensity by the form factor above yields the structure factor S(q). Standard error propagation then yields the uncertainty, but the resulting values (corresponding to the statistical error) are much smaller than the discrepancy between the fitting model and the experimental data (the goodness-of-fit function χ 2 per data point is much larger than one). Recognizing that this discrepancy is mostly due to systematic errors in data acquisition, background subtraction etc. and also possibly to the inadequacy of the theoretical models, we assign to each data point a fixed uncertainty σ; in the following, we take for definiteness σ = 0.05, which is a rough estimate of the difference between model and data, insofar it yields χ 2 values of unit order.

Determining the structure factor for systems composed of many particles interacting via a known potential is one of the fundamental endeavors of liquid state theory.

Two cases must be considered:

• The situation when the particles are dissolved in a simple solvent will act as a "reference state"; more precisely, it will allow us to check whether the "naked" particles can be described simply by a hard core repulsion or whether an additional term must be considered and what its characteristics are.

• The most interesting situation is of course that of particles inserted within the surfactant bilayers, which will be the main result of our study. Due to the anisotropic nature of the matrix, the analysis must account for the difference between the interaction between particles contained within the same layer and that from one layer to the next.

A. Structure factor of an isotropic solution

In solution, the interaction between particles is taken as isotropic, described by a two-body potential depending only on the distance between their centers: V 3D (r) = V 3D (r). Consequently, the structure factor S(q) is equally isotropic, depending only on the magnitude of the wave vector q.

Using as input the interaction potential, the structure factor was determined numerically by solving iteratively the Ornstein-Zernicke equation with the Percus-Yevick closure condition. The algorithm is based on the method of Lado [START_REF] Lado | Pressure-consistent integral equation for classical fluids: Hard-sphere solutions[END_REF][START_REF] Lado | Equation of state of the hard-disk fluid from approximate integral equations[END_REF] ; see Appendix B for the details. The in- teraction potential is described by a hard core with diameter 2R = 9.4 Å and an additional "soft" component, modelled either as a Gaussian centered at the origin [START_REF] Constantin | Repulsion between inorganic particles inserted within surfactant bilayers[END_REF] :

V (r) = U 0 exp -1 2 r ξ 2
r ≥ 2R or as a decreasing exponential with contact value U c and decay range ξ:

V (r) = U c exp [-(r -2R)/ξ] r ≥ 2R.
B. Structure factor in the lamellar phase

The phase (and hence the structure factor) is now anisotropic. We assume that there is no in-plane ordering (the particles form a 2D liquid in the plane of the bilayer), so that S only depends on the absolute value of the in-plane scattering vector q r = q 2

x + q 2 y and on the scattering vector along the membrane normal q z . For an ideal gas (no interaction) S = 1, while if the particles interact only in the plane of the bilayer and not from one bilayer to the next, S = S(q r ) only. In the general case, the structure factor can be developed as [START_REF] Yang | Supramolecular structures of peptide assemblies in membranes by neutron off-plane scattering: Method of analysis[END_REF] :

S(q r , q z ) = S 0 (q r ) + 2 ∞ m=1 cos(q z dm)S m (q r ) (2)
with each partial structure factor S m describing the interaction between particles situated m layers away:

S m (q r ) = δ 0m + 2πρ ∞ 0 rdrJ 0 (rq r ) [g m (r) -1] (3)
where m ≥ 0 and g m (r) is the (normalized) probability of finding a particle at a distance r in bilayer m, knowing that a particle is present at the origin in bilayer 0. This expansion merely reflects the discrete nature of the stack along the z direction. It is especially convenient since generally the interaction does not extend very far along z. Note also that Equation ( 2) is only valid for the geometrically perfect case when the distance between neighboring layers is rigorously d. In practice, this distance varies due to thermal fluctuations and to frozen-in defects, leading to a smearing of the diffraction pattern at high q z values. We describe this effect phenomenologically by a Lorentzian factor. Finally, we use:

S(q r , q z ) = S 0 (q r ) + 2 cos(q z d) 1 + (q z σ) 2 S 1 (q r ) (4)
where the disorder parameter σ has units of length. For thermal fluctuations, [START_REF] Dan | Membraneinduced interactions between inclusions[END_REF] , where η is the Caillé parameter. As we will see below, the structure factor S 1 (q r ) due to nearest-bilayer interaction has a localized peak. For convenience, we describe it by a Gaussian function:

σ 2 = (z 1 -z 0 ) 2 ≃ η(d/π)
S 1 = A 1 exp - (q -q max ) 2 2 ∆q 2 (5)
The partial structure factors S 0 (q) and S 1 (q) are similar to those describing a two-dimensional binary mixture AB, with the formal identification S 0 = S AA = S BB and S 1 = S AB . As we will see below, these two components are enough to describe the experimental data, so there is no need to go beyond m = 1 (interaction between adjacent bilayers).

The (experimentally determined) functions S 0 (q) and S 1 (q) can be described -using well-established results in liquid state theory-in terms of the interaction potentials V 0 (r) (between particles within the same layer) and V 1 (r) (between adjacent layers). The derivation of the model and the implementation details are given in Appendix B.

IV. RESULTS

A. Interaction in solution

The 3D structure factors of BuSn12 particles in ethanol are shown in Figure 1

B. Interaction between layers

In order to quantify the interaction between the layers (along the smectic director ẑ), one needs access to the complete structure factor S(q r , q z ). This is achieved by using the experimental configuration described in Figure 2, as first discussed in Ref. 20: the incoming X-ray beam (with wave vector k i ) is incident upon the flat capillary at an angle α to its normal n, which coincides with the smectic director ẑ. A point on the 2D detector uniquely defines an outgoing wave vector k o for the scattered signal, and thus a scattering vector given by: q = k o -k i . One can therefore assign to each pixel values for q z = q • ẑ and q r = |qq z ẑ|.

The raw scattering image is then re-gridded to an appropriate region of the reciprocal space, (q r , q z ). For each pixel in the target space, the algorithm identifies the corresponding point in the starting image, and the closest 9 pixels are averaged to yield the final intensity value. It is easily shown that the accessible range in reciprocal space is triangle-shaped, with an angle 2α at the origin. Finally, the resulting image is divided by the form factor of the particle to yield the structure factor.

As an illustration, we show in Figure 2 re-gridded data at four different incidence angles α for the sample with n = 1.97610 -3 Å-2 . A modulation along q z is clearly visible, indicating the presence of the S 1 component and hence of an interaction between layers. The images also display a linear slope in q z probably due to imperfect background subtraction. This artefact is removed prior to further treatment.

Although the physically relevant functions are S 0 (q r ) and S 1 (q r ), the most easily accessible quantity is the "equatorial" cut through reciprocal space S eq = S(q r , q z = 0) = S 0 +2S 1 obtained under normal incidence (α = 0 in Figure 2). We therefore subtract S eq from the re-gridded data (an example is shown in Figure 3, left) and fit it with an appropriately modified version of Equation (4). The best fit for the data in Figure 3 is given in the right image.

The vertical streak on the right in the experimental data is due to an imperfectly subtracted diffuse scatter- ing ring due to a kapton window in the beam (we found that proper background subtraction is quite difficult for tilted samples, i.e. for α = 0). Otherwise, the agreement between data and the fit is good, showing that the model is accurate (in particular, there is no need to include higher-order partial structure factors).

This treatment yields S 1 (q r ), which is then subtracted from the equatorial structure factor S eq (q r ) to yield S 0 (q r ). Figure 4 shows these functions for the most concentrated sample (extracted from the data in Figure 2). In the case of S 1 , the solid line represents the best fit to scattering data taken for α = 60 • (Figure 3), while the shaded area contains the curves obtained for α = 45, 60 and 75 • , providing a rough estimate of the experimental uncertainty.

C. Interaction within the layer

Once the experimental data for the partial structure factors S 0 (q r ) and S 1 (q z ) (or, equivalently, S eq and S 1 ) is available, it can be described in terms of the interaction potentials V 0 (r) and V 1 (r) via well-known methods in the theory of liquids, presented in Appendix B. The main goal of this work is determining the membrane-mediated potential V 0 (r). The interbilayer interaction V 1 (r), although needed for a complete analysis, is inessential and, furthermore, is probably very sensitive to the swelling of the phase. We therefore use a single functional form for it, given by a linear decrease from a maximum U 1 when the particles have the same in-plane position r to 0 when they are a distance ξ 1 away. After extensive tests we found that, for various shapes of the in-plane potential V 0 (r), the best fit for V 1 (r) is described by an amplitude U 1 around 1 k B T and a range ξ 1 of about 25 Å.

In the following, we concentrate on the parameters describing V 0 (r). We considered several functional forms, such as an exponential or a linear decrease, a Gaussian centered at the origin or at the contact and a complementary error function. The best result for each of them is shown in Figure 5, except for the linear decrease which yields much worse fits.

We obtained very good results for the Gaussian centered at contact:

V 0 (r) = U c exp - 1 2 r -2R ξ 2
(solid line in Figure 5; the corresponding structure factors are plotted in Figure 6, along with the experimental data) with fit parameters U c = 2.2 ± 0.2 k B T and ξ = 10.8±0.8 Å while the hard disk radius R is fixed at 4.7 Å. The shaded area corresponds to the fit uncertainty, estimated as in Appendix B 3.

D. Hydrophobic matching

In the treatment above, we considered that the interaction potential is independent of the particle concentration. However, one can infer from a simple elastic theory (and confirm by experimental investigations [START_REF] Constantin | Membrane-mediated repulsion between gramicidin pores[END_REF] ) that the amplitude of the interaction decreases with the concentration of inclusions; briefly, inserting a new particle is easier when the bilayer is already deformed by existing inclusions. Accounting for this effect requires introduc- ing a new fit parameter, the concentration n 0 , for which the interaction decreases significantly. We will consider the simplest model, whereby the membrane deformation is linear (and thus the elastic energy is quadratic) in the concentration: V 0 (n, r) = V 0 (0, r)(1n/n 0 ) 2 for n < n 0 and V 0 (n, r) = 0 for n ≥ n 0 . We take V 1 (r) as independent of the concentration. As before, the in-plane potential is described by a Gaussian centered at contact. The fit is significantly improved, yielding χ 2 = 1.00, for parameters U c = 5.0 ± 0.6 k B T, ξ = 9.3 ± 0.7 Å and n 0 ≃ 5 10 -3 Å-2 (roughly twice the highest experimental value available). We emphasize that the values presented above correspond to the highly diluted limit V 0 (0, r), i.e. to the interaction between two particles in the absence of any other; this is the relevant case for comparison with most theoretical and numerical models.

One can see this decrease of the interaction potential between two particles due to the concentration (i.e. to the presence of other particles) as a crude way of accounting for many-body effects. However, this approximation is not controlled and its accuracy can only be verified by more complete theoretical models or by numerical simulations [START_REF] Harroun | Theoretical analysis of hydrophobic matching and membrane-mediated interactions in lipid bilayers containing gramicidin[END_REF] .

V. DISCUSSION AND CONCLUSION

The main results of this work are the in-plane interaction potential of the BuSn12 particles, V 0 (r) in the lamellar phase and V (r) in solution; they are summarized in . One can see straightaway that the uncertainty is small compared to the amplitude of the curves and to the difference between them, except very close to contact. The qualitative conclusions to be drawn are similar to those already obtained in Reference 24 using a simplified model, although the amplitudes are somewhat different. The interaction potential in solution is shortranged (about 5 Å from contact), but at short distances it becomes stronger than the in-plane repulsion estimated with the constant-amplitude model (4 k B T , as opposed to 2.2 k B T at contact). Notwithstanding the uncertainty and the non-additivity of the potentials, it is tempting to speculate on a possible non-monotonic dependence of the "effective" membrane-mediated interaction, defined by subtracting the interaction in solution from that in the lamellar phase.

0.3 0.2 0.1 q r [Å -1 ] 3.0 2.8 2.6 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 S eq (q r ) = S 0 + 2•S 1 0.5 0.4 0.3 0.2 0.1 0.0 q r [Å

Appendix A: Reproducibility

Since the signals to be measured are fairly weak and unstructured, the question of their reproducibility is of the utmost importance. We therefore measured a few samples twice (about one year apart), in the same experimental conditions (see Section II), and in normal incidence, so that only the equatorial structure factor S eq (q r ) is available. The results are shown in Figure 8 (left). Gray lines: earlier data, used in the analysis of Ref. 24. Black lines: more recent data, measured on the same samples. For clarity, the curves are shifted vertically in steps of 0.25.

Notwithstanding the discrepancy close to the peak position, the agreement is quite good, especially at small angles where S departs significantly from 1. Indeed, this is the range where the interaction potential plays a significant role. For coherence, the comparison is done using the same simplified treatment as in Ref. 24, i. e. neglecting the interaction between layers V 1 (r) and obtaining the Fourier transform V (q r ) of the remaining interaction potential V 0 (r) ≡ V (r) via the RPA approximation; see Ref. 24 for the details.

The two sets of data for V (q r ) are shown in Figure 8 (right). The solid line and the shaded area correspond to the results of Ref. 24 (average value and uncertainty) while the symbols are obtained directly from the more recent values of S eq (q r ) (black lines in the left panel).

Although the new data seems to yield a stronger interaction than the old, the difference is within the estimated uncertainty range. We conclude that the measurement method is reliable and that the samples did not age significantly over a period of one year.

Appendix B: Interaction model

Analytical development

The interaction is treated in the framework of the Ornstein-Zernike relation, with the Percus-Yevick closure approximation. The structure factors S m (q) are computed using the method introduced by Lado 31,32 . An implementation for the case of no interbilayer interaction has already been used [START_REF] Constantin | Membrane-mediated repulsion between gramicidin pores[END_REF] . In the present case, we extend the model to the case of several interacting bilayers. As discussed in the main text, the system can be mapped to a multi-component fluid, where each component corresponds to particles inserted within one bilayer and all particles formally occupy the same twodimensional space. Throughout this analysis, only the two-body interaction is taken into account. FIG. 8. Left: comparison between the equatorial structure factor Seq(qr) determined in two different experimental runs (black and gray curves). A noticeable difference appears for n = 0.584, signaled by the double-headed arrow. Right: comparison between the corresponding potentials V (qr), estimated via the RPA approximation.

This system has been studied for a long time [START_REF] Lebowitz | Exact solution of generalized Percus-Yevick equation for a mixture of hard spheres[END_REF][START_REF] Baxter | Ornstein-Zernike relation and Percus-Yevick approximation for fluid mixtures[END_REF] . In the following, we will use the notations of Ref. 35, except for the interaction potential, denoted here by V (instead of φ) for coherence with previous work. We consider n bilayers (indexed by lowercase Greek indices) with periodic boundary conditions (bilayer n is identical to bilayer 0). The relevant equations are:

The definition of the total correlation function h(r) in terms of the radial distribution function g(r):

h αβ (r) = g αβ (r) -1 (B1)
The Ornstein Zernike relation, defining the direct correlation function c(r):

h αβ (r) = c αβ (r) + γ ρ γ d 2 s c αγ (|s|)h γβ (|r -s|) (B2)
where ρ γ = ρ is the number density in the plane of the bilayers and γ runs over all bilayers. Alternatively, the equation above can be rewritten in reciprocal space, yielding:

h αβ (q) = c αβ (q) + γ ρ γ c αγ (q)h γβ (q) (B3)
The Percus-Yevick closure:

c αβ (r) = [1 + h αβ (r)] 1 -exp V αβ (r) k B T (B4)
with V αβ (r) representing the interaction potential between particles situated in bilayers α and β. For simplicity's sake, we consider that the particles only interact within the same bilayer and with particles in the nearestneighbor bilayers (V αβ (r) ≡ 0 for |α -β| > 1). Furthermore, translation symmetry (an identical environment for all bilayers, i.e. periodic boundary conditions) and mirror symmetry (V αβ ≡ V βα , and similar relations for functions c and h) allow us to define all functions with respect to layer 0, i.e. α = 0 throughout. The equations above must be solved numerically; for convergence reasons, it is much better to work with the indirect correlation function (denoted in Ref. 31 by H(r)), γ 0β (r) = h 0β (r)c 0β (r) instead of the total correlation function.

As an example, for n = 3 bilayers (with periodic boundary conditions, 3 ≡ 0), equations (B3) and (B4) can be written out explicitly as:

γ 00 (q) = ρ [c 00 (γ 00 + c 00 ) + 2c 01 (γ 01 + c 01 )] (B5a) γ 01 (q) = ρ [c 00 (γ 01 + c 01 ) + c 01 (γ 00 + c 00 ) + c 01 (γ 01 + c 01 )] (B5b) c 00 (r) = [1 + γ 00 ][exp(-V 00 /k B T ) -1] (B5c) c 01 (r) = [1 + γ 01 ][exp(-V 01 /k B T ) -1] (B5d)
where we omitted the argument of the functions on the right-hand side. One starts by applying in real space equations (B5c-d) with reasonable initial guesses, and then solving (B5a-b) in reciprocal space. The procedure is iterated until stability is reached. The experimentally relevant structure factors are simply: S 0 (q) = 1 + h 00 (q) and S 1 (q) = h 01 (q). For completeness, the solution of (B5a-b) is:

γ 00 (q) = -c 2 00 (c 00 + c 01 -1/ρ) + 2c 2 01 (c 00 + 1/ρ) (c 00 -1/ρ) 2 + c 01 (c 00 -1/ρ) -2c 2 01 (B6a) γ 01 (q) = c 01 -(c 00 + c 01 )(c 00 -1/ρ) + c 00 /ρ + 2c 2 01 (c 00 -1/ρ) 2 + c 01 (c 00 -1/ρ) -2c 2 01 (B6b)
We checked that, for relevant values of the interaction parameters, the structure factors obtained with n = 5 and n = 7 are almost superposed, and slightly different from those obtained for n = 3 (with periodic boundary conditions in all cases). All results presented above are obtained with n = 7. We also checked that the higher order structure factors are negligible, viz. S 3 (q) ≪ S 2 (q) ≪ S 1 (q), in agreement with the experimental observations.

Software implementation

The routine discussed above is implemented in Igor Pro 6 as an "all-at-once" function, i.e. all points of the output vector are returned at the same time. This approach is particularly useful for such iterative procedures, where, in order to get its value at one point, the entire function needs to be calculated anyway. The Fourier transform was implemented as a matrix operation with pre-computed coefficients over a fixed equidistant grid.

Wherever possible, the wave operations were done using the MatrixOp command.

For each combination of fit parameters, the structure factors S 0 and S 1 are calculated for all densities. Comparison with the experimental data yields the goodnessof-fit function χ 2 , which is minimized using the Optimize operation with the simulated annealing method.

Confidence range

Once a minimum is found, the goodness-of-fit function χ 2 is plotted as a function of the parameters U 0 (or U c ) and ξ.

For the 3D case (interaction in solution), one notices that the minimum is in fact an extended valley, covering a wide range of parameters, and roughly defined as the locus of the points where the integral of the potential is constant. This is understandable, since the fit is mainly sensitive to the low values of the scattering wave vector q, where it is affected by V (q → 0), i.e. by the integral of V (r), rather than by its finer details. We therefore choose a number of points along this valley, with close to minimum χ 2 , and trace the corresponding V (r) instances. The confidence range (shown as a shaded area in Figure 1, bottom) is chosen by manually adjusting the parameters to yield lower and higher envelopes to this sheaf of curves.

In the 2D case (interaction within the layer), the minimum is well localized in the (U c , ξ) plane. However, a standard statistical estimation of the confidence range is irrelevant, since it would yield extremely tight confidence ranges. Indeed, the discrepancy between model and data is mainly due to systematic errors and (presumably) to an inexact functional form for the interaction potential. In order to account for the latter effect we took as an acceptable increase ∆χ 2 the difference between the Gaussian centered at contact (χ 2 = 1.34) and that centered at the origin (χ 2 = 1.40), both plotted in Figure 5. The uncertainty on U c and ξ given in Section IV C is determined based on this numerical value. For comparison, we also show the exponential decay with χ 2 = 1.52. We follow a similar procedure for the variable amplitude model discussed in Section IV D, with ∆χ 2 = 0.1.

  (top) for five different concentrations (lines and symbols). They are fitted with the Gaussian interaction model discussed in Section III A, yielding U 0 = 15 k B T and ξ = 5.1 Å, such that the potential at contact is U c = 3.5 k B T . The fits are shown as solid lines. The corresponding interaction potential is shown in Figure 1 (bottom) as solid line. The shaded area represents the uncertainty (see Appendix B 3 for details). For comparison, we show in the same figure the potential obtained in the RPA approximation for the current data (solid line and symbols) and for that in Ref. 24 (dashed line). The exponential model yields very similar results, in terms of goodness of fit and in potential shape for parameter values U c = 4.5 k B T and ξ = 2.1 Å.

FIG. 2 .

 2 FIG. 2. Left: Sketch of the experimental configuration. Right: Re-gridded scattering data for various values of the incident angle α (indicated within the images).

FIG. 3 .

 3 FIG. 3. Comparison between the experimental data for S(qr, qz) -S(qr, qz = 0) derived from the image taken at α = 60 • in Figure 2 and the model (4).

1 FIG. 4 .

 14 FIG. 4. Partial structure factors for the sample with n = 1.97610 -3 Å-2 , extracted from the data in Figure 2. Symbols: the equatorial structure factor S eq (experimental data). Lines: the structure factors S0 and S1. For S1, the solid line represents the estimate extracted from the configuration with α = 60 • , while the shaded area around the curve represents the uncertainty (see text).

33 FIG. 5 .

 335 FIG.5. Best results for some of the functional forms used for V0(r) and corresponding χ 2 values. For the Gaussian centered at contact, we also plot the uncertainty range obtained as described in Appendix B 3.

1 FIG. 6 .

 16 FIG.6. Best fits (lines) to the experimental data (symbols) with the model described in Section IV C. V0(r) is a Gaussian centered at contact (solid line in Figure5.)

Figure 7

 7 Figure 7. One can see straightaway that the uncertainty is small compared to the amplitude of the curves and to the difference between them, except very close to contact.The qualitative conclusions to be drawn are similar to those already obtained in Reference 24 using a simplified model, although the amplitudes are somewhat different. The interaction potential in solution is shortranged (about 5 Å from contact), but at short distances it becomes stronger than the in-plane repulsion estimated with the constant-amplitude model (4 k B T , as opposed to 2.2 k B T at contact). Notwithstanding the uncertainty and the non-additivity of the potentials, it is tempting to speculate on a possible non-monotonic dependence of the "effective" membrane-mediated interaction, defined by subtracting the interaction in solution from that in the lamellar phase.

FIG. 7 .

 7 FIG. 7.Best results for the in-plane interaction potential V0(r) in the lamellar phase, with constant amplitude (solid line -see Section IV C) and with variable amplitude (dotted line -see Section IV D) as well as for the interaction potential V (r) in solution (dashed line -see Section IV A.) The shaded areas represent the uncertainty, determined as explained in Appendix B 3.

ACKNOWLEDGMENTS

The ESRF is gratefully acknowledged for the award of beam time (experiment 02-01-756) and we thank C. Rochas for competent and enthusiastic support. F. Ribot is acknowledged for providing the BuSn12 particles and for edifying discussions.