Houssem Haddar
email: houssem.haddar@inria.fr

Zixian Jiang
email: zixian.jiang@polytechnique.edu

Armin Lechleiter
email: lechleiter@math.uni-bremen.de

Artificial boundary conditions for axisymmetric eddy current probe problems

Keywords: axisymmetric eddy current equations, artificial boundary conditions, DtN maps for non self-adjoint problems

We study different strategies for the truncation of computational domains in the simulation of eddy current probes of elongated axisymmetric tubes. For axial fictitious boundaries, an exact Dirichlet-to-Neumann (DtN) map is proposed and mathematically justified via the analysis of a nonselfadjoint spectral problem. We rely of results from the spectral theory of compact perturbations of selfadjoint operators. Convergence for truncated DtN map is also analyzed. We validate our theoretical results through numerical experiments for a realistic physical setting inspired by eddy current probes of nuclear reactor core tubes.

Introduction

Industrial context

The present work is motivated by eddy current non-destructive testing of steam generators (SG), see Figure 1, which are critical components in nuclear power plants. Heat produced in a nuclear reactor core is transferred into the primary loop of a steam generator, consisting of tubes in U-shape, and boils coolant water in the secondary loop on the shell side into steam. This steam is then delivered to turbines generating electrical power. Conductive magnetic deposits usually observed on the shell side of the U-tubes could, however, affect the power production and even the structure security. The upper part of the tubes of U-shape is accessible to normal inspection from the top of the steam generator. But it is difficult to reach the lower part of the tubes, which are straight and long, without disassembling the SG. Therefore, a non-destructive examination procedure using signals of eddy currents is applied to detect the presence, the shape and/or the physical nature of deposits on the lower part of the U-tubes.

In eddy current testing (ECT), we use a probe consisting of two coils of wire. Each of these coils is connected to a current generator producing an alternating current and to a voltmeter measuring the voltage change across the coil. Once the probe is introduced in the lower part of some U-tube, the generator coil excited by the current creates a primary electromagnetic field which in turn induces a current flow in the electrically conductive material nearby, such as the tube itself. The presence of conductive magnetic deposits will distort the flow of the eddy currents. They induce a current change in the receiver coil which is measured in terms of impedance and is called the ECT signal (c.f. [START_REF] Auld | Review of advances in quantitative eddy current nondestructive evaluation[END_REF][START_REF] Cartz | Nondestructive Testing: Radiography, Ultrasonics, Liquid Penetrant, Magnetic Particle, Eddy Current[END_REF]).

In order to simulate an eddy current testing experiment, one needs to solve the forward problem for any probe position one wants to incorporate into the measurements. For an iterative inversion method based on the exploitation of this forward problem, the number of required simulation is also proportional to the number of iterations. Given the large number of tubes to be probed, one easily understands the crucial importance of designing a fast (and reliable) numerical simulation of the forward problem. We consider here the eddy current problem under axisymmetric assumption (see for instance [START_REF] Bermúdez | Numerical analysis of a finite-element method for the axisymmetric eddy current model of an induction furnace[END_REF]) and investigate strategies to bound the computational domain. While for the radial direction, truncation with brute model for the boundary condition such as Neumann boundary condition would be sufficient due to the conductivity of the tube and the decay of the solution, in the axial direction this strategy requires some fictitious boundaries far from the sources. We rather propose to compute the exact Dirichlet-to-Neumann (DtN) operator for the region outside the source term and apply it as an exact boundary condition on the fictitious boundaries. This would allow the latter to be as close as needed to the source term. DtN boundary conditions for domain cut-off are widely studied in waveguides and gratings [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF][START_REF] Bonnet-Bendhia | Guided waves by electromagnetic gratings and nonuniqueness examples for the diffraction problem[END_REF][START_REF] Dobson | The time-harmonic Maxwell equations in a doubly periodic structure[END_REF][START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF][START_REF] Hazard | Determination of scattering frequencies for an elastic floating body[END_REF][START_REF] Redon | Non-reflecting boundary conditions for acoustic propagation in ducts with acoustic treatment and mean flow[END_REF]. The main difficulty here is in the justification of the DtN analytical expansion using spectral decomposition of a non-self adjoint operator. We shall rely on results from perturbation theory for the spectrum of compactly perturbed selfadjoint operators [START_REF] Kato | Perturbation theory for linear operators[END_REF]. We also study the error due to truncation in the expression of the DtN operator and relate this to the regularity of the problem parameters. Indeed the latter is important from the computational point of view since this truncation is needed in practice. The DtN expansion relies on some eigenvalues and eigenfunctions that are not known analytically and should be numerically approximated. This may be expensive if a high degree of precision is required. However these calculations can be done off-line and therefore would not affect the speed of solving the problem.

There is a large literature on eddy current problems and without being exhaustive we may refer to [START_REF] Alonso Rodríguez | Eddy current approximation of Maxwell equations[END_REF] for a recent survey on the problem, including an introduction to the eddy current phenomenon, the mathematical justification of the eddy current approximation and different formulations and numerical approaches for the three-dimensional problems. For axisymmetric configurations we refer to the work of [START_REF] Assous | Theoretical tools to solve the axisymmetric Maxwell equations[END_REF] for the study of the theoretic tools for the Maxwell's equations in three dimensions, and to the works of [START_REF] Bermúdez | Numerical analysis of a finite-element method for the axisymmetric eddy current model of an induction furnace[END_REF][START_REF] Chaboudez | Numerical modeling in induction heating for axisymmetric geometries[END_REF] for the discussion of the eddy current problem with bounded conductive components in the meridian half-plane, the numerical analysis and some numerical experiments applied to the induction heating system.

The paper is organized as follows. In Section 2, we briefly recall the eddy current model in the cylindrical coordinate system corresponding to the rotational symmetry with respect to the axis of the tube (see Figure 2) and discuss existence and uniqueness of solution to this problem in its equivalent variational formulation in properly defined weighted function spaces. We then introduce truncations of the domain in the radial-direction by introducing some local boundary conditions (see Section 2.1) and then in the axial-direction by constructing the DtN boundary operator (see Section 3). We validate our analytical theory by several numerical tests that are motivated by ECT experiments as done in practice and present these numerical results in Section 4.

Axisymmetric model

Let us briefly outline the origin of the considered model. We consider the time-harmonic Maxwell's equations for the electric field E and the magnetic field

H curl H + (iωǫ -σ)E = J in R 3 , curl E -iωµH = 0 in R 3 , ( 1
)
where J is the applied electric current density with compact support and satisfies div J = 0, and ω, ǫ, µ, σ respectively denote the frequency, the electrical permittivity, the magnetic permeability and the conductivity. In an axisymmetric (i.e., rotationally invariant) setting, for a vector field a we denote by a m = a r e r + a z e z its meridian component and by a θ = a θ e θ its azimuthal component ). The solution to the first system vanishes if J is axisymmetric. Substituting H m in the second system yields the second-order equation for

E θ = E θ e θ , ∂ ∂r
1 µr ∂ ∂r (rE θ ) + ∂ ∂z 1 µ ∂E θ ∂z + ω 2 (ǫ + iσ/ω)E θ = -iωJ θ in R 2 + , (2)
with R 2 + := {(r, z) : r > 0, z ∈ R}. The applied current J θ has compact support in R 2 + and we denote in the sequel J = J θ . The eddy current approximation corresponds to low frequency regimes and high conductivities: ωǫ ≪ σ. From (2) and the above assumption we get the eddy current model ∂ ∂r

1 µr ∂ ∂r (rE θ ) + ∂ ∂z 1 µ ∂E θ ∂z + iωσE θ = -iωJ in R 2 + , (3)
with a Dirichlet boundary condition at r = 0 due to symmetry E θ | r=0 = 0 and, roughly speaking, a decay condition E θ → 0 as r 2 → ∞ to ensure the finite energy of the electric field. From now on, we denote u = E θ . We introduce operators ∇ := (∂ r , ∂ z ) t and div := ∇• on the half-plane R 2 + and the axis of symmetry Γ 0 := {(r, z) : r = 0, z ∈ R}. Then the axisymmetric eddy current model reads

         -div 1 µr ∇(ru) -iωσu = iωJ in R 2 + , u = 0 on Γ 0 , u → 0 as r 2 → ∞. (4)
We shall assume that µ and σ are in L ∞ (R 2 + ) such that µ ≥ µ 0 > 0 on R 2 + and that σ ≥ 0 and σ = 0 for r ≥ r 0 sufficiently large. For λ > 1 and Ω ⊂ R 2 + , we define the weighted function spaces L 2 1 /2,λ (Ω), H 1

1 /2,λ (Ω) and the norms

L 2 1 /2,λ (Ω) := {v : r 1 /2 (1 + r 2 ) -λ /2 v ∈ L 2 (Ω)}, H 1 1 /2,λ (Ω) := {v ∈ L 2 1 /2,λ (Ω) : r -1 /2 ∇(rv) ∈ L 2 (Ω)}, v L 2 1 /2,λ (Ω) = r (1 + r 2 ) λ v L 2 (Ω) , v 2
H 1 1 /2,λ (Ω) = v 2 L 2 1 /2,λ (Ω) + r -1 /2 ∇(rv) 2 L 2 (Ω)
.

The following lemma gives a Poincaré-type inequality related to functions in H 1 1 /2,λ (R 2 + ). The proof uses classical arguments and is given in Appendix A.1 for the convenience of the reader. Note that for the axisymmetric problems, similar but slightly different weighted spaces have been studied (see for example [START_REF] Mercier | Résolution d'un problème aux limites dans un ouvert axisymétrique par éléments finis en r, z et séries de fourier en θ[END_REF] and [7, Section II.1]).

Lemma 2.1. If v is in H 1 1 /2,λ (R 2 + ) for λ > 1, then lim 0<r→0 v(r, •) L 2 (R) = 0. Moreover, there exists a constant C λ > 0 such that for all v in H 1 1 /2,λ (R 2 + ), v H 1 1 /2,λ (R 2 + ) ≤ C λ r -1 /2 ∇(rv) L 2 (R 2 + ) . (5)
One easily verifies by integration by parts that u in

H 1 1 /2,λ (R 2 + ) is a variational solution of the prob- lem (4) if and only if u satisfies α(u, v) := R 2 + 1 µr ∇(ru) • ∇(rv) dr dz - R 2 + iωσuvr dr dz = R 2 + iωJ vr dr dz ∀v ∈ H 1 1 /2,λ (R 2 + ). (6)
Note that the homogeneous Dirichlet boundary condition on Γ 0 of ( 4) is already included in the space

H 1 1 /2,λ (R 2 + ) by Lemma 2.1. Proposition 2.2. Assume that J ∈ L 2 1 /2,λ (R 2 + ) has compact support. Then the variational problem (6) admits a unique solution u in H 1 1 /2,λ (R 2 + ) for all λ > 1.
Proof. The proof is a direct application of the Lax-Milgram Theorem thanks to (5) which yields the coercivity of the sesquilinear form on the left of (6):

ℜα(v, v) = R 2 + 1 µr |∇(rv)| 2 ≥ 1 µ ∞ C 2 λ v 2 H 1 1 /2,λ (R 2 + ) ,
where C λ is the constant given in (5).

Remark 2.3. The source J has compact support bounded away from Γ 0 in R 2 + in the real problem. We have in particular that J vanishes for r > r 0 and |z| > z 0 , where r 0 > 0 and z 0 > 0 are large enough.

Asymptotic behavior for large r

We are interested here in a more precise evaluation of the decay of the solution u for large argument r. We shall assume in addition to the hypothesis from Proposition 2.2 that the source J and the conductivity σ vanish and that the permeability µ is constant for r ≥ r 0 where r 0 > 0 is some constant. One then gets from (4) that

r 2 ∂ 2 u ∂r 2 + r ∂u ∂r -u + r 2 ∂ 2 u ∂z 2 = 0 for r ≥ r 0 .
We then apply the Fourier transform with respect to the variable z and get

r 2 ∂ 2 u ∂r 2 + r ∂ u ∂r -(1 + 4π 2 ξ 2 r 2 ) u = 0, where u(•, ξ) := +∞ -∞ u(•, z)e -2πiξz dz, ξ ∈ R. (7)
The fundamental solutions of [START_REF] Bernardi | Spectral methods for axisymmetric domains[END_REF] for fixed ξ are the two modified Bessel functions I 1 (2π|ξ|r) and K 1 (2π|ξ|r) when ξ = 0, or the functions r and 1/r when ξ = 0. Since u ∈ L 2 1 /2,λ (R 2 + ) for all λ > 1, the asymptotic behavior for large argument of the modified Bessel functions (see [1, (9.6.1),(9.7.1),(9.7.2)]) implies that u has the following expression for r > r 0 ,

u(r, ξ) =      u(r 0 , ξ) K 1 (2π|ξ|r) K 1 (2π|ξ|r 0 ) ξ = 0, u(r 0 , 0) r 0 r ξ = 0. (8)
Let us also quote that z → u(r, z)

∈ H 1 /2 (R) for r > 0 since u ∈ H 1 1 /2,λ (R 2 + ). Proposition 2.4. The solution u ∈ H 1 1 /2,λ (R 2 + ) to (6) satisfies u(r, •) L 2 (R) ≤ r 0 r u(r 0 , •) L 2 (R) and u(r, •) H 1 /2 (R) ≤ r 0 r u(r 0 , •) H 1 /2 (R) ∀r ≥ r 0 .
Proof. By the Plancherel theorem and the Cauchy-Schwarz inequality, we have

u(r, •) 2 L 2 (R) = u(r, •) 2 L 2 (R) ≤ u(r 0 , •) 2 L 2 (R) K 1 (2π| • |r) K 1 (2π| • |r 0 ) 2 L ∞ (R)
.

We note that K 1 (x) ∼ 1/x as 0 < x → 0 (see [1, (9.6.9)]). Therefore,

g d (ξ; r 0 , r) := K 1 (2πξr) K 1 (2πξr 0 ) → r 0 r as 0 < ξ → 0. (9)
On the other hand, the derivative of g d with respect to ξ is

g ′ d (ξ; r 0 , r) = 2πrK ′ 1 (2πξr)K 1 (2πξr 0 ) -2πr 0 K 1 (2πξr)K ′ 1 (2πξr 0 ) K 2 1 (2πξr 0 ) = 2π [-rK 0 (2πξr)K 1 (2πξr 0 ) + r 0 K 1 (2πξr)K 0 (2πξr 0 )] K 2 1 (2πξr 0 )
, where the last equality follows from the recurrence formulas for Bessel functions (see [1, (9.6.26)]). From the integral representation (see [17, (2.1)] and the references therein)

xK 0 (x) K 1 (x) = 4 π 2 +∞ 0 x 2 x 2 + t 2 t -1 dt J 2 1 (t) + Y 2 1 (t)
x > 0, [START_REF] Brezis | Analyse fonctionnelle[END_REF] and one concludes that the function xK 0 (x)/K 1 (x) is increasing in x > 0, which implies in our case that

r 0 K 0 (2πξr 0 ) K 1 (2πξr 0 ) ≤ rK 0 (2πξr) K 1 (2πξr)
and therefore g ′ d (ξ; r 0 , r) ≤ 0, ξ > 0.

Consequently

g d (ξ; r 0 , r) ≤ g d (0+; r 0 , r) = r 0 r , ∀ξ > 0,
which gives the first inequality of the Proposition using the Plancherel theorem. The second one can be proved with the same arguments.

Proposition 2.5. The trace on vertical lines of the solution u ∈ H

1 1 /2,λ (R 2 + ) to (6) satisfies |∂ r (ru)(r, •)| H -1 /2 (R) := R |ξ| -1 |∂ r (rû)(r, ξ)| 2 dξ 1 /2 < ∞ ∀r ≥ r 0 , (11)
and one has the following bound for r > r 0 sufficiently large

|∂ r (ru)(r, •)| H -1 /2 (R) ≤ C 1 √ r u(r 0 , •) L 2 (R) (12)
for some constant C > 0 independent of r and u.

Proof. We use the following integral representations of modified Bessel's functions [1, (9.6.23)]

K 0 (x) = √ π Γ( 1 2 ) +∞ 1 e -xt dt √ t 2 -1 , K 1 (x) = 1 2 √ π Γ( 3 2 ) +∞ 1 xe -xt t 2 -1 dt = 1 2 √ π Γ( 3 2 ) +∞ 1 e -xt t dt √ t 2 -1 .
They imply in particular

K 0 (a) K 0 (b) ≤ e -(a-b) ∀ a > b > 0 and K 0 (b) K 1 (b) ≤ 2Γ( 3 2 ) Γ( 1 2 ) ∀ b > 0. Therefore K 0 (a) K 1 (b) = K 0 (a) K 0 (b) K 0 (b) K 1 (b) ≤ Ce -(a-b) ∀ a > b > 0. ( 13
)
Differentiating the Fourier representation [START_REF] Bonnet-Bendhia | Propagation dans les guides d'ondes[END_REF] with respect to r then applying the obtained expression to the couple (r, r/2) instead of (r, r 0 ), implies for sufficiently large r,

∂ r (r u)(r, ξ) =    u(r/2, ξ)(-2π|ξ|r) K 0 (2π|ξ|r) K 1 (π|ξ|r) ξ = 0, 0 ξ = 0. (14)
Therefore, using the definition of the norm [START_REF] Cartz | Nondestructive Testing: Radiography, Ultrasonics, Liquid Penetrant, Magnetic Particle, Eddy Current[END_REF], one gets

|∂ r (ru)(r, •)| 2 H -1 /2 (R) = R (2πr) 2 |ξ||û(r/2, ξ)| 2 K 0 (2π|ξ|r) K 1 (π|ξ|r) 2 dξ ≤ u(r/2, •) 2 L 2 (R) g n L ∞ (R) , with g n (ξ) := (2πr) 2 |ξ| K 0 (2π|ξ|r) K 1 (π|ξ|r) 2 .
Using [START_REF] Dautray | Spectral theory and applications, With the collaboration of Michel Artola and Michel Cessenat[END_REF] we have

g n (ξ) ≤ C(2πr) 2 |ξ|e -π|ξ|r ≤ C 4πr e .
Therefore, for r/2 > r 0 ,

|∂ r (ru)(r, •)| H -1 /2 (R) ≤ C 4πr e u(r/2, •) L 2 (R) ≤ 2r 0 C 4π er u(r 0 , •) L 2 (R) ,
where the last inequality is due to Proposition 2.4. This proves the bound (12).

Radial cut-off for eddy current simulations

The decay in radial direction suggests that reasonable accuracy can be obtained by truncating the computational domain at r = r * sufficiently large. In fact, for the application we are interested in, this is also justified by the high conductivity of the tube that would absorb most of the energy delivered by the coil inside the tube (and therefore the value of r * that would be convenient is practice is not too large). We shall analyze in the sequel the error resulting from radial truncation independently from the absorption. It turns out in this case that the boundary conditions that lead to reasonable error estimates are Neumann boundary conditions. The case of Dirichlet boundary conditions lead to slower convergence rates that will be confirmed by our numerical examples. We present in this section only the case of Neumann boundary conditions.

For R ≥ 0 we denote

B R := {(r, z) : 0 < r < R, z ∈ R} and Γ R = {(r, z) : r = R, z ∈ R},
and shall use the short notation

L 2 1 /2 (Ω) := L 2 1 /2,0 (Ω) = {v : v √ r ∈ L 2 (Ω)}, H 1 1 /2 (Ω) := H 1 1 /2,0 (Ω) = {v ∈ L 2 1 /2 (Ω) : r -1 ∇(rv) ∈ L 2 1 /2 (Ω)}.
We remark already that the semi-norm of H 1 1 /2,λ (Ω) written as r -1 /2 ∇(rv) L 2 (Ω) is independent of λ. Then from the proof of Lemma 2.1, this implies that for any v ∈ H 1

1 /2 , we can define its trace on Γ 0 with the same way and the trace vanishes in L 2 (R) norm. Moreover, with H s (R) denoting the usual Sobolev space on R and for sufficiently regular function v defined in a neighborhood of Γ R we set

v H s (ΓR) := v(R, •) H s (R) .
Let r * > 0 be sufficiently large such that the support of the source term J is included in B r * . Then the truncated problem with Neumann boundary conditions on Γ r * consists into seeking

u n ∈ H 1 1 /2 (B r * ) satisfying        -div 1 µr ∇(ru n ) -iωσu n = iωJ in B r * , ∂ ∂r (ru n ) = 0 on Γ r * . (15)
Note again that the function space H 1 1 /2 (B r * ) implicitly accounts for homogeneous boundary condition on Γ 0 . The well-posedness of this problem is guaranteed thanks to the following lemma which will also be useful in quantifying error estimates. The proof of this Lemma is given in Appendix A.2. Lemma 2.6. Let r * > 0. For all v ∈ H 1 1 /2 (B r * ) we have the Poincaré-type inequalities

v L 2 (Br * ) ≤ √ r * r -1 ∇(rv) L 2 1 /2 (Br * ) and v L 2 1 /2 (Br * ) ≤ r * √ 2 r -1 ∇(rv) L 2 1 /2 (Br * ) , (16)
and a trace estimate

|v| H 1 /2 (Γr * ) ≤ r -1 /2 * r -1 ∇(rv) L 2 1 /2 (Br * ) , (17)
where the semi-norm

| • | H 1 /2 is defined by, using the Fourier transform v(r * , •) of v(r * , •), |v| H 1 /2 (Γr * ) := | • | -1 /2 v(r * , •) L 2 (R) . ( 18
)
Proposition 2.7. Assume that the source J ∈ L 2 (R 2 + ) has compact support and let r * > 0 be sufficiently large so that the support of J is included in B r * . Then problem (15) has a unique solution u n ∈ H 1 1 /2 (B r * ). Assume in addition that there exists 0 < r 0 < r * such that J and the conductivity σ vanish and the permeability µ is constant for r > r 0 . Then there exists a constant C that depends only on J, r 0 , µ and σ such that

r -1 ∇(r(u n -u)) L 2 1 /2 (Br * ) ≤ Cr -1 * and u n -u L 2 (Br * ) ≤ Cr -1 /2 * ,
where u is the solution to (6).

Proof. The proof of the first part is similar to the proof of Proposition 2.2 thanks to Lemma 2.6. Let us set

w n := u -u n ∈ H 1 1 /2 (B r * ) such that Br * 1 µr ∇(rw n ) • ∇(rv) -iωσw n vr = Γr * 1 µ ∂ ∂r (ru)v dz ∀v ∈ H 1 1 /2 (B r * ),
where the integral on Γ r * should be understood as a H -1 /2 -H 1 /2 duality pairing. Taking v = wn , we obtain

Br * 1 µr |∇(rw n )| 2 -iωσ|w n | 2 r = Γr * 1 µ ∂ ∂r (ru) wn dz ≤ 1 µ(r * ) ∂ ∂r (ru) H -1 /2 (Γr * ) |w n | H 1 /2 (Γr * ) .
Using [START_REF] Ismail | Complete monotonicity of modified Bessel functions[END_REF] we deduce

1 µ ∞ 1 r ∇(rw n ) 2 L 2 1 /2 (Br * ) ≤ Br * 1 µr |∇(rw n )| 2 -iωσ|w n | 2 r ≤ 1 µ(r * ) ∂ ∂r (ru) H -1 /2 (Γr * ) 1 r * 1 r ∇(rw n ) L 2 1 /2 (Br * )
.

Therefore,

1 r ∇(rw n ) L 2 1 /2 (Br * ) ≤ µ ∞ µ(r * ) 1 r * ∂ ∂r (ru) H -1 /2 (Γr * )
.

The first estimate in the Lemma then follows from Proposition 2.5 and the second one can be deduced using the first inequality in [START_REF] Hazard | Determination of scattering frequencies for an elastic floating body[END_REF].

Remark 2.8. As indicated in the beginning of this section, if one uses Dirichlet boundary conditions on Γ r * then one loses half an order of magnitude for the convergence rate in terms of 1/r * (see Appendix C). This is in fact corroborated by our numerical experiments in Section 4.1.

DtN operator and cut-off in the longitudinal direction

We discuss in this section the truncation of the domain in the longitudinal direction, i.e., the z-direction, whenever a truncation has been applied before in the radial direction. We therefore consider the solution u n of ( 15) and in order to shorten notation we abusively denote this solution by u. Recall that the variational formulation of problem ( 15) is to find

u ∈ H 1 1 /2 (B r * ) such that Br * 1 µ 1 r ∇(ru) • ∇(rv) dr dz - Br * iωσuvr dr dz = Br * iωJ vr dr dz ∀v ∈ H 1 1 /2 (B r * ), (19)
where r * > 0 is as in Proposition 2.7. The idea how to truncate the domain in the z-direction is to explicitly compute the DtN map for the regions above and below the source and inhomogeneities in the coefficients µ and σ using the method of separation of variables. The main difficulty to cope with here is to prove that this is feasible even though the main operator is not selfadjoint. We truncate the domain by two horizontal boundaries Γ ± := {z = ±z * } for some z * > 0 large enough such that the source is compactly supported in

B r * ,z * := {(r, z) ∈ B r * : |z| < z * }.
We then assume in addition that µ and σ only depends on the variable r in the complementary region

B ± r * ,z * := {(r, z) ∈ B r * : z ≷ ±z * }. Since in B ± r * ,z * it holds that -div 1 µr ∇(ru) -iωσu = 0, a solution of the form u(r, z) = ρ(r)ζ(z) has to satisfy 1 ζ d 2 ζ dz 2 = -µ 1 ρ d dr 1 µr d dr (rρ) -iωµσ = ν, (20)
where ν ∈ C is some eigenvalue that we will estimate. For the first equation, we obtain

d 2 ζ dz 2 -νζ = 0,
which has as solutions ζ(z) = c exp(± √ νz) if ν = 0, while for the second equation we are led to consider the eigenvalue problem

S σ µ ρ := -µ d dr 1 µr d dr (rρ) -iωµσρ = νρ in I := {r ∈ R : 0 < r < r * }, ρ(0) = 0, d dr (rρ) r=r * = 0. (21)
We first formally observe from (21) (after multiplication by r ρ and integration by parts) that ℑ(ν) ≤ 0 and ℜ(ν

) > 0. Choosing √ ν such that ℜ √ ν > 0, we get that ζ(z) = c exp(∓ √ νz) on B ± r
* ,z * are the only admissible solutions due to their boundedness at infinity. The only missing point, that would allow the construction of analytic expansions of solutions to [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] restricted to B ± r * ,z * , is to prove that the set of eigenfunctions associated with (21) forms a complete set for the traces on Γ ± of solutions to problem [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF].

Analysis of the non-selfadjoint eigenvalue problem

We consider the spaces

L 2 1 /2 (I) := {φ : φ √ r ∈ L 2 (I)}, H 1 1 /2 (I) := {φ ∈ L 2 1 /2 (I) : 1 r ∂ ∂r (rφ) ∈ L 2 1 /2 (I)}.
For convenience, we shall denote in the sequel by (•, •) µ the weighted scalar product

(a, b) µ := (µ -1 a, b) 1 , where (•, •) 1 is the L 2 1 /2 (I) scalar product. Since µ -1 > 0 is bounded (•, •) µ is equivalent to (•, •) 1 .
We also denote by L 2,µ 1 /2 (I) the L 2 1 /2 (I) space equipped with the weighted scalar product (•, •) µ and by H 1,µ 1 /2 (I) the space defined similarly to H 1 1 /2 (I) but with L 2 1 /2 (I) replaced with L 2,µ 1 /2 (I). We remark that the norm in H 1,µ 1 /2 (I) is also equivalent to the norm in H 1 1 /2 (I). The main results of this section (Lemma 3.2 and Proposition 3.3) will be proved under the assumption that µ is of the following form, which corresponds with the practical problem we are interested in,

µ = µ(r) = µ 0 0 < r < r 1 and r 2 < r < r * , µ t r 1 < r < r 2 , (22)
where 0 < r 1 < r 2 correspond with the two horizontal positions of the tube interfaces and µ 0 and µ t are constants and respectively refer to the magnetic permeabilities of the vacuum and the tube.

Lemma 3.1. The embedding H 1 1 /2 (I) ֒→ L 2 1 /2 (I) is dense and compact. Any φ ∈ H 1 1 /2 (I)
is continuous in the closure of I and satisfies φ = 0 at r = 0. Moreover, the following Poincaré-type inequalities hold,

φ L 2 (I) ≤ √ r * 1 r ∇(rφ) L 2 1 /2 (I) and φ L 2 1 /2 (I) ≤ r * √ 2 1 r ∇(rφ) L 2 1 /2 (I) ∀φ ∈ H 1 1 /2 (I). ( 23
)
Proof. The proof of the compact embedding is a simple application of [START_REF] Brezis | Analyse fonctionnelle[END_REF]Corollaire IV.26]. For the detailed proof, see Appendix A.3. The proof of the property φ(0) = 0 and the inequalities is the same as for Lemma 2.6.

This implies that the embedding H 1,µ 1 /2 (I) ֒→ L 2,µ 1 /2 (I) is dense and compact. Then one can define the unbounded operator

A µ : D(A µ ) ⊂ L 2,µ 1 /2 (I) → L 2,µ 1 /2 (I),
where

D(A µ ) = u ∈ H 1,µ 1 /2 (I) : A µ u = µ d dr 1 µr d dr (ru) ∈ L 2,µ 1 /2 (I), d dr (ru) r=r * = 0 .
For φ ∈ D(A µ ), A µ φ is defined by

(A µ φ, ψ) µ = r * 0 1 µr d dr (rφ) d dr (r ψ) dr ∀ψ ∈ H 1,µ 1 /2 (I). (24)
It is clear from this definition that A µ is closed and selfadjoint and according to Lemma 3.1, it has a compact resolvent. Moreover, the second inequality in Lemma 3.1 and the property of µ show that

(A µ φ, φ) µ ≥ c φ 2 H 1,µ 1 /2 (I) ∀φ ∈ H 1 1 /2 (I)
for some positive constant c independent of φ.

We then deduce (see for instance [13, Section VIII.2, Theorem 7]) that A µ has positive eigenvalues

{λ k } k∈N * with corresponding L 2,µ 1 /2 -complete orthonormal (under (•, •) µ scalar product) eigenprojectors {P k } k∈N * such that 0 < c ≤ λ 1 ≤ λ 2 ≤ • • • ≤ λ k → ∞ (k → ∞),
and ∀φ ∈ P k (H 1,µ 1 /2 (I)) (A µ φ, ψ) µ = λ k (φ, ψ) µ ∀ψ ∈ H 1,µ 1 /2 (I). (25)
Since S σ µ is formally only a compact perturbation of A µ by using the perturbation theory one can relate the spectrum of S σ µ to the spectrum of A µ . We first need to have estimates on the eigenvalues {λ k } h∈N * . For that purpose we introduce the following lemma. Lemma 3.2. Let µ = µ(r) as described in [START_REF] Redon | Non-reflecting boundary conditions for acoustic propagation in ducts with acoustic treatment and mean flow[END_REF]. If µ t /µ 0 is sufficiently close to 1, for instance |µ t /µ 0 -1| < 18 , then the eigenvalues {λ k } k∈N * of A µ are simple and grow like O(k 2 ) as k → ∞. Moreover, the difference λ k+1 -λ k → +∞ as k → ∞.

The proof of this Lemma with involved calculations is given in Appendix B. Let us note that in practice, for steam generator tubes, the condition |µ t /µ 0 -1| < 1 8 is satisfied. In fact µ t /µ 0 -1 ≃ 0.01. We also note that all subsequent analysis remains true if we assume that µ ∈ L ∞ (I) with 0 < µ inf ≤ µ ≤ µ sup < +∞ and the corresponding selfadjoint operator A µ has simple eigenvalues {λ k } k∈N * growing up as O(k 2 ) and λ k+1 -λ k → ∞ as k → ∞. Now let us consider the operator S σ µ = A µ + M σ µ defined in [START_REF] Pedreira | A method for computing guided waves in integrated optics. II. Numerical approximation and error analysis[END_REF]. Since the multiplication operator

M σ µ : φ → -iωµσφ ∀φ ∈ L 2,µ 1 /2 (I),
is bounded on L 2,µ 1 /2 (I), the theory for perturbed selfadjoint operators [18, Theorem V-4.15a and Remark V-4.16a] implies: Proposition 3.3. Under the assumptions of Lemma 3.2, the unbounded operator S σ µ : L 2,µ 1 /2 (I) → L 2,µ 1 /2 (I) is closed with compact resolvent and its eigenvalues and eigenprojectors can be indexed as {ν 0j , ν k } and {Q 0j , Q k } respectively, where j = 1, . . . , m < ∞ and k = n + 1, n + 2, . . . with n ≥ m ≥ 0 such that the following results hold:

1. the sequence |ν k -λ k | is bounded as k → ∞.
2. there exists a bounded operator W on L 2,µ 1 /2 (I) with bounded inverse W -1 such that

Q 0 := m j=1 Q 0j = W -1   k≤n P k   W and Q k = W -1 P k W for k > n. (26)
Moreover, {Q 0j , Q k } is a complete family in the sense that

m j=1 Q 0j + k>n Q k = 1. ( 27
)
We note that although the eigenvalues ν k may be simple (this is for instance the case for piece-wise constant coefficients) we can still have m < n in the previous proposition. This would correspond with the presence of Jordan's blocs for eigenvalues decompositions of matrices. However, this case has not been observed numerically (as mentioned in the numerical section).

Spectral decomposition of the DtN operator

We are now in position to provide explicit expression for the DtN operator that will be used to truncate the domain in the z-direction. We first need to specify the space of traces on Γ ± of functions in H 1 1 /2 (B r * ). From the definition of the spectral decomposition of A µ we immediately deduce that for φ ∈ H 1,µ 1 /2 (I)

φ 2 L 2,µ 1 /2 (I) = k≥1 P k φ 2 L 2,µ 1 /2 (I)
and φ 2

H 1,µ 1 /2 (I) = k≥1 (1 + λ k ) P k φ 2 L 2,µ 1 /2 (I) . For θ ∈ [0, 1], we define H θ,µ 1 /2 (I) as the θ interpolation space [H 1,µ 1 /2 (I), L 2,µ 1 /2 (I)] θ (see [19, Définition 2.1] for interpolation spaces) with norm φ 2 H θ,µ 1 /2 (I) := k≥1 (1 + λ k ) θ P k φ 2 L 2,µ 1 /2 (I) (28)
and define H θ,µ Γ ± with I using the obvious isometry. Let v be a regular function of B r * . We denote the trace mapping by

γ ± : v → v| Γ± .
Theorem 3.4. The trace mapping γ ± can be extended to a continuous and surjective mapping from

H 1 1 /2 (B r * ,z * ) onto H 1 /2,µ 1 /2 (Γ ± ) and from H 1 1 /2 (B ± r * ,z * ) onto H 1 /2,µ 1 /2 (Γ ± ). Proof.
Obviously we have the equivalent definition

H 1 1 /2 (B r * ,z * ) = H 1,µ 1 /2 (B r * ,z * ) = v : v ∈ L 2 ((-z * , z * ); H 1,µ 1 /2 (I)), ∂v ∂z ∈ L 2 ((-z * , z * ); L 2,µ 1 /2 (I)) ,
with the same norm. Therefore the trace mapping properties for H Let us also mention the following result that will be useful later

Lemma 3.5. Let θ ∈ [0, 1]. The norm φ H θ,µ 1 /2 := ∞ k=1 (1 + λ k ) θ P k φ 2 L 2,µ 1 /2 (I)
defines an equivalent norm on H θ,µ 1 /2 (I). Proof. From interpolation theory, it is sufficient to prove the result for θ = 0 and θ = 1. The statement is obviously true for θ = 1, while for θ = 0 one needs only to observe that

H 0,µ 1 /2 (I) = L 2,µ 1 /2 (I). Let φ ± in H 1 /2,µ
1 /2 (Γ ± ) and denote by µ ± and σ ± the restrictions of µ and σ to B ± r * ,z * . Thanks to Theorem 3.4 one can uniquely define

u ± ∈ H 1 1 /2 (B ± r * ,z * ) solution of              -div 1 µ ± r ∇(ru ± ) -iωσ ± u ± = 0 in B ± r * ,z * , 1 µ ± ∂ ∂r (ru ± ) = 0 on Γ r * , u ± = φ ± on Γ ± . (29)
The construction of u ± can be done for instance by using some continuous lifting linear operators R ± :

H 1 /2,µ 1 /2 (Γ ± ) → H 1 1 /2 (B ± r * ,z * ) such that γ ± R ± (φ) = φ (these operators exist according to Theorem 3.4). The H 1 1 /2 (B ± r * ,z * ) norm of u ± indeed continuously depends on the H 1 /2,µ
1 /2 (Γ ± ) norm of the boundary data φ ± (respectively). Definition 3.6. We define the DtN operators

T ± : H 1 /2,µ 1 /2 (Γ ± ) → H 1 /2,µ 1 /2 (Γ ± ) * by T ± φ ± , ψ ± = B ± r * ,z * 1 µ ± r ∇(ru ± )∇(rR ± ψ ± ) dr dz + B ± r * ,z * div 1 µr ∇(ru ± ) R ± ψ ± r dr dz for all ψ ± ∈ H 1 /2,µ 1 /2 (Γ ± ), where u ± ∈ H 1 (B ± r * ,z * )
is the unique solution of problem (29) and where •, • denotes the

H 1 /2,µ * 1 /2 -H 1 /2,µ 1 /2
duality product that coincides with (•, •) µ for L 2,µ 1 /2 functions. Indeed

T ± φ ± , ψ ± = B ± r * ,z * 1 µ ± r ∇(ru ± )∇(rR ± ψ ± ) dr dz - B ± r * ,z * iωσ ± uR ± ψ ± r dr dz (30)
and therefore, from the definition of R ± and the continuity property for the solutions u ± ,

T ± φ ± , ψ ± ≤ C φ ± H 1 /2,µ 1 /2 (Γ±) ψ ± H 1 /2,µ 1 /2 (Γ±)
for some constant C independent from φ and ψ. This proves that T ± : H

1 /2,µ 1 /2 (Γ ± ) → H 1 /2,µ
1 /2 (Γ ± ) * are well-defined and are continuous. We remark that for sufficiently regular u ± , we have (using Green's formula)

T ± φ ± = ∓ 1 µ ± ∂u± ∂z | Γ± . (31)
Since γ ± -γ ± R ± γ ± = 0, we also observe, using Green's formula (and a density argument) that

T ± φ ± , γ ± v = B ± r * ,z * 1 µ ± r ∇(ru ± )∇(rv) dr dz + B ± r * ,z * div 1 µr ∇(ru ± ) vr dr dz ± for all v ∈ H 1 1 /2 (B ± r * ,z *
). Therefore we also have

T ± φ ± , γ ± v = B ± r * ,z * 1 µ ± r ∇(ru ± )∇(rv) dr dz - B ± r * ,z * iωσ ± u ± vr dr dz for all v ∈ H 1 1 /2 (B ±
+ T + γ + u, γ + v + T -γ -u, γ -v = Br * ,z * iωJ vr dr dz ∀v ∈ H 1 1 /2 (B r * ,z * ). (32)
We immediately get the following equivalence result.

Proposition 3.7. A function u ∈ H 1 1 /2 (B r * ) is solution of (19) if and only if u| Br * ,z * ∈ H 1 1 /2 (B r * ,z *
) and is solution of (32) and u = u ± on B ± r * ,z * where u ± ∈ H 1 1 /2 (B ± r * ,z * ) are solution of (29) with φ ± = γ ± (u| Br * ,z * ).

Formulation (32) is the one that we would like to use in practice. Proposition 3.7 and the wellposedness of [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] show that (32) is also well-posed. To be numerically effective one needs explicit expressions for T ± . We shall use for that purpose Proposition 3.3. We are then led to consider the spectral decompositions of S σ + µ + and S σ - µ -that correspond to the one in Proposition 3.3 for (µ, σ) = (µ + , σ + ) and (µ, σ) = (µ -, σ -) respectively. Since the treatment of both cases is the same and in order to simplify the notation we shall use the same notation for the spectral decomposition of S σ + µ + and S σ - µ -. For φ ± ∈ H 1 /2,µ 1 /2 (Γ ± ) we have the spectral decomposition

φ ± = m j=1 Q 0j (φ ± ) + k>n Q k (φ ± ).
By definition of Q 0j and Q k the functions u ± defined on B ± r * ,z * by

u ± (r, z) = m j=1 Q 0j (φ ± )(r) exp(∓ √ ν 0j (z ∓ z * )) + k>n Q k (φ ± )(r) exp(∓ √ ν k (z ∓ z * )) in B ± r * ,z * , (33)
(the square root is determined as the one with positive real part) formally satisfy (29). In order to rigorously prove this, one only needs to verify that this function is in

H 1 1 /2 (B ± r * ,z * ). Since the eigenfunctions Q 0j (φ ± ) and Q k (φ ± ) are in H 1,µ 1 /2 (I), one easily checks that u ± N (r, z) := m j=1 Q 0j (φ ± )(r) exp(∓ √ ν 0j (z ∓ z * )) + N k=n+1 Q k (φ ± )(r) exp(∓ √ ν k (z ∓ z * )) in B ± r * ,z * , (34)
is in H 1 1 /2 (B ± r * ,z * ) and verifies (29) with boundary data on Γ ± equal

φ ± N = m j=1 Q 0j (φ ± ) + N k=n+1 Q k (φ ± ).
We then can conclude using the following lemma.

Lemma 3.8. Let φ ∈ H 1 /2,µ 1 /2 (I) and set for N > n, φ N = m j=1 Q 0j (φ) + N k=n+1 Q k (φ). Then, φ N -φ H 1 /2,µ 1 /2 (I) → 0 as N → ∞.
The proof of this Lemma is itself a straightforward consequence of the following result since, using the notation of the Lemma below,

φ N -φ 2 H 1 /2,µ 1 /2 (I) = k>N (1 + |ν k |) 1 /2 P k W φ 2 L 2,µ 1 /2 (I) .
Lemma 3.9. Let θ ∈ [0, 1] and let ν * ∈ R. The norm defined by

φ 2 H θ,µ 1 /2 (I) := ∞ k=1 (1 + |ν k |) θ P k W φ 2 L 2,µ 1 /2 (I)
where ν k , k > n, are the eigenvalues of S σ µ as defined in Proposition 3.3 and ν k = ν * for k ≤ n, defines an equivalent norm on H θ,µ 1 /2 (I).

Proof. We first observe that thanks to Lemma 3.5, the result is obvious for θ = 0 since H 0,µ 1 /2 (I) = L 2,µ 1 /2 (I) and

φ H 0,µ 1 /2 (I) = W φ 2 H 0,µ 1 /2 (I)
where W : L 2,µ 1 /2 (I) → L 2,µ 1 /2 (I) is an isomorphism. Using interpolation theory one then only needs to prove the result for θ = 1. The case of θ = 1 will also be proved using interpolation theory since, using again Lemma 3.5 and the definition of A µ , we have

H 1,µ 1 /2 (I) = [D(A µ ), L 2,µ 1 /2 (I)]1 /2 . Therefore it is sufficient to prove that φ 2 D(Aµ) = φ 2 L 2,µ 1 /2 (I) + A µ φ 2 L 2,µ 1 /2 (I) = ∞ k=1 (1 + λ 2 k ) P k φ 2 L 2,µ 1 /2 (I)
is equivalent to

φ 2 H 2,µ 1 /2 (I) = ∞ k=1 (1 + |ν k | 2 ) P k W φ 2 L 2,µ 1 /2 (I) .
Using the identity P k (W S σ µ φ) = ν k P k W φ for k > n, we observe that

φ 2 H 2,µ 1 /2 (I) = φ 2 H 0,µ 1 /2 (I) + S σ µ (φ -Q 0 φ) 2 H 0,µ 1 /2 (I) + |ν * | 2 Q 0 φ 2 H 0,µ 1 /2 (I) . (35)
Since

S σ µ (Q 0 φ) 2 L 2,µ 1 /2 (I) ≤ C Q 0 φ 2 L 2,µ 1 /2 (I) with C = sup{|ν 0j | 2 , j = 1, m}, then S σ µ (I -Q 0 ) = A µ + M 0 , where M 0 := M σ µ -S σ µ Q 0 is a bounded operator on L 2,µ
1 /2 (I). Therefore, with C denoting a constant independent of φ but whose value may change from a line to another, and using the first part of the proof,

φ 2 D(Aµ) ≤ 2 S σ µ (φ -Q 0 φ) 2 L 2,µ 1 /2 (I) + 2 M 0 φ 2 L 2,µ 1 /2 (I) + φ 2 L 2,µ 1 /2 (I) ≤ C φ 2 H 2,µ 1 /2 (I)
and

φ 2 H 2,µ 1 /2 (I) ≤ 2 A µ φ 2 H 0,µ 1 /2 (I) + 2 M 0 φ 2 H 0,µ 1 /2 (I) + (1 + |ν * | 2 ) φ 2 H 0,µ 1 /2 (I) ≤ C φ 2 D(Aµ) ,
which proves the desired equivalence of norms and concludes the proof.

The expression of u ± N in B ± r * ,z * (34) yields,

∓ ∂u ± N ∂z Γ± = m j=1 √ ν 0j Q 0j (φ ± ) + n<k≤N √ ν k Q k (φ ± ).
Therefore, using (31) and letting N → ∞ we obtain (explicitly specifying in the notation the dependence on Γ ± on the spectral decomposition)

T ± (φ ± ) = m ± j=1 ν ± 0j Q ± 0j φ ± + k>n ± ν ± k Q ± k φ ± . ( 36
)

On the analysis of spectral error truncation

For numerical simulations the spectral representation of operators T ± should be truncated. We shall give here some estimates on the error due to this truncation. For N > n ± , we define the projectors

Q ± N := m ± j=1 Q ± 0j + n ± <k≤N Q ± k ,
and the truncated DtN operators

T ± N := T ± Q ± N = m ± j=1 ν ± 0j Q ± 0j + n ± <k≤N ν ± k Q ± k . (37)
According to Lemma 3.9,

Q ± N : H 1 /2,µ 1 /2 (Γ ± ) → H 1 /2,µ 1 /2 (Γ ± ) * are continuous and therefore T ± N : H 1 /2,µ 1 /2 (Γ ± ) → H 1 /2,µ 1 /2 (Γ ± ) * are also continuous. We are interested in considering u N ∈ H 1 1 /2 (B r * ,z * ) solving Br * ,z * 1 µ ∇(ru N ) • 1 r ∇(rv) dr dz - Br * ,z * iωσu N vr dr dz + T + N γ + u N , Q + N γ + v + T - N γ -u N , Q - N γ -v = Br * ,z * iωJ vr dr dz ∀v ∈ H 1 1 /2 (B r * ,z * ). (38)
This variational problem is well-posed for all N as indicated in the following. Using the Riesz representation theorem, we introduce A 0 :

H 1 1 /2 (B r * ,z * ) → H 1 1 /2 (B r * ,z * ), A : H 1 1 /2 (B r * ,z * ) → H 1 1 /2 (B r * ,z * ) and A N : H 1 1 /2 (B r * ,z * ) → H 1 1 /2 (B r * ,z * ) defined by (A 0 w, v) H 1 1 /2 (Br * ,z * ) = Br * ,z * 1 µ ∇(rw) • 1 r ∇(rv) -iωσwvr dr dz, (Aw, v) H 1 1 /2 (Br * ,z * ) = (A 0 w, v) H 1 1 /2 (Br * ,z * ) + T + γ + w, γ + v + T -γ -w, γ -v , and
(A N w, v) H 1 1 /2 (Br * ,z * ) = (A 0 w, v) H 1 1 /2 (Br * ,z * ) + T + Q + N γ + w, Q + N γ + v + T -Q - N γ -w, Q - N γ -v ,
for all v ∈ H 1 1 /2 (B r * ,z * ), respectively. We recall that the operator A 0 is coercive, and more precisely

ℜ(A 0 w, w) H 1 1 /2 (Br * ,z * ) ≥ a 0 w 2 H 1 1 /2 (Br * ,z * )
for some positive constant a 0 independent of w. We observe from (30) that ℜ T ± φ ± , φ ± ≥ 0, and therefore

ℜ(A N w, w) H 1 1 /2 (Br * ,z * ) ≥ a 0 w 2 H 1 1 /2 (Br * ,z * ) .
This means in particular, thanks to the Lax-Milgram theorem that A N is bijective and also

A -1 N ≤ 1/a 0 .
Consequently problem (38) has a unique solution u N ∈ H 1 1 /2 (B r * ,z * ) that continuously depends on J with a modulus of continuity independent of N .

From the continuity of T ± we easily obtain

(A -A N )w H 1 1 /2 (Br * ,z * ) ≤ C γ + w -Q + N γ + w H 1 /2,µ 1 /2 (Γ+) + γ -w -Q - N γ -w H 1 /2,µ 1 /2 (Γ-) (39)
for some constant C independent of N and w ∈ H 1 1 /2 (B r * ,z * ). Therefore, using Lemma 3.9,

lim N →∞ (A -A N )w H 1 1 /2 (Br * ,z * ) = 0 ∀w ∈ H 1 1 /2 (B r * ,z * ).
With u ∈ H 1 1 /2 (B r * ,z * ) denoting the solution of (32), we observe that

Au = A N u N . Therefore, u -u N = A -1 N (A N u -Au).

This proves in particular that

u -u N H 1 1 /2 (Br * ,z * ) ≤ 1/a 0 (A -A N )u H 1 1 /2 (Br * ,z * ) → 0 as N → ∞. (40)
We can summarize these results in the following proposition Proposition 3.10. Under the same assumptions as in Proposition 2.7, the variational problem (38) has a unique solution

u N ∈ H 1 1 /2 (B r * ,z * ). Moreover, if u ∈ H 1 1 /2 (B r * ,z * ) is the solution of (32), then u -u N H 1 1 /2 (Br * ,z * ) → 0 as N → ∞.
We shall now give some indication on the rate of convergence under some additional regularity assumptions on the source term J and the coefficients µ and σ. Obviously, from (39) and ( 40)

u -u N H 1 1 /2 (Br * ,z * ) ≤ C/a 0 γ + u -Q + N γ + u H 1 /2,µ 1 /2 (Γ+) + γ -u -Q - N γ -u H 1 /2,µ 1 /2 (Γ+) . (41)
Therefore the speed of convergence will depend on the regularity of γ ± u. Considering problem [START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF] satisfied by u in the unbounded domain B r * and differentiating the equations with respect to z (i.e. considering the equation satisfied (u(r, z + ∆z) -u(r, z))/∆z, then letting ∆z → 0) one easily observes from the well-posedness of problem [START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF] that if in addition

∂ m J ∂z m ∈ L 2 (B r * ), ∂ m σ ∂z m ∈ L ∞ (B r * ) and ∂ m µ -1 ∂z m ∈ L ∞ (B r * ),
for some integer m ≥ 0, then

∂ m u ∂z m ∈ H 1 1 /2 (B r * ).
Consequently, if this holds with m = 2, then the first equation in (15) yields

γ ± S σ ± µ ± u = - 1 µ ± γ ± ∂ 2 u ∂z 2 ∈ L 2,µ 1 /2 (Γ ± ).
With the help of Lemma 3.9, we can then estimate,

γ ± u -Q N ± γ ± u 2 H 1 /2,µ 1 /2 (Γ+) ≤ C 1 ∞ k=N +1 (1 + |ν ± k |) 1 /2 P ± k W ± φ 2 L 2,µ 1 /2 (I) ≤ C 2 (1 + |ν ± N |) 3 /2 ) ∞ k=N +1 (1 + |ν ± k |) 2 P ± k W ± φ 2 L 2,µ 1 /2 (I) ≤ C 3 (1 + |ν ± N |) 3 /2 ) φ 2 H 0,µ 1 /2 (Γ±) + S σ ± µ ± u 2 H 0,µ 1 /2 (Γ±)
where the constants C 1 , C 2 and C 3 are independent of N and where we used in the second inequality the fact that |ν k | → ∞ as k → ∞. According to Proposition 3.3 and Lemma 3.2,

|ν N | ≥ CN 2
for some constant C > 0 independent from N . From the discussion above we then can deduce the following theorem.

Theorem 3.11. Under the assumptions as in Proposition 2.7 and the additional assumptions that

∂ m J ∂z m ∈ L 2 (B r * ), ∂ m σ ∂z m ∈ L ∞ (B r * ) and ∂ m µ -1 ∂z m ∈ L ∞ (B r *
), for m = 0, 1, 2, there exists a constant C that only depends on J, µ, σ, r * and z * such that

u -u N H 1 1 /2 (Br * ,z * ) ≤ C N 3/2
, where u ∈ H 1 1 /2 (B r * ,z * ) and u N ∈ H 1 1 /2 (B r * ,z * ) are the respective solutions of (32) and (38). We end this section with a two remarks. The first one is on exponential convergence rates and was brought to our attention by one of the referees. Remark 3.12. If we change the variational formulation into

u N ∈ H 1 1 /2 (B r * ,z * ) solving Br * ,z * 1 µ ∇(ru N ) • 1 r ∇(rv) dr dz - Br * ,z * iωσu N vr dr dz + T + N γ + u N , γ + v + T - N γ -u N , γ -v = Br * ,z * iωJ vr dr dz ∀v ∈ H 1 1 /2 (B r * ,z * ), (42)
then, following the approach in [START_REF] Pedreira | A method for computing guided waves in integrated optics. II. Numerical approximation and error analysis[END_REF] for waveguides (see also [START_REF] Bonnet-Bendhia | Propagation dans les guides d'ondes[END_REF]Chapter 3]), one can prove convergence at exponential rates by exploiting the exponential decay of the (truncated) analytic expression (33) in the z direction. However, for this formulation, one is able to prove well-posedness of (42) only for N sufficiently large.

The second one is related to the case of Dirichlet boundary conditions.

λ k = ρ 1 (t(k)) r * 2
where ρ 1 (t) is a positive zero of the cylinder function C 1 (•; t).

Numerical validation

We recall the two-dimensional geometric representation of the eddy current testing procedure in the Orz plan from Figure 2 or, more precisely, Figure 3a. In the following examples, the two coils involved are represented by two rectangles with 0.67mm in length (radial direction) and 2mm in height (longitudinal direction). They are located 7.83mm away from the z-axis and have a distance of 0.5mm between them. The SG tube measures 9.84mm in radius for the interior interface and 11.11mm for the exterior interface. We assume some deposit with a rectangular shape on the shell side of the tube with 2mm in length and 6cm in height. The probe coils and the deposit are placed symmetrically with regard to the r-axis. The permeability and conductivity of different materials are given in Table 1. The background permeability µ 0 is the permeability of vacuum.

vacuum tube deposit permeability µ v = µ 0 µ t = 1.01µ 0 µ d = 10µ 0 conductivity (in S • m -1 ) σ v 0 σ t = 1 × 10 3 σ d = 1 × 10 4
Table 1: Values of the physical parameters for the numerical examples.

To approximate solutions to the original eddy current problem ( 6

Truncation error in the radial direction

Next we truncate the computational domain much closer to the tube at r = r * , see Figure 3b (Br * ,Z ) and in the norm • L 2 (Br * ,Z ) with higher rates than the relative error of the Dirichlet problem (61) does. (Figure 4a). This observation precisely corresponds to our theoretical results (see Remark 2.8). Note that it is reasonable that the convergence rates observed in numerical tests are better than theoretical ones which give just a lower bound of estimate. The advantage of truncating the computational domain in the radial direction using a Neumann instead of a Dirichlet boundary condition is clearly confirmed by these examples.

In eddy current testing, one is interested in particular in measurements of impedances, which only depend on the solution inside the deposit domain Ω d . To this end, we also compare the relative error due to truncation in the radial direction on Ω d . From Figures 4c and4d, at a truncation position r * = 50mm (log 10 r * ≈ -1.3), the relative error issued from the Neumann problem in the semi-norm of H 1

1 /2 (Ω d ) is less than 0.1% and that in the norm of L 2 (Ω d ) are less than 1%. Therefore we conclude that simulations computed in a domain truncated at r = r * = 50mm using Neumann boundary conditions are sufficiently precise for iterative reconstruction algorithms, since the noise level in the measurements would most probably be higher that the numerical error. Concerning the finite element space on B r * ,Z , this truncation reduced the degrees of freedom in our experiments to about 16000.

-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 log10 truncation position log10 relative error H 1/2 1 semi-norm Dirichlet Neumann (a) | • | H 1 1/2 (B r * ,Z ) -semi-norm -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 log10 truncation position log10 relative error L 2 norm Dirichlet Neumann (b) • L 2 (B r * ,Z ) -norm -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -8 -7 -6 -5 -4 -3 -2 -
(c) | • | H 1 1/2 (Ω d ) -norm -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -8 -7 -6 -5 -4 -3 -2 -
.c., r -1 in H 1 1 /2 semi-norm and r -1 /2 in L 2 norm.

Error introduced by the DtN maps

In the following, we denote by u exact a reference solution for the eddy current problem computed on the truncated infinite band To build the DtN maps, we first discretize the interval I = (0, r * ) (that has the same length as Γ ± ) into 5000 equi-length segments and use an eigenvalue solver (more precisely, the function eigs in Matlab) to compute the first eigenvalues {ν ± 0j , ν ± k } and the corresponding eigenvectors corresponding to the physical parameters µ ± and σ ± . In our numerical example, we observe that the eigenvalues are simple. This means in particular that no Jordan block is present for this set of eigenvalues since otherwise the function eigs would return an eigenvalue with a multiplicity equal to the dimension of the Jordan bloc sub-space (and also return linearly dependent eigenvectors). The eigenprojections {Q ± 0j , Q ± k } are approximated by matrices which we obtain by interpolating the complex conjugate of these eigenvectors with the boundary elements on Γ ± (141 elements on each boundary), using the bi-orthogonality between the eigenvectors of the operator and those of its adjoint. The sum of these matrices weighted with the square root of the corresponding eigenvalues yields numerical approximations of the truncated DtN maps introduced in (37).

Figure 5 illustrates the relative errors of u DtN using different truncation parameters N for the DtN operator T ± N , see (37), with respect to u exact in the

• L 2 1 /2 (Br * ,z * ) -norm and the |•| H 1 1 /2
(Br * ,z * ) -semi-norms. The relative error decreases as the truncation order N increases before saturating at about N = 16. For N = 20, the errors are sufficiently small. Let us note that in the case σ t = 0, i.e. the operator S σ µ = A µ is selfadjoint, using N = 1 would be sufficient to achieve the same accuracy. (Br * ,z * ) -semi-norm compared to the reference solution. Again, one clearly observes that using the DtN maps for the horizontal truncation introduces a reasonably small error compared to the reference solution u exact while truncating using Dirichlet-or Neumann boundary conditions on horizontal boundaries close to the coils and the deposit yields unacceptable errors. In particular, merely pre-computed DtN maps can ensure fast simulations of non-destructive eddy current measurements when many forward problems need to be solved. As mentioned in the introduction, such fast simulations are crucial for, e.g., iterative solution methods for the inversion of these measurements.

A Some properties of the weighted spaces

We recall the definition of the following spaces and corresponding norms. Let Ω ⊂ R 2 + := {(r, z) : r > 0, z ∈ R} an open set. For λ > 1,

L 2 1 /2,λ (Ω) := {v : r 1 /2 (1 + r 2 ) -λ /2 v ∈ L 2 (Ω)}, H 1 1 /2,λ (Ω) := {v ∈ L 2 1 /2,λ (Ω) : r -1 /2 ∇(rv) ∈ L 2 (Ω)}, v L 2 1 /2,λ (Ω) = r (1 + r 2 ) λ v L 2 (Ω) , v 2
H 1 1 /2,λ (Ω) = v 2 L 2 1 /2,λ (Ω) + r -1 /2 ∇(rv) 2 L 2 (Ω)
.

For λ = 0, we define

L 2 1 /2 (Ω) := L 2 1 /2,0 (Ω) = {v : v √ r ∈ L 2 (Ω)}, H 1 1 /2 (Ω) := H 1 1 /2,0 (Ω) = {v ∈ L 2 1 /2 (Ω) : r -1 ∇(rv) ∈ L 2 1 /2 (Ω)}.
We shall also use the short notation

|v| 2 H 1 1 /2 (Ω) = r -1 /2 ∇(rv) 2 L 2 (Ω)
.

For r * > 0 and an interval I = {r ∈ R : 0 < r < r * } we define

L 2 1 /2 (I) := {φ : φ √ r ∈ L 2 (I)}, H 1 1 /2 (I) := {φ ∈ L 2 1 /2 (I) : r -1 ∂ r (rφ) ∈ L 2 1 /2 (I)}.
A.1 Proof of Lemma 2.1

Proof. First of all, we prove that the trace on Γ 0 as the limit of traces on Γ ǫ exists and its L 2 -norm vanishes. Given 0 < ǫ < r * , we set B ǫ r * := {(r, z) ∈ B r * : r ≥ ǫ} and I ǫ := {r ∈ R : ǫ < r < r * }.

One easily observes that

L 2 1 /2,λ (B ǫ r * ) = L 2 (B ǫ r * ) and H 1 1 /2,λ (B ǫ r * ) = H 1 (B ǫ r * ) ⊂ L 2 (R, H 1 ((ǫ, r * ))
. Since H 1 ((ǫ, r * )) ⊂ C((ǫ, r * )), for 0 < ǫ < r < r ′ < r * and for almost all z ∈ R, we can write for v ∈ H 1 1 /2,λ (B ǫ r * ),

|r ′ v(r ′ , z) -rv(r, z)| = r ′ r ∂ ∂s (sv(s, z)) ds ≤ |r ′ -r| 1 /2 √ r * |v(•, z)| H 1 1 /2 (I ǫ ) , R |r ′ v(r ′ , z) -rv(r, z)| 2 dz ≤ |r ′ -r|r * R |v(•, z)| 2 H 1 1 /2 (I ǫ ) dz ≤ |r ′ -r|r * |v| 2 H 1 1 /2 (B ǫ r * ) .
Thus, for r n → 0 (n → ∞), {r n v(r n , •)} n∈N * is a Cauchy sequence in L 2 (R). Since L 2 (R) is complete, the sequence converges and we denote the L 2 (R)-norm of its limit by l ≥ 0. Now we will show that l = 0.

If not, due to the continuity of rv on r for almost all z, one should have

∃δ > 0 ∀0 < r < δ R |rv(r, z)| 2 dz ≥ l 2 2 .
For 0 < ǫ < δ < r * , with Fubini's theorem,

v 2 L 2 1 /2,λ (B ǫ r * ) ≥ v 2 L 2 1 /2,λ (B ǫ δ ) = R δ ǫ 1 r(1 + r 2 ) 1 /2,λ |rv(r, z)| 2 dr dz ≥ l 2 2 1 (1 + δ 2 ) λ δ ǫ 1 r dr ǫ→0 ---→ ∞,
which contradicts the fact that v ∈ L 2 1 /2,λ (R 2 + ). Hence we have proved that l = 0. Therefore, for almost all z ∈ R and v ∈ H

1 1 /2,λ (B ǫ r * ) ⊂ L 2 (R, H 1 ((ǫ, r * )), |v(r, z)| 2 = 1 r 2 |rv| 2 = 1 r 2 r 0 ∂ ∂s (sv(s, z)) ds 2 ≤ r 0 1 √ s ∂ ∂s (sv(s, z)) 2 ds ∀r > 0. (43)
For any r > 0, we integrate the above inequality for z over R and get

v(r, •) 2 L 2 (R) = R |v(r, z)| 2 dz ≤ R r 0 1 √ s ∂ ∂s (sv(s, z)) 2 ds dz. (44)
We conclude that

lim r→0 v(r, •) 2 L 2 (R) = 0. (45)
Now using (44),

R 2 + r (1 + r 2 ) λ |v| 2 dr dz ≤ ∞ -∞ ∞ 0 r (1 + r 2 ) λ dr ∞ 0 1 √ r ∂ ∂r (rv(r, z)) 2 dr dz = ∞ 0 r (1 + r 2 ) λ dr |v| 2 H 1 1 /2 (R 2 + ) .
Therefore, the inequality is proved by setting (I) . We integrate the above inequality for r over I and obtain

C = 1 + ∞ 0 r (1 + r 2 ) λ dr.
v(•, z) 2 L 2 (I) ≤ r * 0 dr|v(•, z)| 2 H 1 1 /2 (I) = r * |v(•, z)| 2 H 1 1 /2 (I) , v(•, z) 2 L 2 1 /2 (I) ≤ r * 0 r dr|v(•, z)| 2 H 1 1 /2 (I) = r 2 * 2 |v(•, z)| 2 H 1 1 /2 (I) .
Therefore, integration over z ∈ R yields

v 2 L 2 (Br * ) = v(•, z) L 2 (I) 2 L 2 (R) ≤ r * |v(•, z)| 2 H 1 1 /2 (I) 2 L 2 (R) ≤ r * |v| 2 H 1 1 /2 (Br * ) , v 2 L 2 1 /2 (Br * ) = v(•, z) L 2 1 /2 (I) 2 L 2 (R) ≤ r 2 * 2 |v(•, z)| 2 H 1 1 /2 (I) 2 L 2 (R) ≤ r 2 * 2 |v| 2 H 1 1 /2 (Br * ) .
By setting C p = √ r * in the first inequality and C p = r * / √ 2 in the second inequality, the Poincaré-type inequalities [START_REF] Hazard | Determination of scattering frequencies for an elastic floating body[END_REF] are proved. We now prove [START_REF] Ismail | Complete monotonicity of modified Bessel functions[END_REF]. We recall that

|v| 2 H 1 /2 (Γr * ) = R |ξ||v(r * , ξ)| 2 dξ = 1 r * R |ξ|r * |v(r * , ξ)| 2 dξ .
For ǫ > 0,

r * |v(r * , ξ)| 2 -ǫ|v(ǫ, ξ)| 2 = r * ǫ ∂ r (r|v| 2 ) dr = r * ǫ 2ℜ(v∂ r (rv)) -|v| 2 dr ≤ 2 r * ǫ |v||∂ r (rv)| dr.
Multiplying the inequality with |ξ|, integrating for |ξ| < M , letting first ǫ → 0 (use (45) to observe that the ǫ terms in the left hand side goes to zero) and then letting M → ∞, yields

|v| 2 H 1 /2 (Γr * ) ≤ 2 r * R |ξ| r * 0 |v||∂ r (rv)| dr. Therefore |v| 2 H 1 /2 (Γr * ) 2 r * R r * 0 |ξ| r 1 /2 v r -1 /2 ∂ r (rv) dr dξ ≤ 1 r * R r * 0 |ξ| 2 r 1 /2 v 2 + r -1 /2 ∂ r (rv) 2 dr dξ = 1 r * |v| 2 H 1 1 /2 (Br * ) .
This proves [START_REF] Ismail | Complete monotonicity of modified Bessel functions[END_REF].

A.3 Proof of Lemma 3.1

Proof. We suppose B is a unit ball in H 1 1 /2 (I). To prove the compactness of B in L 2 1 /2 (I), it is sufficient to show that B := {φ(•) √ • : φ ∈ B} is compact in L 2 (I). We use [START_REF] Brezis | Analyse fonctionnelle[END_REF]Corollaire IV.26]. We suppose for arbitrary η > 0 small enough, ω ⊂]η, r * -η[ is strictly included in I, written as ω ⊂⊂ I. We note τ h the translation operator:

(τ h φ)(r) = φ(r + h).
First of all, we shall show

∀h ∈ R with |h| < η and ∀ψ = φ(•) √ • ∈ B, τ h ψ -ψ L 2 (ω) h→0 ---→ 0.
For r ∈ ω, we have |h| < η < r. For h > 0,

|ψ(r + h) -ψ(r)| 2 = |φ(r + h) √ r + h -φ(r) √ r| 2 ≤ |φ(r + h)(r + h) -φ(r)r| 2 1 r + |φ(r + h)| 2 ( √ r + h - √ r) 2 r + h r = r+h r d ds (sφ(s)) ds 2 1 r + |φ(r + h)| 2 h √ r + h + √ r 2 r + h r ≤ h r r+h r d ds (sφ(s)) 2 ds + |φ(r + h)| 2 h 2 √ r 2 2r r ≤ 2h|φ| 2 H 1 1 /2 (I) + h 2 |φ(r + h)| 2 ,
thus by Lemma 2.6 and the fact that φ ∈ B

τ h ψ -ψ 2 L 2 (ω) ≤ 2hr * |φ| 2 H 1 1 /2 (I) + h 2 φ 2 L 2 (I) ≤ 2hr * |φ| 2 H 1 1 /2 (I) + h 2 r * |φ| 2 H 1 1 /2 (I) ≤ 5h 2 r * h→0 ---→ 0.
For h < 0, we note always h > 0 but we calculate

|ψ(r -h) -ψ(r)| 2 = |φ(r -h) √ r -h -φ(r) √ r| 2 ≤ |φ(r)r -φ(r -h)(r -h)| 2 1 r + |φ(r -h)| 2 ( √ r - √ r -h) 2 r -h r = r r-h d ds (sφ(s)) ds 2 1 r + |φ(r -h)| 2 h √ r + √ r -h 2 r -h r ≤ h r r r-h d ds (sφ(s)) 2 ds + |φ(r -h)| 2 h 2 √ r -h 2 r -h r ≤ h|φ| 2 H 1 1 /2 (I) + h 4 |φ(r -h)| 2
again by Lemma 2.6 we have

τ -h ψ -ψ 2 L 2 (ω) ≤ hr * |φ| 2 H 1 1 /2 (I) + h 4 r * |φ| 2 H 1 1 /2 (I) ≤ 5h 4 r * h→0 ---→ 0.
It remains to prove that

∀ǫ > 0 ∃ω ⊂⊂ I such that φ(•) √ • L 2 (I\ω) < ǫ ∀φ ∈ B.
If we take ω = (η, r * -η), then by Lemma 2.6

φ(•) √ • 2 L 2 (I\ω) = φ 2 L 2 1 /2 ((0,η)) + φ 2 L 2 1 /2 ((r * -η,r * )) ≤ η 2 2 |φ| 2 H 1 1 /2 ((0,η)) + η 2 2 |φ| 2 H 1 1 /2 ((r * -η,r * )) ≤ η 2 2 .
By setting η small enough we obtain the result.

[φ] = µ -1 d dr (rφ) = 0 r = r 1 , r 2 , (46a)
φ(0) = 0 and d dr (rφ) r * = 0. (46b)
Here [•] is the jump operator. By setting ζ = √ λ, we rewrite (46a) in the form of a Bessel's equation

r 2 φ ′′ + rφ ′ + (ζ 2 r 2 -1)φ = 0.
Using the first condition in (46c), we get φ in the following form up to a constant factor

φ(r) = φ(ζr) :=      J 1 (ζr) 0 < r < r 1 , aC 1 (ζr; s) r 1 < r < r 2 , bC 1 (ζr; t) r 2 < r < r * ,
where J 1 (•) is the first kind Bessel function of order 1, C 1 (•; s) is the cylinder function of order 1. In general, a cylinder function is defined as

C ν (x; s) := J ν (x) cos(πs) + Y ν (x) sin(πs), (47)
where Y ν is the second kind Bessel function of order ν. Due to the second boundary condition in (46c),

ζ = ρ 0 (t) r *
where ρ 0 (t) is a zero of cylinder function C 0 (•; t). Finally, a, b ∈ R and s, t > 0 are the four constants to be determined by the four jump conditions (46b) at r = r 1 and r = r 2 . To simplify the notation, we will note for the moment

x 1 = ρ0(t) r * r 1 , x 2 = ρ0(t) r * r 2 .
Considering the recurrence relations of cylinder functions [1, (9.1.27)], the jumps conditions (46b) write

J 1 (x 1 ) = aC 1 (x 1 ; s), 1 µ 0 J 0 (x 1 ) = a µ t C 1 (x 1 ; s), (48a)
aC 1 (x 2 ; s) = bC 1 (x 2 ; t), a µ t C 0 (x 2 ; s) = b µ 0 C 0 (x 2 ; t). (48b)
Using the definition of cylinder functions (47) and the Wronskians of Bessel functions [1, (9.1.16)], one obtains from (48a)

tan(πs) = - ( µt µ0 -1)J 0 (x 1 )J 1 (x 1 ) 2 πx1 -( µt µ0 -1)J 0 (x 1 )Y 1 (x 1 ) , (49)
and from (48b) tan(πt) = -

( µt µ0 -1)J 0 (x 2 )J 1 (x 2 ) + ( µt µ0 -1)J 0 (x 2 )Y 1 (x 2 ) -2 πx2 tan(πs) 2 πx2 + ( µt µ0 -1)Y 0 (x 2 )J 1 (x 2 ) + ( µt µ0 -1)Y 0 (x 2 )Y 1 (x 2 ) tan(πs) . ( 50
)
If we denote ( µt µ0 -1) by δ, we get from (49) -( 50)

tan(πt) = - δ J 0 (x 2 )J 1 (x 2 ) 2 πx1 -δJ 0 (x 1 )Y 1 (x 1 ) + 2 πx2 -δJ 0 (x 2 )Y 1 (x 2 ) J 0 (x 1 )J 1 (x 1 ) 4 π 2 x1x2 + δ 2 πx1 Y 0 (x 2 )J 1 (x 2 ) -2 πx2 J 0 (x 1 )Y 1 (x 1 ) -δ 2 Y 0 (x 2 )Y 1 (x 2 )J 0 (x 1 )J 1 (x 1 )
.

Hence, if λ = ζ 2 = ( ρ0(t) r * ) 2 is an eigenvalue of A µ , then t is a zero of the following function

ϑ(t) := tan(πt) + χ( ρ0(t) r * ) with χ(ζ; δ) := δ J 0 (ζr 2 )J 1 (ζr 2 ) 2 πζr1 -δJ 0 (ζr 1 )Y 1 (ζr 1 ) + 2 πζr2 -δJ 0 (ζr 2 )Y 1 (ζr 2 ) J 0 (ζr 1 )J 1 (ζr 1 ) 4 π 2 ζr1ζr2 + δ 2 πζr1 Y 0 (ζr 2 )J 1 (ζr 2 ) -2 πζr2 J 0 (ζr 1 )Y 1 (ζr 1 ) -δ 2 Y 0 (ζr 2 )Y 1 (ζr 2 )J 0 (ζr 1 )J 1 (ζr 1 ) . (52)
From [1, Section 9.5], if we fix ρ 0 (0) = 0, then we can write ρ 0 (t) a continuous and increasing function of the continuous variable t and we have the McMahon's expansions for large t (see [1, (9.5.12)]):

ρ 0 (t) = β + 1 8β + O(β -3 ) (t → +∞), (53)
where Therefore, we get the existence of a zero of ϑ(t) on (k - By (54) and the fact that r 1 , r 2 < r * , we have | d dt χ( ρ0(t) r * )| < π. So we conclude that there is a unique zero of ϑ(t) on (k -1 2 , k + 1 2 ). Without losing the generality, we denote by λ k = ( ρ0(t k ) r * ) 2 the k th eigenvalue of A µ with t k ∈ (k -1 2 , k + 1 2 ). For |δ| < min{δ 1 , δ 2 } = 1 8 , we have from the above discussion that there exists K > 0 such that

β = t - 1 4 π, B ′ 0 (ζ) = - 8 π 2 ζ 3 r 1 r 2 , B ′ 1 (ζ) = - 4 πζ 2 r 1 Y 0 (ζr 2 )J 1 (ζr 2 ) + 2r 2 πζr 1 Y 0 (ζr 2 )J 0 (ζr 2 ) -Y 1 (ζr 2 )J 1 (ζr 2 ) + 4 πζ 2 r 2 J 0 (
B 0 (ζ) = 4 π 2 ζ 2 r 1 r 2 , B ′ 0 (ζ) = - 8 π 2 ζ 3 r 1 r 2 , B 1 (ζ) ∼ 2 π 2 ζ 2 r 1 r 2 , B ′ 1 (ζ) ∼ - 8 π 2 ζ 3 r 1 r 2 , B 2 (ζ) ∼ 2r 1 π 2 r 2 ln(ζ), B ′ 2 (ζ) ∼
t k ∈ k -1 4 , k + 1 4 ∀k > K.
Using McMahon's expansion (53), one gets that λ k grows like O(k 2 ),

λ k ∼ π 2 r 2 * k 2 as k → ∞.
Moreover,

λ k+1 -λ k = 1 r 2 * ρ 2 0 (t k+1 ) -ρ 2 0 (t k ) > 1 r 2 * k + 1 - 1 4 - 1 4 2 π 2 -k + 1 4 - 1 4 2 π 2 + O(k -1 ) = 1 r 2 * 1 2 2k - 1 2 π 2 + O(k -1 ) → +∞ as k → ∞.
We conclude that the difference λ k+1 -λ k → +∞ as k → ∞.

C Dirichlet boundary conditions for radial domain truncation

Definition C.1. We build a lifting operator R r * : H 1 /2 (R) → H 1 1 /2 (B r0 ) such that its Fourier transform satisfies (R r * φ)(r, ξ) = I 1 (2π|ξ|r) I 1 (2π|ξ|r * ) φ(ξ),

where I 1 is the modified Bessel function.

We verify easily that (R r * φ)| r=0 = 0, (R r * φ)| r=r * = φ and -div r -1 ∇(rR r * φ) = 0 in B r * . By multiplying the previous divergence-free condition with rR r * φ and integrating by parts, we get r -1 ∇(rR r * φ) 2

Considering the asymptotic behavior of I 0 , I 1 with big argument, we have

C(r * ) ∼ O( √ r * ) r * → ∞. ( 60
)
So the lifting operator R r * grows with a rate of (r * ) 1 /2 when r * tends to infinity.

Remark C.2. Lemma 2.6 and (60) show that the lifting R r * is "minimal" in the sense that its norm grows with the least rate, i.e. as (r * ) 1 /2 when r * tends to infinity.

Figure 1 :

 1 Figure 1: Sketch of a steam generator.

2 Figure 2 :

 22 Figure 2: Three-and two-dimensional geometric representations of a steam generator tube covered with deposits and a probe consisting of two coils.

Remark 3 . 13 .

 313 The results and proofs of this section apply also to the case where the Neumann boundary conditions on r = r * are replaced with Dirichlet boundary conditions. The only modification would be the replacement of H 1 1 /2 (B) by H 1 1 /2,0 (B) := {u ∈ H 1 1 /2 (B); u = 0 on r = r * } where B stands for B r * or B r * ,z * . The eigenvalues λ k are in this case

2 +

 2 ) on the unbounded domain R by numerical simulations, we use a domain B R,Z with very large truncation parameters R = 300mm, Z = 100mm and we set Neumann conditions on these boundaries. These values of R and Z are large enough to ensure that the corresponding reference solution is close enough to the true solution to be able to study the (non-)convergence of the different domain truncations presented above. All numerical examples are done using the open-source finite element software FreeFem++. The computation of the reference solution uses a mesh that is adaptively refined with respect to this solution with a maximum edge size h max = 2mm as well as P1 finite elements on the mesh. The degrees of freedom of the finite element space are about 49900 for B R,Z with R = 300mm and Z = 100mm.

Figure 3 :

 3 Figure 3: Truncation in the radial and longitudinal directions.

  by setting Dirichlet or Neumann boundary conditions on Γ r * = {r = r * }. Using the same physical parameters as above and setting again Z = 100mm to approximate solutions to the truncated problem on B r * on the domain B r * ,Z ; again, the value for Z gave sufficient numerical accuracy in our tests. In Figure 4 we show the numerical results corresponding to the convergence results of Proposition 2.7 and Proposition C.3. As r * increases, the relative error of the the Neumann problem (15) tend to zero in the semi-norm |•| H 1 1 /2

Figure 4 :

 4 Figure 4: Relative errors with different truncation positions in the radial direction and different boundary conditions. Theoretical convergence rates: with Dirichlet b.c., r -1 /2 in H 1 1 /2 semi-norm and no convergence in L 2 norm; with Neumann b.c., r -1 in H 1 1 /2 semi-norm and r -1 /2 in L 2 norm.

  B r * with r = r * = 50mm using Neumann boundary conditions on Γ r * , compare Figure 3b. To compute u exact numerically, we resolve the problem in a domain bounded B r * ,Z with Z = 100mm, as explained above. Then B r * ,Z is truncated into the bounded domain B r * ,z * with Γ ± = {0 < r < r * , z = ±z * = ±5mm}, compare Figure 3c. The degrees of freedom of the P1 finite element space reduced by this truncation to about 3500 elements. We set different boundary conditions -Dirichlet, Neumann or DtN boundary conditions -on the top and bottom boundaries Γ ± and solve the corresponding variational problems again using the finite element software package FreeFem++. The solutions are denoted by u Dirichlet , u Neumann and u DtN in the following.

Figure 5 :

 5 Figure 5: Relative errors for eddy current simulations using DtN maps with different truncation orders N .

Figure 6 2 1 / 2 (

 622 Figure 6 illustrates real and imaginary parts of solutions for the three different horizontal truncation techniques (Dirichlet, Neumann, and DtN) we investigated above. It shows in particular that the truncation using DtN maps constructed with the first 20 eigenvalues and eigenprojections approaches the most the exact model. Moreover, Table 2 indicates the relative errors of the truncated models in the • L 2 1 /2 (Br * ,z * ) -norm and the | • | H 1 1 /2

  ℜ(uDirichlet) (d) ℑ(uDirichlet) (e) ℜ(uNeumann) (f) ℑ(uNeumann) (g) ℜ(uDtN) (h) ℑ(uDtN)

Figure 6 :

 6 Figure 6: Real and imaginary parts of u fields of the truncated model with different boundary conditions. DtN maps of truncation order N = 20.

A. 2 2 H 1 1 / 2

 2212 Proof of Lemma 2.6 Proof. Let r ∈ I = (0, r * ). From inequality (43) in the proof of Lemma 2.1, one has |v(r, z)| 2 ≤ |v(•, z)|

χπ 2 ζ 2 r 1 r 2 (r 1 + r 2 π 2 ζ 2 r 1 r 2 r 1 π 4 ζ 4 r 2 1 r 2 2 (r 1 + r 2 ), C 2 (ζ) ∼ 8 π 4 ζ 4 r 2 1 r 2 2 D 2 C 3 (ζ) ∼ 4 π 4 ζ 4 r 2 1 r 2 2 D 3 r 2 2 D 4 + 3δ 2 1 4 -8δ 1 -δ 2 1 < 1

 21212122232324111 ′ (ζ; δ) = 0. (56)Using the principal asymptotic forms of Bessel functions for large arguments[1, (9.2.1)-(9.2.2)], one has as ζ → +∞,A 1 (ζ) ∼ -2 π 2 ζ 2 r 1 r 2 cos(2ζr 1 ) + cos(2ζr 2 ) , A 2 (ζ) ∼ -1 π 2 ζ 2 r 1 r 2 cos(2ζr 1 ) + cos(2ζr 2 ) + sin(2ζ(r 1 + r 2 )) , B 1 (ζ) ∼ 2 π 2 ζ 2 r 1 r 2 2 + sin(2ζr 1 ) -sin(2ζr 2 ) , B 2 (ζ) ∼ 1 π 2 ζ 2 r 1 r 2 cos(2ζr 1 ) cos(2ζr 2 ), sin(2ζr 1 ) + r 2 sin(2ζr 2 ) -(r 1 + r 2 ) cos(2ζ(r 1 + r 2 )) , B ′ 1 (ζ) ∼ 4 π 2 ζ 2 r 1 r 2 r 1 cos(2ζr 1 ) -r 2 cos(2ζr 2 ) , B ′ 2 (ζ) ∼ -2 π 2 ζ 2 r 1 r 2 r 1 sin(2ζr 1 ) cos(2ζr 2 ) + r 2 cos(2ζr 1 ) sin(2ζr 2 ) .Therefore, as ζ → +∞,C 1 (ζ) ∼ -16 (ζ; r 1 , r 2 ), (ζ; r 1 , r 2 ), C 4 (ζ) ∼ 2 π 4 ζ 4 r 2 1 (ζ; r 1 , r 2 ).whereD 2 (ζ; r 1 , r 2 ) := r 1 sin(2ζr 1 ) + r 2 sin(2ζr 2 ) -(r 1 + r 2 ) cos(2ζ(r 1 + r 2 )) -(r 1 + r 2 ) 2 + sin(2ζr 1 ) -sin(2ζr 2 ) + cos(2ζr 1 ) + cos(2ζr 2 ) r 1 cos(2ζr 1 ) -r 2 cos(2ζr 2 ) , D 3 (ζ; r 1 , r 2 ) := r 1 sin(2ζr 1 ) + r 2 sin(2ζr 2 ) -(r 1 + r 2 ) cos(2ζ(r 1 + r 2 )) 2 + sin(2ζr 1 ) -sin(2ζr 2 )+ cos(2ζr 1 ) + cos(2ζr 2 ) + sin(2ζ(r 1 + r 2 )) r 1 cos(2ζr 1 ) -r 2 cos(2ζr 2 )-(r 1 + r 2 ) cos(2ζr 1 ) cos(2ζr 2 )-cos(2ζr 1 ) + cos(2ζr 2 ) r 1 sin(2ζr 1 ) cos(2ζr 2 ) + r 2 cos(2ζr 1 ) sin(2ζr 2 ) , D 4 (ζ; r 1 , r 2 ) := r 1 sin(2ζr 1 ) + r 2 sin(2ζr 2 ) -(r 1 + r 2 ) cos(2ζ(r 1 + r 2 )) cos(2ζr 1 ) cos(2ζr 2 )-cos(2ζr 1 ) + cos(2ζr 2 ) + sin(2ζ(r 1 + r 2 )) r 1 sin(2ζr 1 ) cos(2ζr 2 ) + r 2 cos(2ζr 1 ) sin(2ζr 2 ) .Therefore,D 2 (•; r 1 , r 2 ) ∞ ≤ 8(r 1 + r 2 ), D 3 (•; r 1 , r 2 ) ∞ ≤ 10(r 1 + r 2 ), D 4 (•; r 1 , r 2 ) ∞ ≤ 6(r 1 + r 2 ).We have|χ(ζ; δ)| ≤ 4|δ| + 3|δ| 2 4 -8|δ| -|δ| 2 as ζ → ∞, (57) |χ ′ (ζ; δ)| ≤ 4|δ| + 16|δ| 2 + 10|δ| 3 + 3|δ| 4 4 -16|δ| -2|δ| 2 -4|δ| 3 (r 1 + r 2 ) as ζ → ∞,(58)Take for instance δ 1 = 1 5 . From (55) and (57), χ(ζ; δ) is uniformly bounded for ζ > 0 and there exists Z 1 > 0 such that |χ(ζ; δ)| < 2 4δ 1 ∀ζ > Z 1 ∀|δ| < δ 1 .

2 ( 2 L 2 1 / 2 (

 2222 Br * ) = R ∂ ∂r (rR r * φ) (r * , z) φ(z) dz = R 2πr * I 0 (2π|ξ|r * ) I 1 (2π|ξ|r * ) |ξ|| φ(ξ)| 2 dξ ≤ 2πr * I 0 (2π|ξ|r * ) I 1 (2π|ξ|r * ) R 1 + |ξ| 2 1 /2 | φ(ξ)| 2 dξ = 2πr * I 0 (2π|ξ|r * ) I 1 (2π|ξ|r * ) φ 2 H 1 /2 (R) . This implies R r * φ ∈ H 1 1 /2 (B r * ),andr -1 ∇(rR r * φ) Br * ) ≤ C(r * ) φ H 1 /2 (R) with C(r * ) = 2πr * I 0 (2π|ξ|r * ) I 1 (2π|ξ|r * ) .

  . A vector field a is called axisymmetric if, in the sense of distributions, ∂ θ a vanishes. According to [3, Propostion 2.2], the Maxwell equations (1) decouple into two systems, a first one for (H θ , E m ), and a second one for (H m , E θ

  1,µ 1 /2 (B r * ,z * ) is a direct application of classical theory for trace spaces: [19, Théorème 3.2]. Similar considerations apply for H 1 1 /2 (B ± r * ,z * ).

  r * ,z * ). Then it becomes clear from the variational formulation (19) that u| Br * ,z * ∈ H 1 1 /2 (B r * ,z * ) and satisfies

	Br * ,z *	1 µ	∇(ru) •	1 r	∇(rv) dr dz -	Br * ,z *	iωσuvr dr dz

Table 2 :

 2 Errors of longitudinal domain truncation with different boundary conditions. The DtN maps are truncated at N = 20.

	P P P P P P P P norm b. c.	Dirichlet	Neumann	DtN
	1 /2 • L 2 | • | H 1 1 /2 (Br * ,z * ) (Br * ,z * )	55.73% 28.79%	181.57% 49.26%	0.15% 1.86%

  ζr 1 )Y 1 (ζr 1 ) + 2r 1 πζr 2 Y 1 (ζr 1 )J 1 (ζr 1 ) -Y 0 (ζr 1 )J 0 (ζr 1 ) ,B ′ 2 (ζ) =r 2 Y 2 1 (ζr 2 ) -Y 2 0 (ζr 2 ) J 0 (ζr 1 )J 1 (ζr 1 ) + r 1 J 2 1 (ζr 1 ) -J 2 0 (ζr 1 ) Y 0 (ζr 2 )Y 1 (ζr 2 ) + 2 ζ Y 0 (ζr 2 )Y 1 (ζr 2 )J 0 (ζr 1 )J 1 (ζr 1 ),With the limiting forms of Bessel functions for small arguments [1, (9.1.7)-(9.1.9)], we gets as ζ → 0,

	A 1 (ζ) →	1 π	r 2 r 1	+	r 1 r 2	,	A ′ 1 (ζ) ∼ -	4 πζ	r 2 r 1	+	r 1 r 2	,
	A 2 (ζ) →	1 π	r 2 r 1	+	r 1 r 2	,	A ′ 2 (ζ) ∼ -	2r 1 r 2 π	ζ ln(ζ),

  1 4 , k + 1 4 ). Let δ 2 = 1 8 . From (56) and (58), χ ′ (ζ; δ) is uniformly bounded for ζ > 0 and there exists Z 2 > 0 such that

	|χ ′ (ζ; δ)| < 2	4|δ 2 + 16δ 2 2 + 10δ 3 2 + 3δ 4 2 4 -16δ 2 -2δ 2 2 -4δ 3

2 (r 1 + r 2 ) < 1 2 (r 1 + r 2 ) ∀ζ > Z 2 ∀|δ| < δ 2 .

/2 (I) * as the dual space of H θ,µ 1 /2 (I) with pivot space L

2,µ 1 /2 (I) under (•, •) µ scalar product. The definition of the spaces H θ,µ 1 /2 (Γ ± ) and H θ,µ 1 /2 (Γ ± ) * are obtained from H θ,µ 1 /2 (I) and H θ,µ 1 /2 (I) * by identifying

Acknowledgement

The authors would like to thank the two referees for theirs careful reading of the paper and their remarks that helped correcting some mistakes in the original version and greatly improved the quality of paper.

Now we shall show that for all k ∈ N * , there is one and only one zero of ϑ in the interval (k -1 2 , k + 1 2 ). Firstly, we note that t → tan(πt) is bijective from (k -1 2 , k + 1 2 ) to R. Therefore the existence can be obtained by the uniform boundedness of the function t → χ( ρ0(t) r * ) for all t > 0. Considering the fact that ρ 0 (•) is continuously increasing and (53), we only need to prove that the function χ is uniformly bounded for all ζ > 0. Secondly, we note that

If we can prove

then the uniqueness can be concluded from the monotonicity of ϑ(t) on these intervals.

Using the McMahon's expansion (53), we get

We rewrite the function χ given by (51) in the following form

We also have

We calculate the derivatives of these functions using the properties of Bessel functions and get

The problem with Dirichlet boundary condition on

Proposition C.3. Let r * > 0 be sufficiently large so that the support of the source term J ∈ L 2 (R 2 + ) is included in B r * . Then problem (61) has a unique solution u d ∈ H 1 1 /2,0 (B r * ). Assume in addition that there exists positive r 0 < r * such that the source J and the conductivity σ vanish and the permeability µ is constant for r > r 0 . Then there exists a constant C that depends only on J, r 0 , µ and σ such that

where u is the solution of (6) (in Proposition 2.2).

Proof. The proof of the first part is similar to the proof of Proposition 2.2 thanks to Lemma 2.6. Let us set

Using the lifting operator R r * , then wd := w d -R r * (u| Γr * ) satisfies

A similar argument to the proof of Proposition 2.2 thanks to Lemma 2.6 yields the the existence and uniqueness of wd ∈ H 1 1 /2,0 (B r * ), thus the existence and uniqueness of the solution w d ∈ H 1 1 /2 (B r * ). By taking v = wd in the variational formulation (63), we get the following estimates

The last inequality is due to the fact that σ vanishes for r > r 0 . Thanks to the second Poincaré-type inequality ( 16), we have

Considering (59), (60) and Proposition 2.4, we have


A Sarafianos

R Llido

O Gagliano

V Serradeil

M Lisart

V Goubier

J M Dutertre

A Tria

V Pouget

D Lewis

preference for oral

This paper presents the electrical model of a PMOS transistor in 90nm technology under 1064nm Photoelectric Laser Stimulation. The model was built and tuned from measurements made on test structures. It permits to simulate the effect of a continuous wave laser on a PMOS transistor by taking into account the laser's parameters (i.e. spot size and location, or power) and the PMOS' geometry and bias. It offers a significant gain of time by comparison with experiments and makes possible to build 3D photocurrent cartographies generated by the laser on the PMOS.

Introduction

Failure analysis (FA) methodologies make an extensive use of laser stimulation techniques. These techniques are very expensive and time consuming. In this paper, we present an electrical model of Photoelectric Laser Stimulation (PLS) of a PMOS transistor in 90nm technology. This electrical model makes it possible to simulate the response of a PMOS to PLS in a very small amount of time by comparison with real experiments on a laser equipment. Thus, the correct behaviour of a device under test may be predicted by simulation. It may serve as a "golden" reference to detect any defect by comparison and measurements. This approach permits a significant gain of time by reducing the need to use a laser source for characterization purposes. Obviously the obtained models are technology and layout dependent. Then, a new model has to be built for any technology improvement. Such models are intended to reduce the time passed on laser equipments.

This work presents the electrical model, (thanks to ELDO which is a SPICE simulator from Mentor Graphics [START_REF]ELDO user's Manual -Software Version 6.6_1 Release[END_REF]), of the photoelectric effect applied to a PMOS transistor in 90nm CMOS technology under backside continuous PLS. This model is valid for small laser power (~mW). In order to validate and correlate the model, measurements on PMOS transistors under laser illumination at a wavelength of 1064nm were performed with an I-phemos Hamamatsu equipment [START_REF]Inverted Emission Microscope iPHEMOS Series Hamamatsu[END_REF].

The novelty of the model presented in this paper is that it takes into accounts the laser's parameters (i.e. spot size and power), spatial parameters (i.e. the spot location and the PMOS' geometry), and the PMOS' bias. Therefore, the model is calibrated for each laser conditions (power, spot size, location, etc.) [START_REF] Sanchez | Delay variation mapping induced by dynamic laser stimulation[END_REF].

This article is organized as follows. Section 2 reports the behaviour under PLS of the two kinds of PN junctions encountered in PMOS transistors: P+/Nwell and Nwell/Psubstrate junctions. Measurements made on silicon permit to build an accurate electrical model. These models are used in turn to build a complete electrical model of a PMOS transistor under PLS as described in section 3. It offers the ability to draw 3D cartographies of PLS of PMOS for FA purposes as reported in section 4. Finally, our findings are summarized in the concluding section 5 with some perspectives.

Electrical modelling of the PN junctions found in PMOS transistors

In a PMOS transistor there are three PN junctions which are susceptible to give rise to a photoelectric effect if exposed to a laser beam (Drain/Nwell, source/Nwell, and Nwell/Psubstrate). Therefore to well model a PMOS transistor under laser illumination, it is first of all necessary to understand the phenomenon involved when a PN junction is stimulated by a photoelectric laser. corresponding author: alexandre.sarafianos@st.com, fax: +(33) 442.688.729, phone: +(33) 442.688.536 http://dx.doi.org/10.1016/j.microrel.2012.06.047

PN junctions measurements under PLS

The Device Under Test (DUT) is a Pwell on Nwell junction embedded in a test structure of a STMicroelectronics 90nm technology. The test series were performed with a 20X objectives and a laser spot located in the middle of the diode and focused on the active zone. The laser power at the output of the chosen objective is adjustable between 0 and 41mW.

The I(V) characteristics depicted on figure 1 were obtained for various laser power according the aforementioned settings. A strong induced photocurrent was measured for the PN junction in reverse-biasing. The wide Space Charge Region of the junction in reverse bias explains that result. At the given power settings, any variation in the bias has no effect on the photocurrent provided it is maintained in reverse-biasing. However, when the diode is forward-biased (over its threshold voltage ~0.7V) the SCR is tightened and the induced photocurrent is almost negligible by comparison with the current flow in the diode [START_REF] Glowacki | Expanding the scope of laser stimulation techniques for functional analysis and reliability of semiconductor devices by in-depth investigation of the optical interaction with the devices[END_REF].

In first approximation, and only with small laser power on this technology, it is possible to highlight the fact that in reverse mode the diode generates a constant photoelectric current independently of the voltage applied to its two electrodes. Where I laser is the current generated in the PN junction by the laser in Ampere and P laser is the laser power in W/cm².

Study of the spatial dependence of the induced photocurrent

The distance between the laser spot and the PN junction has a strong impact on the value of the generated photocurrent. In the following, this effect is studied for both P+/Nwell and Nwell/Psubstrate diodes. That is why it is important to study this effect on P+/Nwell and Nwell/Psubstrate diodes.

When the laser beam is centered on the PN junction, the photocurrent is maximum, and when the laser moves outside the junction the value of the photocurrent deacreases.

Study of the P+/Nwell junction

In order to model the spatial dependence effect on a drain/bulk diode (P+/Nwell junction), a PMOS transistor with long channel (W=L=10µm) was used. The measurements were made with the drain at 1.2V, the bulk grounded, and the two other electrodes left floating. The laser spot was moved on a line from the center of the junction to a distance of 300µm. For each step, the photocurent was measured. This experiment was conducted for the three objectives of our laser I-phemos equipement (2.5X, 20X, and 50X) [START_REF]Inverted Emission Microscope iPHEMOS Series Hamamatsu[END_REF].The obtained photocurent versus distance curves are depicted on figure 2 after normalization according the maximum photocurrent (Iph2.5X, Iph20X, and Iph50X). Their shapes exhibit a Gaussian-like behaviour. For each I-phemos objective, the coefficients a, b, c 1 and c 2 were found different (see Table I). The three functions approximated mathematically (Igauss2.5X, 20X, 50X) are also plotted (See Fig. 2).

Study of the Nwell/Psubstrate junction

For the Nwell/Psubstrate junction, the same methodology was used than for the Drain/Bulk junction. The Gaussian profiles of the laser effect on silicon is not the same than for the drain/bulk The main difference is due to the important area of the diode Nwell/Psubstrate. (20µm x 20µm). Table 2 presents the coefficient of equation ( 2) in order to model Gaussian effect of laser on Nwell/Psubstrate junction. Therefore wherever the laser beam is located on the studied device, the model is able to give the value of the photocurrent generated on the electrodes of the PMOS transistor.

Electrical model of a PN junction

In order to simulate a PN junction under PLS, the solution proposed in this paper is to define two sub circuits: Subckt Iphp to emulate the source and drain photocurrent generation, and Subckt Iph_nw_psub to simulate the laser effect on the diode Nwell/Psubstrate. These sub circuits contain a current source whose intensity value is defined for both by the parameter Iph_val (3):

gauss Ilaser S val Iph * * _ ( 3
)
Where S is the junction area in µm 2 , I laser is the value of the photocurrent generated on the PN junction (as defined in (1)) in Amper, and α gauss the spatial dependency coefficient (cf. ( 2)). Note that α gauss has different values for the two kind of diode (P+/Nwell and Nwell/Psubstrate).

The photocurrent generated on source and drain electrodes takes into account the value of the coefficient c which depends on the distance parameter d decided by the user for each junction of the PMOS, as noticed in equation [START_REF]Inverted Emission Microscope iPHEMOS Series Hamamatsu[END_REF]. The coefficient a, b, c 1 and c 2 are those define in Table 1.

For the photocurrent generated on the diode Nwell/Psubstrate the same process is used. But, this time the simulator uses the parameter from Table 2.

Electrical model of a PMOS transistor under PLS

Building the electrical model of the PMOS under PLS

After the calibration step of the two different subcircuit (Iphp and Iph_nw_psub) it is possible to model a 90nm technology PMOS transistor under PLS on the basis of manufacturer technological MOS models. On the four MOS electrodes, two sub circuits which model the P+/Nwell junctions under PLS (Sub circuits Iphp) are called by the main program. In the same way the sub circuit Iph_nw-Psub is called. (see Fig. 3).

When the user calls these sub circuits in the main netlist, for each junction he has to inform the value S (area of the junction), and d (shortest distance between the centre of the laser spot to the edge of the junction). Values of the 3 resistances of the model are defined as following (with geometrical parameters in µm defined on figure 3):

W L R R R nw nwc nwc 2 * 2 1 (5) W L R R nw nw nw * (6) W L L L R R Psub nw nw Psub psub 2 * (7)
Where R nw and R psub are the value of the resistance of each different type of silicon for one square. Therefore the model takes into account the diodes area, the different possible places of the laser beam on the device, and different PMOS' bias chosen by the user.

Static laser positioning

This electrical model is able to simulate the effect of PLS for the different states of the PMOS transistor (On/Off). Graphs Fig. 5a and 5b highlight the good correlation between ELDO simulation and measurements on a turned OFF PMOS transistor W=10µm/L=0.09µm (See Fig. 5a) and a turned ON PMOS W=10µm/L=1.2µm (See fig. 5b), when the laser beam is parked in the middle of the PMOS transistor.

Current cartographies

In this section, it is proposed to draw 3D current cartographies (extracted from the ELDO simulator) of photocurrents generated on a PMOS transistor W=L=10µm [START_REF] Llido | Improving defect localization techniques with laser beam with specific analysis and set-up modules[END_REF]. The user has the possibility to define the meshing density of cartographies thanks to a parameter called stepXY. This parameter is defined in function of the parameter distance d. It is also possible to extract 3D cartographies of photocurrents generated by the three PN junctions of a PMOS. The 3D cartographies of the photocurrent generated in the two P+/Nwell junctions highlight the Gaussian shape of the laser amplitude versus the spot location (See Fig. 6).

Due to the important area of the Nwell, 3D current cartography extracted from the ELDO simulation, (and confirmed by measures) does not have the same profile. The value of the maximal photocurrent is higher (1.45mA) (See Fig. 7).

Conclusion and perspectives

The electrical ELDO simulation (based on SPICE language) of the interaction between laser and silicon on a PMOS transistor in 90nm technology seems to be an extremely reliable, fast and also economical tool. In failure analysis it could have a great utility in order to compare the simulated response of a golden circuit under PLS with actual measurements. We have studied in a first place the modelling of PMOS transistors under PLS, because it is more complex to build than for NMOS (hence, a third PN junction between the N and P wells is encountered in PMOS). The validity of our approach is assessed by the very good correlation obtained between simulations and measurements. This work will be extended to NMOS transistors. Therefore, as a perspective, it will permit us to simulate the behaviour under PLS of complex logic gates made of PMOS and NMOS transistors.

Fig. 1 .

 1 Fig. 1. I(V) characteristics of a junction Pwell on Nwell under PLS for different laser power.

Fig. 2 .

 2 Fig. 2. Spatial dependence of the induced photocurrent (measurements and models).It makes it possible to extract a mathematical model based on the sum of two Gaussian functions (2), where d is the distance between the spot and the junction expressed in µm:

Fig 3 .

 3 Fig 3. ELDO model of the PMOS under PLS.

Fig. 5 .

 5 Fig. 5. ELDO simulation for a spot parked in the center of the device, at different laser power, with a spot size of 3.25µm. (a) PMOS transistor W=10µm/L=0.09µm turned OFF (b) PMOS transistor W=10 µm/L=1.2µm turned ON

Fig. 6 .

 6 Fig. 6. 3D cartography extracted from ELDO of: (a) the source under PLS. (b) the drain under PLS.

Fig. 7 .

 7 Fig. 7. 3D cartography of the photoelectric contribution of the Psubstrate under PLS extracted from ELDO.

Table 1

 1 Coefficients of the Gaussian function (2) for the different Iphemos objectives for the drain/bulk study.

	Coefficient	2.5X	20X	50X
	a	0.4	0.6	0.7
	b	0.6	0.4	0.3
	c 1	2.5	23.8	1000
	c 2	55	654	15000

Table 2

 2 Coefficients of the Gaussian function (2) for the different Iphemos objectives for the Nwell/Psubstrate junction study.

	Coefficient	2.5X	20X	50X
	a	0.5	0.6	0.7
	b	0.5	0.4	0.3
	c 1	2	1000	1500
	c 2	48	15000	17000