J Ér

Émi Dard É

AND Niami Nasr

Lisl Weynans

IMMERSED BOUNDARY METHOD FOR THE COMPLETE ELECTRODE MODEL IN ELECTRICAL IMPEDANCE TOMOGRAPHY

Keywords: Mathematics Subject Classification. 65N21, 65N25, 65N85, 35R30. 1

We propose an immersed boundary scheme for the numerical resolution of the Complete Electrode Model in Electrical Impedance Tomography, that we use as a main ingredient in the resolution of inverse problems in medical imaging. Such method allows to use a Cartesian mesh without accurate discretization of the boundary, which is useful in situations where the boundary is complicated and/or changing. We prove the convergence of our method, and illustrate its efficiency with two dimensional direct and inverse problems.

Introduction

This work is dedicated to the numerical resolution of direct and inverse problems related to Electrical Impedance Tomography (EIT) using an immersed boundary method. The aim of EIT is to reconstruct the electrical conductivity distribution inside a domain by imposing electrical currents on the boundary of this domain, and measuring the resulting voltages on the same boundary. It has several applications in the medical field, in particular in lung monitoring and stroke detection [START_REF] Agnelli | Classification of stroke using neural networks in electrical impedance tomography[END_REF][START_REF] Candiani | Computational Framework for Applying Electrical Impedance Tomography to Head Imaging[END_REF][START_REF] Isaacson | Reconstructions of chest phantoms by the D-bar method for electrical impedance tomography[END_REF]. Mathematically, the problem, known as Calderón problem or inverse conductivity problem, is a severely ill-posed inverse problem. We refer to [START_REF] Borcea | Electrical impedance tomography[END_REF][START_REF] Cheney | Electrical impedance tomography[END_REF][START_REF] Uhlmann | Electrical impedance tomography and Calderón's problem[END_REF] and the references therein for an overview on the inverse conductivity problem.

In practical experiments, the currents are driven inside the body of interest through a collection of surface electrodes, no current being driven between the electrodes. For each current pattern, the potential differences between the electrodes are measured. This practical setting is accurately modeled by the Complete Electrode Model (CEM) [START_REF] Cheng | Electrode models for electric current computed tomography[END_REF][START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF]. It takes into account the shape of the electrodes as well as the shunting effect, that is the thin resistive layer that appears at the interface between the electrodes and the object during the measurements. The CEM is known to correctly predict experimental data, and therefore is widely used in the numerical resolution of both direct and inverse problems related to EIT.

Usually, Finite Element Methods (FEM) are used to compute solutions of the CEM. This leads to fast and efficient numerical methods in simple geometries. However, when the geometry is complicated, the FEM can become very expensive, especially in three-dimensional problems. This is due to the need of an accurate discretization of the geometry, in particular on the electrodes where the potential varies rapidly in space, imposing the use of very refined meshes.

Another difficulty could arise when the geometry is not perfectly known, which usually leads to highly incorrect reconstructions [START_REF] Gersing | Influence of changing peripheral geometry on electrical impedance tomography measurements[END_REF][START_REF] Nissinen | Compensation of Modelling Errors Due to Unknown Domain Boundary in Electrical Impedance Tomography[END_REF]. Iterative Newton-type algorithms that reconstruct both the conductivity and the geometry have been proposed to tackle this difficulty [START_REF] Dardé | Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography[END_REF][START_REF] Dardé | Simultaneous recovery of admittivity and body shape in electrical impedance tomography: an experimental evaluation[END_REF]. In such algorithms, the geometry changes from one iteration to the other. Therefore, a new mesh is created at each step, leading again to a costly numerical implementation.

In the present paper, our objective is to propose an alternative to the numerical resolution of CEM with classical mesh-adapted numerical methods, by using Immersed Boundary Methods (IBM), with the aim of reducing the cost of computation in complex and moving geometries. In such approach, the domain of interest is included in a larger domain with a simple geometry, typically a square in two-dimensional settings, or a cube in three-dimensional settings. The geometry of the domain of interest is encoded as the zero level-set of a known function, and is not discretized exactly. Instead, the larger domain is discretized using typically a Cartesian mesh. Of course, this approach requires a particular attention to include the boundary conditions of the initial problem.

The first method on Cartesian grids for elliptic problems was designed by Mayo in 1984 [START_REF] Mayo | The fast solution of Poisson's and the biharmonic equations on irregular regions[END_REF], and developed further in [START_REF] Mayo | The rapid evaluation of volume integrals of potential theory on general regions[END_REF] and [START_REF] Mayo | Fast parallel iterative solution of Poisson's and the biharmonic equations on irregular regions[END_REF]. In that work an integral equation was derived to solve elliptic interface problems with piecewise coefficients to second order accuracy in maximum norm. Then LeVeque and [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] devised the very well known Immersed Interface Method (IIM), that relies on Taylor expansions of the solution on each side of the interface, with a local coordinate transformation near the interface to express the jump conditions in an appropriate frame. The elliptic operator is discretized on each grid point near the interface with formulas accounting for the jumps across the interface. In order to find these formulas, a linear system with six unknowns needs to be solved for each of the concerned grid points. Numerous developments of the IIM have been performed thereafter, among them [START_REF] Berthelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF][START_REF] Li | Maximum principle preserving schemes for interface problems with discontinuous coefficients[END_REF][START_REF] Wiegmann | The explicit-jump immersed interface method: Finite difference methods for PDEs with piecewise smooth solutions[END_REF]. Another class of Cartesian method was introduced by Zhou et al. is the Matched Interface and Boundary (MIB) method [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF]. The solution on each side of the interface is extended on fictitious points on the other side. These fictitious values are computed by iteratively enforcing the lowest order interface jump conditions. Other classes of Cartesian methods also exist, less accurate in the case of interface problems, but probably simpler to implement: Gibou et al. [START_REF] Gibou | A second-order-accurate symmetric discretization of the Poisson equation on irregular domains[END_REF][START_REF] Gibou | A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem[END_REF] developed methods inspired by Fedkiw's Ghost-Fluid Method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] for multiphase flows. These methods are second order accurate for Dirichlet boundary conditions on arbitrary domains, but only first order accurate for interface problems. Then new methods were proposed to increase the accuracy in case of interface problems [START_REF] Bochkov | Solving elliptic interface problems with jump conditions on Cartesian grids[END_REF][START_REF] Egan | xGFM: recovering convergence of fluxes in the ghost fluid method[END_REF].

In this paper, we approximate the CEM with a finite-difference scheme on a Cartesian mesh, and we follow the methodology of [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF], where additional unknowns were defined at the intersections of the interface with the grid. These additional unknowns were used in the discretization of the elliptic operator near the interface, which avoided deriving specific finite differences formulas containing jump terms, or corrective terms. In order to solve the interface unknowns, the flux jump conditions were discretized and added to the linear system to solve. This use of additionnal unknowns at the interface between two domains, applied in the context of electrical impedance tomography, simplifies the discretization because the values at the boundary of the domain are directly involved in the equations modeling the input currents on the electrodes.

As stated above, the main difficulties for the discretization lie in the limit conditions. In particular, the limit conditions related to the input currents are integrals on the electrodes, which are non-standard from a finite differences point of view, and therefore need particular meticulousness. Another difficulty appears at the boundary of the electrodes, where the limit condition changes abruptly, leading to singularities in the solution. We tackle this difficulty by using the smoothened Complete Electrode Model [START_REF] Hyvönen | Smoothened complete electrode model[END_REF], in which the transition between the conditions on electrodes and the condition in between electrodes in smooth, leading to smooth solutions.

Finally, as the problem is typically of Neumann type in order to accurately model the physical measurements setting, one has to deal with a compatibility condition on the data which is not easy to verify at the numerical level. We circumvent this last problem with a slight modification of the model, which allows to get rid of the compatibility condition without changing the solution of interest: this might be of interest on its own, as example to simplify numerical implementation of the CEM model by finite elements.

1.1. Outline of the paper. The paper is organized as follows. In Section 2, we present the Complete Electrode Model and its variations of interest for our study. In Section 3, we present the numerical scheme, and prove a linear rate of convergence. We present numerical results in a two-dimensional setting in Section 4, for both the direct problem of Electrical Impedance Tomography, the inverse conductivity problem, and the electrodes location problem. Some needed and technical computations are presented in the Appendix A.

The Complete Electrode Model for Electrical Impedance Tomography

2.1. Geometrical and functional settings. Here and in the following, we suppose that all functions and vectors are real-valued. Everything can be straightforwardly extended to the complex case.

Let Ω be a bounded domain of R d , d = 2 or 3, with a smooth boundary ∂Ω. We denote by ν the exterior unit normal of ∂Ω. Let (E m ) m=1,...,M be M ≥ 2 mutually disjoint connected subdomains of ∂Ω with positive Lebesgue-measure. We set E = M m=1 E m , and E c = ∂Ω \ E. We denote R M the subspace of mean-free vectors of R M ,

R M = V ∈ R M , M m=1 V m = 0 ,

and R M

+ the subset of R M of vectors with positive components,

R M + = V ∈ R M , V m > 0 ∀m ∈ {1, . . . , M } .
For s ∈ R, we introduce the quotient space

H s = (H s (Ω) × R M )/R. Two elements (u, U ) and (v, V ) of H s (Ω) × R M are identified in H s if
they differ by an additive constant, in the following sense:

(2.1) (u, U ) ∼ (v, V ) ⇔ ∃c ∈ R, u = v + c, U = V + c[1, . . . , 1].
Classicaly, H s endowed with the norm

u, U H s = inf c∈R ( u -c H s (Ω) + U -c[1, . . . , 1] R M ), is a Banach space.
For all m in {1, . . . , M }, let e m ⊂ E m be of positive Lebesgue measure. We set

b : ((u, U ), (v, V )) ∈ (H 1 (Ω) × R M ) 2 → Ω ∇u • ∇v dx + M m=1 em (u -U m )(v -V m ) ds(x),
and, for all (v, Proof. This result is standard, but we recall its proof for the reader's convenience. Firstly, we note that if (u

V ) ∈ H 1 (Ω) × R M , (v, V ) 2 b = b ((v, V ), (v, V )) .
1 , U 1 ) ∼ (u 2 , U 2 ) and (v 1 , V 1 ) ∼ (v 2 , V 2 ) in the sense of (2.1), then b ((u 1 , U 1 ), (v 1 , V 1 )) = b ((u 2 , U 2 ), (v 2 , V 2
)), so the Lemma makes sense. Clearly b is bilinear, symmetric and positive in

H 1 ×H 1 . Let (v, V ) ∈ H 1 be such that v, V b = 0. Then ∇v = 0 in Ω, and there exists c ∈ R such that v = c a.e. in Ω. Furthermore em (v -V m ) 2 ds(x) = 0, ∀m ∈ {1, . . . , M } , we obtain that V = c [1, . . . , 1]
, for the same constant c. Therefore, (v, V ) ∼ (0, 0), and b is a scalar product on H 1 .

Using the continuity of the trace application from

H 1 (Ω) to L 2 (∂Ω), it is readily seen that . b . H 1 .
Conversely, let (v n , V n ) be a sequence in H 1 such that v n , V n H 1 = 1 and v n , V n b goes to zero as n goes to infinity. Necessarily, there exists a sequence (ṽ n , Ṽn ) in

H 1 (Ω) × R M such that • (v n , V n ) ∼ (ṽ n ,
Ṽn ) in the sense of (2.1), and therefore ṽn , Ṽn b goes to zero,

• 1 ≤ ṽn 2 H 1 (Ω) + Ṽn 2 R M ≤ 2, for all n ∈ N.
As (ṽ n ) is a bounded sequence in H 1 (Ω) and ( Ṽn ) is a bounded sequence in R M , there exists a subsequence (still denoted (ṽ n , Ṽn )) and an element (ṽ ∞ , Ṽ∞ ) ∈ H 1 (Ω) × R M such that ṽn weakly converges to ṽ∞ in H 1 (Ω) and Ṽn strongly converges to Ṽ∞ in R M . As ∇ṽ n strongly converges to 0, ṽn actually strongly converges to ṽ∞ in H 1 (Ω), ∇ṽ ∞ = 0 and there exists a constant c ∞ such that ṽ∞ = c ∞ a.e. in Ω. Now, as for all m ∈ {1, . . . , M },

0 ≤ em (ṽ n -Ṽn,m ) 2 ds(x) ≤ ṽn , Ṽn 2 b ,
taking the limit as n goes to infinity leads to Ṽ∞,m

= ṽ∞|em = c ∞ , or equivalently Ṽ∞ = c ∞ [1, . . . , 1]. But as (v n , V n ) ∼ (ṽ n , Ṽn ), we have 1 = v n , V n H 1 = ṽn , Ṽn H 1 ≤ ṽn -c ∞ H 1 (Ω) + Ṽn -c ∞ [1, . . . , 1] R M ,
leading to a contradiction. Therefore . H 1 . b , and the result follows.

2.2. Generalized Complete Electrode Model. We are now in position to introduce the generalized Complete Electrode Model we will use in the following. This model differs slightly from the standard Complete Electrode Model appearing in the context of Electrical Impedance Tomography, as we allow interior and boundary source terms, and we also allow spatially varying surface admittivities on the electrodes. Both features will be useful in the following study.

In our geometrical setting, Ω is the domain of interest, and each open set E m represents one of the electrodes. We denote σ ∈ L ∞ (Ω) the conductivity of the domain, which is assumed to verify the standard ellipticity condition σ ≥ c > 0 a.e. in Ω.

In the standard Complete Electrode Model, the thin resistive layer appearing at each electrode-body interface is represented by the vector of contact impedances z ∈ R M + : one contact impedance z m per electrode E m . The corresponding contact admittivity on electrode E m is then constant, equal to 1 zm . In our generalized model, we allow spatially varying admittivities. To do so, we consider ξ m ∈ Z m , with

Z m = {ξ ∈ L ∞ (E m ), ξ ≥ 0, ξ ≡ 0} .
These varying admittivities have been introduced in [START_REF] Hyvönen | Smoothened complete electrode model[END_REF] as a way to obtain regularity on the solutions of the Complete Electrode Model. The standard Complete Electrode Model corresponds to the choice ξ m = 1 zm , which is indeed in Z m . Let I ∈ R M represent the input currents imposed on the electrodes, and (f, g) ∈ L 2 (Ω) × L 2 (∂Ω). The boundary value problem corresponding to the Complete Electrode Model reads: find (u,

U ) ∈ H 1 such that (2.2)            -∇ • (σ∇u) = f in Ω, σ∇u • ν = g on E c , σ∇u • ν + ξ m (u -U m ) = g on E m , m = 1, . . . , M, Em σ∇u • ν ds(x) = I m , m = 1, . . . , M.
Before proving that problem (2.2) is well-posed, we make several comments:

• A simple application of the divergence theorem shows that for problem (2.2) to admit a solution, it is necessary that I, f and g satisfy the compatibility condition

(2.3) M m=1 I m + Ω f dx + Ec g ds(x) = 0.
Of course, this is simply the well-known current conservation law applied to our system of equations. In the case of practical Electrical Impedance Tomography, where the source terms f and g are null, we retrieve the mean-free condition on the input currents

M m=1 I m = 0 ⇔ I ∈ R M .
• Suppose that (u, U ) ∈ H 1 (Ω) × R M satisfies the system of equations (2.2). Then for any constant c ∈ R, (u + c, U + c[1, . . . , 1]) also satisfies (2.2). This is natural as from a physics point of view, only potential differences can be measured. Therefore we work in the quotient space H 1 to identify all elements of H 1 (Ω) × R M that differ up to an additive constant. As we shall see, this restores uniqueness. • In the spirit of [START_REF] Hyvönen | Smoothened complete electrode model[END_REF], the limit conditions (2.2 2 ) and (2.2 3 ) can be written at once: indeed, denoting by 1 Em the operator that maps a function ξ ∈ L ∞ (E m ) to its extension by zero on ∂Ω, and defining the operator

Ξ : (ξ 1 , . . . , ξ m ) ∈ M m=1 L ∞ (E m ) → M m=1 1 Em ξ m ∈ L ∞ (∂Ω),
then (2.2 2 ) and (2.2 3 ) are equivalent to

σ∇u • ν + Ξ(ξ 1 , . . . , ξ M )(u -Ξ(U 1 , . . . , U M )) = g on ∂Ω.
• In practical applications in Electrical Impedance Tomography, the source terms f and g are both null, in which case the compatibility condition (2.3) becomes the usual condition I ∈ R M . Nevertheless, non-null source terms will be used in our study in the numerical tests, to ensure smoothness on u and catch the convergence rates of our numerical scheme (see subsection "Regularity results" below). • In problem (2.2), the current configuration I is imposed, and the potential configuration U is measured. From a mathematical point of view, it would be equivalent to impose U , and measure I. To apply Immersed Boundary Method to our problem, it would be simpler to impose U , as the limit conditions are then standard. Nevertheless, in practice, it is always the current that is imposed and the potential that is measured, which breaks the symmetry as the noise is greater on the measured data then on the imposed one. Therefore we focus our study on problem (2.2).

Well-posedness. We now prove that problem (2.2) is well-posed. We follow the standard steps used in the case f = 0 and g = 0, based on the variational formulation of problem (2.2), that we adapt to our context.

Proposition 2.2. Suppose f , I and g satisfy (2.3). Then problem (2.2) admits an unique solution (u, U ) ∈ H 1 .

Proof. Consider the variational problem: find (u,

U ) ∈ H 1 such that, for all (v, V ) in H 1 , (2.4) B (u, U ), (v, V ) = L (v, V ) , with B (u, U ), (v, V ) = Ω σ∇u • ∇v dx + M m=1 Em ξ m (u -U m )(v -V m )ds(x), and
L (v, V ) = Ω f v dx + M m=1 Em g (v -V m ) ds(x) + Ec g v ds(x) + I • V.
It is clear that B is bilinear and continuous on H 1 × H 1 . By definition of Z m , there exists c > 0 such that for all m in {1, . . . , M }, there exists e m ⊂ E m of positive Lebesgue measure such that ξ m ≥ c on e m . This property combined with the assumption on σ and Lemma 2.1 easily implies the coercivity of B. On the other hand, as for all (v, V ) ∈ H 1 (Ω) × R M and all constant c, we have due to the compatibility condition (2.3)

L ((v, V )) = L ((v -c, V -c[1, . . . , 1])) ,
L is a continuous linear form on H 1 . Therefore, Lax-Migram theorem [6, Corollary 5.8] implies the existence and uniqueness of (u, U ) ∈ H 1 satisfying (2.4). Remains to prove that (2.4) is equivalent to (2.2). The fact that (2.2) implies (2.4) is standard, so we focus on the reverse implication, and prove that (u, U ) ∈ H 1 solution of (2.4)

satisfies (2.2). Choosing (v, V ) ∼ (ϕ, 0 R M ) in (2.4), with ϕ ∈ C ∞ c (Ω), immediately implies that -∇ • (σ∇u) = f in Ω. Choosing then (v, V ) ∼ (v 0 , 0 R M ) with v 0 ∈ H 1 (Ω), we obtain Ω (∇ • σ∇u v 0 + σ∇u • ∇v 0 ) dx = M m=1 Em (g -ξ m (u -U m )) v 0 ds(x) + Ec g v 0 ds(x),
which implies by surjectivity of the trace operator v 0 ∈ H 1 (Ω) → v 0|∂Ω ∈ H 1/2 (∂Ω) and the definition of the conormal derivative σ∇u • ν that u satisfies (2.2) 2 and (2.2) 3 . In particular σ∇u • ν belongs to L 2 (∂Ω).

Finally, choosing (v, V ) ∼ (0, V 0 ) with V 0 ∈ R M yields

I • V 0 = M m=1 Em (g -ξ m (u -U m )) V 0,m ds(x) = M m=1 Em σ∇u • ν V 0,m ds(x) = W • V 0 , with W ∈ R M defined by W m = Em σ∇u • ν ds(x)
, which ends the proof.

Regularity results. The convergence results we obtain for the numerical scheme we develop in the next section require smooth solutions to our problem of interest. Unfortunately, in lots of situations, and in particular in the relevant ones for application, the solutions of (2.2) fail to be smooth. Indeed, suppose that σ is smooth. Due to the mixed boundary conditions (2.2) 2 , (2.2) 3 and (2.2) 4 , the potential u satisfying system (2.2) might not be as smooth as in standard elliptic problem.

More precisely, for the standard Complete Electrode Model in which ξ m is a positive constant, and in practical settings where the source terms f and g are null, it is known that (u, U ) solution of (2.2) belongs to H 2-ε for all ε > 0, but fails to be in H 2 for all input currents I except for the null one [START_REF] Dardé | Electrode modelling: the effect of contact impedance[END_REF][START_REF] Hyvönen | Smoothened complete electrode model[END_REF].

The smoothened Complete Electrode Model has been introduced in [START_REF] Hyvönen | Smoothened complete electrode model[END_REF] precisely to overcome this difficulty. It consists in replacing ξ m constant by ξ m ∈ Z m , compactly supported in E m and smooth. Then standard regularity results apply, leading to a smooth potential u. As in practice the surface admittivity is unknown, and there is no clear evidence supporting the idea that it is constant, such change of model is reasonable. Obviously, its main drawback is more complex parametrization of the contact admittivities. Our convergence results always apply to the smoothened Complete Electrode Model.

Nevertheless, regardless of the choice of ξ m , there is no explicit solution for problem (2.2) when f = g = 0, except in very particular geometric configurations. This is why we introduce the source terms f and g: they allow us to construct explicit solutions to problem (2.2) to numerically test our convergence results.

An alternative formulation well-posed in H

1 (Ω) × R d . As seen previously, problem (2.2) is well posed in H 1 .
In other words, the solution of (2.2) is defined up to an additive constant. From a numerical point of view, one has to fix that constant (or in other words, choose the ground level), which is usually done by adding a constraint on U , such that

U 1 = 0 or U ∈ R M : the solution is then searched in H 1 (Ω) × R, with R = {U ∈ R m , U 1 = 0} or R = R M ,
restoring uniqueness. Theoretically, the constraint can be imposed on the potential u (as an example, one could impose u to be mean free over Ω), but this adds computational issues when solving numerically the problem.

In this study, we propose another approach, consisting in modifying problem (2.2) in a new problem, which will be well-posed in H 1 (Ω)×R d , and which solution is one of the solutions to problem (2.2). Note that by construction no additional constraint is imposed on the new problem to obtain well-posedness.

Setting U 1 = 0. Let ε > 0 be fixed. We consider the modified problem (2.5)

           -∇ • (σ∇u) = f in Ω σ∇u • ν = g on E c σ∇u • ν + ξ m (u -U m ) = g on E m , m = 1, . . . , M, Em σ∇u • ν ds(x) + εδ m1 U m = I m , m = 1, . . . , M.
The boundary equation (2.5) 4 is equivalent to

E1 σ∇u • ν ds(x) + εU 1 = I 1 , Em σ∇u • ν ds(x) = I m , ∀m ∈ {2, . . . , M } .
In other words, only one equation of the initial problem (2.2) is modified. We claim the following Proposition 2.3. Problem (2.5) admits an unique solution (u, U ) ∈ H 1 (Ω) × R d . Furthermore, if I, f and g satisfy the compatibility condition (2.3), U 1 = 0 and (u, U ), seen as an element of H 1 , is the solution of (2.2).

Proof. We endow H 1 (Ω) × R M with its standard scalar product, with corresponding norm

(u, U ) 2 H 1 (Ω)×R M = u 2 H 1 (Ω) + U 2 R M . First, adapting the proof of Proposition 2.2, it is readily seen that problem (2.5) is equivalent to the variationnal formulation (2.6) B 1 (u, U ), (v, V ) = L (v, V ) , ∀(v, V ) ∈ H 1 (Ω) × R M .
with

B 1 (u, U ), (v, V ) = Ω σ∇u • ∇v dx + M m=1 Em ξ m (u -U m )(v -V m )ds(x) + ε U 1 V 1 and L (v, V ) = Ω f v dx + M m=1 Em g (v -V m ) ds(x) + Ec g v ds(x) + I • V.
Similarly, a minor adaptation of the proof of Lemma 2.1 shows that

(u, U ) 2 b1 = Ω |∇u| 2 dx + M m=1 |u -U m | 2 ds(x) + U 2 1 ,
is a norm on H 1 (Ω) × R M equivalent to the standard one. The assumptions on σ and ξ m immediately imply the coerciveness of the continuous bilinear form B 1 , while the linear form L is clearly continuous on H 1 (Ω) × R M . Therefore, Lax-Milgram theorem implies the existence and uniqueness of the solution of (2.6). We now suppose that I, f and g satisfy (2.3). Choosing v = 1 and V = [1, . . . , 1] leads to

ε U 1 = Ω f dx + Ec g ds(x) + M m=1 I m = 0.
As a consequence, U 1 = 0, which immediately implies that (u, U ) satisfies the system of equation (2.2), which ends the proof.

We see that solving Problem (2.5), which is a very slight variation of our initial problem (2.2), allows to retrieve the unique solution of (2.2) satisfying the additional condition U 1 = 0, without strongly enforcing this condition. Note that this holds for any value of the positive parameter ε, as long as I, f and g are compatible in the sense of (2.3).

As it allows to solve (2.2) without having to fix the ground level, from this point we will always focus on problem (2.5).

Remark 2.4. The focus on the first electrode is completely arbitrary here. A simple renumbering of the electrodes allows to use the same formulation to enforce the zero condition on any electrode.

Remark 2.5. Problem (2.5) is well-posed regardless of the compatibility of the data I, f and g. If the data are not compatible, then (u, U ) solution of (2.5) cannot solve (2.2) as the latter has no solution.

As seen in the proof of Proposition 2.3, if (u, U ) is the solution of (2.2), we always have

εU 1 = Ω f dx + Ec g ds(x) + M m=1 I m . (2.7) As a consequence, if we set Ĩ ∈ R M defined by Ĩ1 = I 1 - Ω f dx + Ec g ds(x) + M m=1 I m , Ĩm = I m , m = 1,
we see that (u, U ) satisfies (2.2) with the original source terms (I, f, g) replaced by ( Ĩ, f, g). This second set of data obviously satisfies the compatibility condition (2.3). In a sense, problem (2.5) correct the data so that they always satisfy the compatibility condition.

In the case of compatible data, the value of parameter ε has no influence on the solution of problem (2.5), and can be freely set equal to any arbitrary value, for instance one. From a numerical point of view, this is not the case. Indeed the numerical compatibility condition has no reason to be exactly satisfied by the finite-difference formulation, but only up to the order of the approximation, even if the continuous data satisfies (2.3). Therefore, due to equation (2.7), the value of ε influences the value of U 1 . This causes an interplay between the two terms in the left-hand side of (2.5) 4 that will have an influence on the amplitude of the numerical error, but not on the convergence order itself.

Setting U ∈ R M . As a final remark, we note that it is equally easy to select (u, U ) solution of (2.2) satisfying the additional mean free value condition on U , that is U ∈ R M , without strongly enforcing this condition (that is, with a problem posed in

H 1 (Ω) × R M and not in H 1 (Ω) × R M .
Indeed, it suffices to solve the following problem:

(2.8)

           -∇ • (σ∇u) = f in Ω σ∇u • ν = g on E c σ∇u • ν + ξ m (u -U m ) = g on E m , m = 1, . . . , M, Em σ∇u • ν ds(x) + ε U = I m , m = 1, . . . , M.
Here we have set

U = M k=1 U k . The corresponding variational form is B (u, U ), (v, V ) = L (v, V ) , ∀(v, V ) ∈ H 1 (Ω) × R M . with B (u, U ), (v, V ) = Ω σ∇u • ∇v dx + M m=1 Em ξ m (u -U m )(v -V m )ds(x) + ε U V ,
the continuous linear form L being unchanged. Following the same line of reasoning used in the previous section, it is not difficult to see that (2.8) is well-posed in H 1 (Ω) × R M , and that if (I, f, g) satisfies (2.3), (u, U ) solution of (2.8) satisfies the original Complete Electrode Model equations (2.2), with furthermore U = 0. Remark 2.6. From a numerical point of view, the additional term U seems not practical, as it is non-local and therefore tends to fill the matrix corresponding to the discretization of (2.8). However, it appears only on equation (2.8) 4 , and affects only the term U ∈ R M , M being the number of electrodes on the boundary of Ω, and therefore typically less than one hundred. As the potential u is usually discretized in spaces containing thousands of degrees of freedom, in practice it does not affect significantly the matrix filling.

2.4. Towards the immersed boundary method. In our numerical scheme, the domain of interest Ω is immersed into a larger domain Ω e , a square in dimension 2 and a cube in dimension 3, such that Ω ⊂ Ω e . The domain Ω e is fully discretized with a Cartesian grid, implying that some discretization points lie in Ω e \Ω, a part of the extended domain where there is no physical problem in our context. In that situation, a first idea is to simply set the solution to zero on the degrees of freedom lying outside of Ω. But this operation is cumbersome, and not appropriate in our approach, as we need smoothness in the solution on either side of the interface ∂Ω to compute quadrature formulas on the electrodes.

We therefore choose a different approach, which consists in solving an additional partial differential equation in Ω e \ Ω, similar to the problem posed in Ω. As a consequence, we introduce the following additional problem (2.9)

   -∆u e = 0 in Ω e \ Ω u e = 0
on ∂Ω e u e = u on ∂Ω, the function u being the electric potential, solution of (2.5). Clearly, such problem is well-posed in H 1 (Ω), and its solution depends continuously on u, and therefore on the EIT source terms (I, f, g). Furthermore, it will automatically inherit the smoothness of u thanks to standard elliptic regularity results.

As a conclusion, our goal is now to define an immersed boundary method to reconstruct (u, U, u e ), with (u, U ) solution of (2.5) and u e solution of (2.9).

An immersed boundary scheme for CEM

In this section we focus on the two-dimensional case, and present the discretization of problems (2.5) and (2.9) with an immersed boundary method on a Cartesian grid. The domain Ω is immersed into a larger square domain Ω e , and the boundary conditions of (2.5) are taken into account by using additional variables located on the boundary of Ω. We prove that the method converges with first order accuracy. The convergence proof is based on a discrete maximum principle, used to provide estimates of the coefficients of the inverse of the discretization matrix.

3.1. Classification of grid points. We consider a uniform Cartesian grid defined on the square Ω e . The grid spacing is denoted h. Each node of this cartesian grid is called a grid point and is denoted M ij = (x i , y j ) = (i h, j h). We denote by u h ij the approximation of the function u at the point (x i , y j ). The set of grid points located inside the domain Ω e is denoted Ω h .

To the purpose of the discretization of the CEM we need to define additional points on the boundary ∂Ω. We define the boundary point I i+1/2,j = (x i+1/2,j , y j ) as the intersection of the boundary ∂Ω and the segment [M ij M i+1j ], if it exists. Similarly, the boundary point I i,j+1/2 = (x i , y i,j+1/2 ) is defined as the intersection of the boundary and the segment [M ij M ij+1 ]. At each boundary point we create an additional unknown u h i+1/2,j or u h i,j+1/2 . The set of boundary points is denoted δΩ h . We say that a grid point is irregular if one of its direct neighbors is a boundary point, see Figure 1. On the contrary, grid points that are not irregular are called regular grid points. The set of irregular grid nodes is denoted Ω * h . The set of electrode values is denoted E h . The grid or boundary points can be denoted with letters such as P or Q rather than with indices such as M i,j if it is convenient. We also denote if more convenient x P and y P the coordinates of a point P .

3.2. Numerical scheme.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 1
. Left: regular points described by circles •, irregular points (belonging to Ω * h ) described by bullets •, right: boundary points (belonging to δΩ h ).

Elliptic operator. To discretize the elliptic operator (2.5) 1 or (2.9) 1 on each grid point, regular or not, we use a five point stencil with the grid point M i,j and its nearest neighbors, boundary or grid points, in each direction. To make this explicit, we denote u h S the value of the solution on the nearest point in the south direction, with coordinates (x S , y S ). Similarly, we define u h N , u h W and u h E and the associated coordinates (x N , y N ), (x W , y W ) and (x E , y E ). The discretization reads

-∇.(σ∇u) h i,j = -σ i+1/2,j u h E -u h ij x E -x i -σ i-1/2,j u h ij -u h W x i -x W 1 h -σ i,j+1/2 u h N -u h ij y N -y j -σ i,j-1/2 u h ij -u h S y j -y S 1 h .
with σ i+1/2,j the value of σ at point (

x i + x E 2
, y j , ). The truncation error of this discretization is second-order accurate on regular grid points, and zeroth-order on irregular grid points.

Flux boundary conditions. At each boundary point we discretize the boundary conditions (2.5) 2 or (2.5) 3 , depending if the boundary point belongs to an electrode or not.

For instance, on a boundary point I i+1/2,j , the discretization reads

σ(∇u • ν) h i+1/2,j + ξ m (I i+1/2,j ) u h i+1/2,j -U m = 0 if I i+1/2,j ∈ E m , (3.1) σ(∇u • ν) h i+1/2,j = 0 otherwise. (3.2)
with (∇u • ν) h i+1/2,j denoting the discretized normal derivative. A similar discretization is applied at each boundary point I i,j+1/2 .

The discretization of the normal derivative (∇u • ν) h i+1/2,j depends of the local geometry of the interface. As illustrated on Figure 2, the first intersection between the normal to the boundary and the grid is located on a segment: either

[M i,j , M i,j-1 ], or [M i,j-1 , M i+1,j-1 ], or [M i,j , M i,j+1 ], or [M i,j+1 , M i+1,j+1
]. The discrete normal derivative is computed as the normal derivative of the linear interpolant of the numerical solution on the triangle composed of the boundary point I i+1/2,j and the aforementioned segment.

If we denote K this triangle, (x 1 , y 1 ), (x 2 , y 2 ) and (x 3 , y 3 ) its vertices, and u 1 , u 2 and u 3 the associated values, the basis functions on the vertices for the linear interpolation write λ j (x, y) = α j x + β j y + γ j , j = 1, 2, 3, with

α j = y k -y i (x j -x k )(y j -y i ) -(x j -x i )(y i -y k ) , β j = x i -x k (x j -x k )(y j -y i ) -(x j -x i )(y i -y k ) , γ j = x k y i -x i y k (x j -x k )(y j -y i ) -(x j -x i )(y i -y k ) ,
(ν x , ν y ) being an approximation of the normal at the interface point. With these notations, the approximation of the normal derivative writes for instance for the interface point I i+1/2,j

(∇u • ν) h i+1/2,j = (u 1 α 1 + u 2 α 2 + u 3 α 3 )ν x + (u 1 β 1 + u 2 β 2 + u 3 β 3 )ν y .
This discretization is first-order accurate because it is based on a linear interpolation. Integral boundary conditions. We now discretize (2.5) 4 . More precisely, using (2.5) 3 , it is equivalent to:

Em g + ξ m (U m -u) ds(x) + ε δ m1 U 1 = I m , m = 1, . . . , M. (3.3)
From now on, we focus on the equivalent boundary condition (3.3).

We discretize this integral boundary condition (3.3) on each electrode with a first-order quadrature formula based on a first-order discrete dirac function [START_REF] Smereka | The numerical approximation of a delta function with application to level set methods[END_REF]. This discrete Dirac function is by construction positive and non-zero only on the irregular grid points, therefore it can be written as

P ∈Ω * h ω P g(P ) + ξ m (P )(U m -u P ) + ε δ m1 U 1 = I m , m = 1, . . . , M. (3.4)
with the coefficients ω P that are the weights of the first-order quadrature formula in [START_REF] Smereka | The numerical approximation of a delta function with application to level set methods[END_REF]. The truncature error of this formula is therefore first-order.

Monotonicity of the discretization matrix.

Here we aim to prove that the discretization matrix of the linear system described before, that we denote A h , is monotone. Let us first prove the following property which will be useful in the reasoning: Proposition 3.1. With the convention used for the normal to the boundary, if the minimum of v is located on a boundary point, then at this boundary point the discrete normal derivative is negative. Moreover, if the approximated normal derivative at this point is zero, then the three points values involved in the stencil are equal.

Proof. The approximation of the normal derivative is constant in space, because it is computed from a linear interpolation on a triangle. Therefore, if the minimum of v is located on a boundary point, then at this boundary point the discrete normal derivative is negative.

Proposition 3.2. The discretization matrix of the linear system described before, that we denote A h , is monotone, that is, A h is invertible and all values of A -1 h are non-negative. Proof. Let v be an array such that all coefficients of A h v are non-negative, which we denote A h v ≥ 0. We aim to prove that all coefficients of v are non-negative. To this purpose we consider the minimum of v in the whole domain: it can be either a grid point value, or a boundary point value, or an electrode potential. We detail all these cases hereafter.

• If the minimum is reached on a grid point:

In this case we denote (i 0 , j 0 ) the indices of the smallest component of v. We assume the grid point is a regular grid point. Otherwise the formula would have slightly different weights, but the reasoning would be the same. Using the elliptic operator inequality on this point we can write:

4v i0,j0 -v i0+1,j0 -v i0-1,j0 -v i0,j0+1 -v i0,j0-1 ≥ 0,
and we deduce that

v i0+1,j0 = v i0-1,j0 = v i0,j0+1 = v i0,j0-1 = v i0,j0 .
Repeating recursively this reasoning on the neighbours of (i 0 , j 0 ), then on the neighbours of the neighbours etc, we deduce that all values of v corresponding to grid or boundary points are equal to v i0,j0 . Therefore this case amounts to the case of the minimum located on a boundary point.

• If the minimum is located on a boundary point not belonging to an electrode: Without loss of generality, we assume that the minimum is located on I i+1/2,j . On this interface point we have both relationships

(σ∇v • n) h i+1/2,j
≤ 0 because the minimum is located on this point,

(σ∇v • n) h i+1/2,j
≥ 0 because all coefficients of A h v are non-negative, which leads to (σ∇v • n) h i+1/2,j = 0 and the values of the grid points involved in the discretization of the normal derivative are also equal to the minimum. Using the reasoning of the previsous subsection, it means that all grid point and boundary point values are equal to this minimum. Therefore this case finally amounts to considering the case of the minimum on a boundary point belonging to an electrode.

• If the minimum is located on a boundary point belonging to an electrode:

Without loss of generality, we assume that the minimum is located on I i+1/2,j . On this interface point we have the relationships

(σ∇v • n) h i+1/2,j ≤ 0, (σ∇v • n) h i+1/2,j + ξ m (I i+1/2,j ) v h i+1/2,j -V m ≥ 0,
the first one because the minimum is located on this point, and the second one because all coefficients of A h v are non-negative. Therefore

v h i+1/2,j -V m ≥ 0, which means that V m = v h
i+1/2,j , so the electrode potential V m is also the minimum of v. Therefore this case amounts to considering that the minimum of v corresponds to an electrode potential V m .

• If the minimum corresponds to an electrode potential V m with m = 1:

On this electrode we have the relationship:

P ∈Ω * h ω P ξ m (P ) V m -v P ≥ 0,
therefore all boundary points belonging to E m are equal to the minimum V m . Using the previous reasoning it means that all values of v corresponding to grid points, boundary points or electrode potentials are equal to the minimum.

• If the minimum corresponds to the electrode potential V 1 :

On this electrode we have the relationship:

P ∈Ω * h ω P ξ 1 (P ) V 1 -v P + ε V 1 ≥ 0, which means that ε V 1 ≥ 0.
The minimum value of v is non-negative, therefore all coefficients of v are non-negative.

Therefore we have proven that if A h v is non-negative, then v is also non-negative. This property has two implications: A h is invertible, and all values of A -1 h are non-negative.

3.4. Discrete Green functions. In the following, the letters P and Q represent indices in the linear system, representing either discretization points (on the grid or on the boundary) or electrode values.

For instance, we denote u(P ) the coefficient of the row of u with the same index than the point P .

Similarly, A h U (P ) represents the coefficient of the P -th row of the array A h U , and A h (P, Q) is the coefficient of the P -th row and Q-th column of the matrix A h . We also define by A h (:, Q) and A h (P, :) the Q-th column and the P -th row of the matrix A h . In the spirit of [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF][START_REF] Weynans | Super-convergence in maximum norm of the gradient for the Shortley-Weller method[END_REF] for each Q ∈ Ω h ∪ δΩ h ∪ E h , let us define the discrete Green's function

G h (:, Q) = G h (P, Q) P ∈Ω h ∪δΩ h ∪E h
as the solution of the discrete problem:

(3.5) A h G h (:, Q)(P ) = 0, P = Q 1, P = Q
The matrix A h being monotone, all values of G h (:, Q) are positive. In fact, the discrete Green functions are columns of the inverse of A h . We can write the solution of the discretized problem as a sum of the source terms multiplied by the values of the discrete Green functions:

u h (P ) = Q∈Ω h ∪δΩ h ∪E h G h (P, Q) (A h u h )(Q), ∀P ∈ Ω h ∪ δΩ h ∪ E h .
We can obtain bounds on the discrete Green functions with the following result:

Proposition 3.3. Let S and S be two subsets of points , W a discrete function, α > 0, β > 0 and i, j ∈ N such that:

   (A h W )(P ) ≥ 0, ∀P ∈ Ω h ∪ δΩ h ∪ E h \ S, (A h W )(P ) ≥ α -i , ∀P ∈ S, (A h W )(P ) ≥ -(β -j ), ∀P ∈ S. Then Q∈S G h (P, Q) ≤ α i W (P ) + α i β -j Q∈ S G h (P, Q), ∀P ∈ Ω h ∪ δΩ h ∪ E h .

Proof. Using the definition of the discrete Green functions, we can write

A h W (P ) ≥ A h α -i Q∈S G h (:, Q) -β -j Q∈ S G h (:, Q) (P ), ∀P ∈ Ω h ∪ δΩ h ∪ E h .

As all coefficients of A -1

h are non-negative, it leads to

W (P ) -α -i Q∈S G h (P, Q) + β -j Q∈ S G h (P, Q) ≥ 0, ∀P ∈ Ω h ∪ δΩ h ∪ E h ,
and finally we obtain the following bound:

Q∈S G h (P, Q) ≤ α i W (P ) + α i β -j Q∈ S G h (P, Q), ∀P ∈ Ω h ∪ δΩ h ∪ E h .
3.5. Estimates of discrete Green functions and convergence result. In this section, we obtain upper bounds for the discrete Green functions corresponding to the matrix arising from the discretization of (2.5)-(2.9), and deduce from them the convergence order of the numerical scheme. We make the assumption that the functions ξ m are all smooth enough so that the solution u has the required regularity for the following analysis, in particular for the truncation error estimation presented in Subsection 3.2 to be valid.

Proposition 3.4. For h small enough, the following upper bounds hold:

Q∈Ω h \Ω * h G h (:, Q) ≤ O(1), Q∈δΩ h ∪E h G h (:, Q) ≤ O(1), Q∈Ω * h G h (:, Q) ≤ O(h).
Proof. We use several discrete functions in the context of proposition 3.3 in order to obtain bounds for the different blocks of the inverse matrix A -1 h .

• We define the array W such that W ≡ 1. Therefore, all expressions of the linear system vanish excepted the quadrature formula involving the electrode E 1 :

-∇.(σ∇W

) h i,j = 0 ∀M i,j ∈ Ω h , (σ∇W • ν) h i+1/2,j -ξ m (I i+1/2,j ) W m -W h i+1/2,j = 0 ∀ I i+1/2,j ∈ δΩ h , (σ∇W • ν) h i,j+1/2 -ξ m (I i,j+1/2 ) W m -W h i,j+1/2 = 0, ∀ I i,j+1/2 ∈ δΩ h , P ∈Em ω p ξ m (P ) W m -WP = 0 ∀ m = 1, P ∈Em ω p ξ m (P ) W m -WP + εW 1 = ε.
G h (:, E 1 ) representing the discrete Green function associated with the electrode E 1 , we obtain

G h (:, E 1 ) = 1 ε . (3.6)
• We consider the exact solution ū of (2.5)-(2.9), with f = 1, g = 1, I m = 1 ∀ m = 1, and I 1 such that the compatibility condition (2.3) is satisfied. We define the array W as the discretisation of ū on the grid points, boundary points and electrode values. The discretization of the elliptic operator and the fluxes is consistent at least with first-order accuracy, excepted for the grid points in Ω * h where it is only zeroth-order accurate. Thus for h small enough, we can write that

-∇.(σ W ) h i,j ≥ 1 2 , ∀M i,j ∈ Ω h \ Ω * h , -∇.(σ W ) h i,j ≥ -C 1 , ∀M i,j ∈ Ω * h , (σ∇ W • ν) h i+1/2,j -ξ m (I i+1/2,j ) Wm -Wi+1/2,j ≥ 1 2 , ∀I i+1/2,j ∈ δΩ h , (σ∇ W • ν) h i,j+1/2 -ξ m (I i,j+1/2 ) Wm -Wi,j+1/2 ≥ 1 2 , ∀I i,j+1/2 ∈ δΩ h P ∈Em ω p ξ m (P ) Wm -WP ≥ 1 2 , ∀ m = 1 P ∈E1 ω p ξ m (P ) W1 -WP + ε W1 ≥ -2|I 1 |,
with C 1 a strictly positive constant. It can also be re-written as:

(A h W )(P ) ≥ 1 2 , ∀P ∈ (Ω h \ Ω * h ) ∪ (E h \ E 1 ), (A h W )(P ) ≥ 1 2 , ∀P ∈ δΩ h , (A h W )(P ) ≥ -C 1 , ∀P ∈ Ω * h , (A h W )(E 1 ) ≥ -2|I 1 |,
and it leads to:

Q∈(Ω h \Ω * h )∪(E h \E1) G h (:, Q) + Q∈δΩ h G h (:, Q) ≤ 2 W + 2C 1 Q∈Ω * h G h (:, Q) + 4|I 1 |G h (:, E 1 ). (3.7)
• We define the array W such that W = 1 for all points in Ω h , and W = 0 for all points in δΩ h ∪ E h . Because of the jump in the values of W between Ω h and δΩ h ∪ E h , when we apply the discretization to W we obtain

-∇.(σ∇ W ) h i,j = 0 ∀M i,j ∈ Ω h \ Ω * h , -∇.(σ∇ W ) h i,j ≥ C 2 h 2 ∀ M i,j ∈ Ω * h , (σ∇ W • ν) h i+1/2,j -ξ m (I i+1/2,j ) Wm -Wi+1/2,j ≥ - C 3 h ∀ I i+1/2,j ∈ δΩ h , (σ∇ W • ν) h i,j+1/2 -ξ m (I i,j+1/2 ) Wm -Wi,j+1/2 ≥ - C 3 h , ∀ I i,j+1/2 ∈ δΩ h P ∈Em ω p ξ m (P ) Wm -WP + εδ m1 W1 = 0 ∀ m,
with C 2 and C 3 trictly positive constants. We thus obtain

C 2 h 2 Q∈Ω * h G h (:, Q) ≤ W + C 3 h Q∈δΩ h G h (:, Q),
and then:

Q∈Ω * h G h (:, Q) ≤ h 2 C 2 W + h C 3 C 2 Q∈δΩ h G h (:, Q). (3.8)
Combining 3.8 with 3.7 we obtain

Q∈(Ω h \Ω * h )∪(E h \E1) G h (:, Q) + Q∈δΩ h G h (:, Q) ≤ 2 W + 2C 1 h 2 C 2 W + 4|I 1 |G h (:, E 1 ) (3.9) + 2C 1 C 3 C 2 h Q∈δΩ h G h (:, Q), (3.10)
and for h small enough it leads to

Q∈(Ω h \Ω * h )∪(E h \E1) G h (:, Q) + Q∈δΩ h G h (:, Q) ≤ C (3.11)
with C a positive constant. Re-injecting this inequality into 3.8 we obtain

Q∈Ω * h G h (:, Q) ≤ h 2 C 2 W + h C 3 C 2 C. (3.12)
The inequalities 3.6, 3.11 and 3.12 provide the bounds of the proposition. Proposition 3.5. If we denote ū the exact solution of (2.5)-(2.9) and u h the numerical solution, the local error |ū(P ) -u h (P )| satisfies for h small enough:

|ū(P ) -u h (P )| ≤ O(h). (3.13)
Therefore the numerical solution converges with first-order accuracy to the exact solution in L ∞ -norm.

Proof. We denote τ the truncation error array arising from the discretization of (2.5)-(2.9). We have

|ū(P ) -u h (P )| = Q∈Ω h ∪δΩ h ∪E h G h (P, Q)τ (Q) ≤ Q∈Ω * h G h (P, Q)τ (Q) + Q∈Ω h \Ω * h G h (P, Q)τ (Q) + Q∈δΩ h ∪E h G h (P, Q)τ (Q) ,
which, according to the previous results, immediately implies

|ū(P ) -u h (P )| ≤ O(1) Q∈Ω * h G h (P, Q) + O(h 2 ) Q∈Ω h \Ω * h G h (P, Q) + O(h) Q∈δΩ h ∪E h G h (P, Q) ≤ O(1)O(h) + O(h 2 )O(1) + O(h)O(1) = O(h).

Numerical experiments

We now study the efficiency of the immersed boundary method for Electrical Impedance Tomography introduced in Section 3.2. More precisely, we test the method on three different two-dimensional problems:

• the direct problem of EIT, which in our context reads: for given source terms f , g and I, compute (u, U ) solution of (2.5) (in the sense of Proposition 2.3). In that case, the geometry (i.e. Ω and the electrodes' position) and the conductivity are supposed to be known parameters. This is a well-posed problem, and our numerical method is designed to obtain a good reconstruction of its solution without accurate discretization of the geometry. • the inverse problem of EIT: for given input currents I ∈ R M , and corresponding measured voltages U , reconstruct the interior conductivity σ. In other word, knowing the couple (I, U ), we search for σ such that (u, U ) = (u(σ), U (σ)) solves (2.5) with source terms f = 0, g = 0 and I. Here again, the geometry is supposed to be known. This of course is the central problem in Electrical Impedance Tomography, and therefore is the main validation of our method. • a geometric inverse problem in EIT: for given input currents I ∈ R M , and corresponding measured voltages U , reconstruct the positions of the electrodes knowing Ω and the conductivity σ.

Obviously, this problem does not make much sense from a practical point of view. But on one hand it can be seen at a first step to tackle the real problem of imperfectly known geometrical setting in EIT, which is known to severely deteriorate the quality of most conductivity's reconstruction algorithms (at the exception of the ones that are specifically designed to cope with this problem, see [START_REF] Agnelli | Simultaneous reconstruction of conductivity, boundary shape, and contact impedances in electrical impedance tomography[END_REF][START_REF] Dardé | Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography[END_REF][START_REF] Kolehmainen | Calderón's inverse problem with an imperfectly known boundary and reconstruction up to a conformal deformation[END_REF][START_REF] Kolehmainen | The inverse conductivity problem with an imperfectly known boundary[END_REF] and the references therein). On the second hand, it is a problem of interest in our context, as the electrodes are not exactly approximated in the immersed boundary method we proposed. Therefore, knowing the high-sensibility of EIT-measurements to electrodes' position, it is a challenging problem. The shape parameters [α 0 , . . . , α 2N ] are chosen so that 0 < r(θ) < 2 for all θ ∈ [0, 2π]. Finally, we set

Ω = ρ[cos(θ), sin(θ)] ∈ R 2 , θ ∈ [0, 2π], 0 ≤ ρ < r(θ 0 ) .
Note that, with such parameterization of Ω, the outward unit vector at a point x = x [cos(θ x ), sin(θ x )] ∈ ∂Ω, denoted ν(x), is given by

ν(x) = 1 ρ(θ x ) [r(θ x )u(θ x ) -r (θ x )v(θ x )] ,
with v(θ) = [-sin(θ), cos(θ)] and ρ(θ) = r(θ) 2 + r (θ) 2 . We also set

τ (x) = 1 ρ(θ x ) [r (θ x )u(θ x ) + r(θ x )v(θ x )] ,
the vector tangential to ∂Ω at the point x, such that (ν(x), τ (x)) is a direct basis of R 2 .

In the following tests, we use three different geometries for Ω: Ω 1 is a disk of radius 1.5, corresponding to α = [1.5], Ω 2 corresponds to the choice α = [1.51, 0.01, 0.05, 0.2, 0.035, 0.01, 0.1], and Ω 3 corresponds to the choice α = [1.6, 0.002, 0.01, 0.003, 0.035, 0.2, 0.15].

The M -electrodes E 1 , . . . , E M are parameterized by two vectors (Θ

1 , Θ 2 ) ∈ R M × R M verifying Θ 1 1 < Θ 2 1 < Θ 2 2 < . . . < Θ 2 M < Θ 1 1 + 2π, such that E m = r(θ)u(θ), θ ∈ [Θ 1 m , Θ 2 m ]
. As in practical application the length of the electrodes is prescribed, we usually set Θ 1 , choose a length size L > 0, and compute Θ 2 by solving the equation In this article we use a specific polar parametrization of the interface, which allows us to compute explicitly all normals to the interface. However, such parametrization implies that the domain of interest is star-shaped, which might not be the case in certain applications.

In such situations, it is usually still possible to use a description of the interface with a level-set function. Such a description amounts to represent implicitly the interface as the zero level set of a function. The level set method was introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. We refer the interested reader to [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Sethian | Evolution, implementation, and application of level set and fast marching methods for advancing fronts[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF].Although any smooth enough function cancelling at the interface location can be used, in practice the signed distance to the interface is very often used:

(4.2) ϕ(x) =    dist ∂Ω (x)
outside of the interface ∂Ω, 0 on the interface, -dist ∂Ω (x) inside of the interface ∂Ω.

We recall here some properties occurring in this case:

• As recalled in [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF], the level-set, being a distance function, is 1-Lipshitz and almost everywhere differentiable. Moreover, if φ is differentiable at a point x, then it satisfies the so-called Eikonal equation at x: ∇ϕ(x) = 1. • The smoothness of the level-set is in fact strongly related to the smoothness of the interface:

as proved in [25, p.355], if the interface is C 2 , then there exists a real r 0 > 0 such that the level-set is C 2 for all x, y such that |ϕ(x, y)| < r 0 . • The outward normal vector of the isoline of ϕ passing on x, denoted ν(x), can be expressed where ϕ is differentiable as

(4.3) ν(x) = ∇ϕ(x) ∇ϕ(x) , 4.2. Direct problems.
First of all, we investigate the ability of the proposed immersed boundary method to efficiently compute solutions of the Electrical Impedance Tomography problem (2.5). We highlight that the convergence results announced in Section 3.5 are indeed observed in practice for solutions smooth enough, and present several examples of reconstruction of the inner potential.

4.2.1.

Convergence of the method. To verify that the convergence rates proved in Section 3.5 are indeed observed in practice, we proceed as follows:

we set Ω e = [-2, 2] × [-2, 2], u(x, y) = sin(x y), U = 0 R 16 ,
and compute the corresponding source terms f , g and I in (2.2). In particular, to obtain the vector of input currents I, we compute the corresponding integrals appearing in (2.2) 4 using a numerical integration scheme.

We then use the immersed boundary method described in Section 3.2 to obtain an approximate solution u h to problem (2.2) for different values of the mesh parameter h. To test also the convergence of the derivatives, we compute an approximation of the gradient du h by numerical differentiation.

In the following, we compute u h and du h for various values of h ranging from 1 25 to 1 225 . The parameter ε in (2.5) 4 is chosen small, of order ∼ 10 -10 . The results are displayed in Figure 4: we observe a linear convergence both for u and its gradient.

Remark 4.2. The linear convergence of u h to u was expected. The linear convergence of the gradient is not surprising, as the construction of the numerical method is based on a first-order discretization of fluxes near the boundary of the domain. Nonetheless, to prove it is not obvious and non-standard within the finite difference framework. We left the study of the convergence of the gradient for future works.

Remark 4.3. We have presented here numerical results for a given value of the parameter ε. Let us remark that changing the value of the parameter does not change the convergence behaviour. However, the value of ε affects the amplitude of the error at fixed h. This is expected since the value of U 1 depends on the numerical compatibility condition, and therefore is not exactly zero. Therefore, there is an interplay between the therm εU 1 and the integral term in the formula (2.5) 4 that influences the value of both terms.

On Figure 5, we present a comparison of numerical errors for the generalized CEM, with smooth spatially varying admittivities, and for the classical CEM with ξ m = 1 z m constant. The geometry considered is a circle of radius R = 0.5 embeded in a square [-1, 1] × [-1, 1], with 4 electrodes equally spaced.The discretization parameter h ranges from 1 50 to 1 250 . The convergence rate displayed on the Figure is computed with the classical CEM. We observe that the amplitude of the numerical errors for both models as well as their rate of convergence is very close. Finally, we present some reconstructions of the inner potential, with Ω e = [-2, 2]×[-2, 2] discretized with a 700 × 700 Cartesian grid, and in different geometries, with different background conductivities and various patterns applied to the electrodes, see Figures 6 and7 To simplify notations, we gather all input currents and electrodes' potentials corresponding to a certain conductivity σ in M × P matrices I = I 1 . . . I P , U = U 1 , . . . , U P .

The measurement matrix U depends linearly on the input currents I and non-linearly on the unknown conductivity σ. We will sometimes make this dependency explicit by writing U (σ, I ). We also note U meas ≈ U (σ * , I ), the matrix which would idealistically corresponds to our actual measurements, but in practice is only close to it, as in applications measurements are always corrupted by noise. The inverse conductivity problem can be formulated in the following concise way: knowing I , U meas reconstruct σ * .

Over time, sophisticated methods have been introduced to reconstruct either the conductivity or the support of inclusions from (continuous or electrodes' like) boundary measurements [START_REF] Brühl | Recent progress in electrical impedance tomography[END_REF][START_REF] Cheney | NOSER: An algorithm for solving the inverse conductivity problem[END_REF][START_REF] Garde | Series reversion in Calderón's problem[END_REF][START_REF] Garde | Convergence and regularization for monotonicity-based shape reconstruction in electrical impedance tomography[END_REF][START_REF] Harrach | Recent progress on the factorization method for electrical impedance tomography[END_REF][START_REF] Harrach | Monotonicity-based shape reconstruction in electrical impedance tomography[END_REF][START_REF] Ikehata | Reconstruction of the support function for inclusion from boundary measurements[END_REF][START_REF] Kaipio | Statistical and computational inverse problems[END_REF][START_REF] Lechleiter | The factorization method applied to the complete electrode model of impedance tomography[END_REF][START_REF] Mueller | The D-bar method for electrical impedance tomographydemystified[END_REF][START_REF] Vauhkonen | Tikhonov regularization and prior information in electrical impedance tomography[END_REF]. In our present study, we will use a comparatively very naive approach to reconstruct σ * , that is a basic least-square minimization approach combined with Tychonov regularization. We choose this on purpose, as our primary objective is to test the ability of our direct method to give sufficiently good approximation of the solution to problem (2.5) to be efficiently used in an standard algorithm for the inverse conductivity problem.

More precisely, we seek to minimize the discrepancy between the solution of the direct problem (2.5), corresponding to a guessed conductivity σ, and the measured data, i.e., we seek σ such that the functional (4.4)

F (σ) = U (σ, I) -U meas 2 2
is minimal. The problem being ill-posed, a regularization is needed. As a consequence, we actually minimize the functional (4.5)

F (σ) = U (σ, I) -U meas 2 R M + 2 R(σ)
where R is a properly chosen regularisation functional and > 0 a regularisation parameter. In other words, we use a deterministic Tychonov regularization of our problem.

From now on, we assume that σ * ∈ L ∞ ∩ H 1 (Ω). We also assume that there exists a known conductivity σ ∈ L ∞ ∩ H 1 (Ω) such that σ * -σ ∈ H 1 0 (Ω). This is typically the case in applications where one searches perturbations of a known background conductivity σ , assuming that the perturbations do not touch the boundary of the domain (in that situation, the pertubations are σ * -σ ). It is then natural to seek to minimize the following functional (4.6)

F (σ) = 1 2 U (σ, I) -U meas 2 R M + 2 σ -σ 2 H 1 (Ω) .
The strategy of minimization is simple: for a conductivity σ and a current pattern I ∈ R M , we denote (u(σ, I), U (σ, I)) the solution of (2.5) with f = g = 0. Then, for a given conductivity σ n , define Finally, we define δσ n as the unique v ∈ H 1 0 (Ω) such that

-∆v + v = ∆(σ n -σ ) -(σ n -σ ) + d x u T n d x w n + d y u T n d y w n .
Then it is not difficult to prove the following result. Lemma 4.4. Suppose that δσ n belongs to L ∞ (Ω). There exists T n > 0 such that for all t ∈ (0, T n ),

F (σ n + tδσ n ) ≤ F (σ n ).
The proof of this point is standard, and is given in the Appendix.

Our algorithm of minimization is then the following:

• Choose σ 0 ∈ L ∞ ∩ H 1 (Ω) (for example, σ 0 = σ ),
• Until the following stopping criteria is reached:

|F (σ n + t n δσ n ) -F (σ n )| < τ |δσ n |, with τ small enough
(1) For a given conductivity σ n , compute δσ n as described above, (2) Define σ n+1 = σ n + t n δσ n with t n = arg min t∈(0,Tn) F (σ n + tδσ n ).

• End.

Note that, at each step n of this algorithm, we have to solve P problems (2.5) to compute d x u n and d y u n , then again P problems (2.5) to compute d w w n and d y w n , and a single Laplace-type problem to compute δσ n , that is 2P + 1 elliptic problems. Luckily, the first 2P problems are exactly (2.5) with a fixed given conductivity σ n , and only the right-hand side changing, meaning that the corresponding system to inverse is the same for the 2P problems. Hence the computation of δσ n is actually very fast.

Our numerical experiments show that, for t close to 0, the function t → F (σ n + tδσ n ) behave like a quadratic function in t with a single local minimum. We therefore search for t n using a Golden-section search.

We now present some conductivity's reconstructions in the domains presented in Figure 3. In each test, the background conductivity σ is equal to one, and

• first case: we place a circular inclusion in the center.

• second case : we place a similar inclusion closer to the boundary.

• third case : we place two inclusions of different amplitudes inside the domain. The measurement matrix U meas is obtained by solving (2.5) using our numerical method, in a very refined mesh and for P = 15 linearly independent current patterns. The mesh used to compute U meas differs from the mesh used to reconstruct the conductivity in order to avoir inverse crime. To obtain the noisy data, we then perturb the measurements with an additive Gaussian noise, scaled so that the relative amplitude of the noise is set to δ ∈ [0, 1]. Finally, we set the regularization parameter to = 10 -4 .

Figures 8 and9 present reconstructions in the different geometrical domains with various levels of added Gaussian noise, discretized by a 400 × 400 Cartesian grid. As can be seen, the reconstructions are correct even with this very naive approach, meaning that our immersed boundary methods can be efficiently used in the inverse conductivity problem. Remark 4.5. In practical applications, the conductance ξ m (or equivalently the contact impedances z m ) are unknown, and should be also reconstructed. This can be done without additional difficulties, using that (u, U ) solution of (2.5) is Frechet differentiable with respect to ξ m , with an explicit expression of the Frechet derivative. 4.3.2. Reconstruction of the electrodes' positions, knowing the conductivity. We now focus on the geometric inverse problem described in the introduction, which can be summarized as follows: from the knowledge of (I , U meas ), find the locations of the electrodes, in other words the parameters Θ 1 and Θ 2 , all the other parameters (α, σ, ...) being known.

As already said, even if this problem makes little sense from a practical point of view, it is interesting and challenging for our method as (2.2) is known to be very unstable with respect to the geometrical setting, and our Boundary Immersed Method is, by design, not precise in terms of geometry of ∂Ω.

Adapting the notations of previous Section 4.3.1, we denote by U = U (Θ 1 , Θ 2 , I) the (second term of the) solution of (2.5) with f = g = 0, I ∈ R M , and electrodes' positions determined by the angles Θ 1 and Θ 2 (see Section 4.1 for the parametrization of the geometry). Now, for a set of input currents (I 1 , . . . , I P ) ∈ (R M ) P , we denote

I = I 1 , . . . , I P , U = U 1 , . . . , U P .
We also set U meas to be the following M × P possibly noisy measurements matrix, approximating the exact measurement matrix corresponding to the exact electrodes positions Θ 1 * and Θ 2 * : U meas ≈ U (Θ 1 * , Θ 2 * , I ). Our strategy is similar as before: we seek to minimize the functional

F (Θ 1 , Θ 2 ) = U (Θ 1 , Θ 2 , I) -U meas 2 2 + 2 Θ 1 -Θ 1 2 R M + 2 Θ 2 -Θ 2 2 R M .
As previously, Θ k has to be think as a good guess for the actual location of the k-th electrode, around which we search the actual position of the electrode. The partial derivatives of the functional with respect to parameters Θ 1 k or Θ 2 k are easily derived from the shape derivative obtained in [START_REF] Dardé | Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography[END_REF], and in particular the sampling formula [15, Corollary 3.4]: Lemma 4.6. Let (ũ, Ũ ) be the solution of (2.5) for some electrode pattern Ĩ ∈ R M , and f = g = 0. For k in {1, 2} and m in {1, . . . , M }, let

x k m = r(Θ k m )[cos(Θ k m ), sin(Θ k m )].
We have for m ∈ {1, . . . , M }:

(4.7) ∂U ∂Θ 1 m • Ĩ = ρ(Θ 1 m )(U m -u(x 1 m ))( Ũm -ũ(x 1 m )) and ∂U ∂Θ 2 m • Ĩ = ρ(Θ 2 m )(U m -u(x 2 m ))( Ũm -ũ(x 2 m )),
We refer to Section 4.1 for the definition of r and ρ, and postpone the proof of Lemma 4.6 to the Appendix. The main difficulty with this sampling formula are the terms u(x k m ) and ũ(x k m ), which corresponds to the value of the interior potential at the end-points of the electrodes, and are difficult to compute accurately for two reasons: first, from a theoretical point of view, this potential is continuous at the end points of the electrodes, but no more than that, and may varies rapidly which makes it difficult to compute accurately. Furthermore, as in our immersed method the electrodes are not exactly discretized, such value is not directly available. We approximate it by linear interpolations from values of the nearest interface points located inside the electrode. Nevertheless, we obtain sufficiently precise results to reconstruct the position of the electrodes from the measurements by minimizing the functional F using a gradient descent algorithm.

We first present reconstructions with Ω a disc of radius 0.5 covered with 4 electrodes.We either perturbed the first angle or the second angle of the first electrode, the other electrodes being unperturbed, but we nevertheless let all the angles be free in our algorithm. As expected, the electrodes E 2 , E 3 and E 4 , being already at their correct positions, are almost not moved by the algorithm, while the position of the first electrode is retrieved as shown in the following table: In our final experiments, we suppose that the length of the electrodes is fixed and known. In that case, Θ 1 and Θ 2 are not independent, as they are linked by relation (4.1). The method is easily adapted to this situation: we simply decide that Θ 2 is now a function of Θ 1 given by relation (4.1), and U depends only on I and Θ 1 . The derivatives are then obtain by simple combination of Lemma 4.6, the chain rule, and the fact that equation (4.1) immediately implies

Searched angles Starting angles

Reconstructions Θ 1 1 Θ 2 1 Θ 1 1 Θ 2 1 Θ 1 1 Θ 2 1 -2.35619 -
∂Θ 2 k ∂Θ 1 l = δ kl ρ(Θ 1 l ) ρ(Θ 2 k )
.

For our final tests, we choose Ω e = [- Therefore, assuming that v belongs to L ∞ (Ω), we can choose δσ = v, which ends the proof of Lemma 4.4.

Remark A.1. Clearly, v solution of (A.3) might fail to be in L ∞ (Ω). However, from a numerical point of view, we will obtain an approximation of v which will, by construction, be in L ∞ (Ω). Such approximation turns out to be a good direction descent, as shown in the presented reconstructions.

Electrode position reconstruction -Proof of Lemma 4.6. Lemma 4.6 is a consequence of a more general result, that is the Frechet derivability of the measurement map with respect to small perturbation of the boundary of Ω, obtained in [START_REF] Dardé | Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography[END_REF]. We start by very briefly recalling the results proved in this paper. For h ∈ C 1 (∂Ω, R d ), we set F [h] : x ∈ ∂Ω → x + h(x) and

∂Ω h = F [h](∂Ω).
For h small enough, ∂Ω h is the boundary of a smooth domain Ω h , which of course is a perturbation of Ω. For h still small enough, Ω h is covered by M well separated electrodes E m,h defined by

E m,h = {x + h(x), x ∈ E m } .
As a consequence, for any I ∈ R M , we can denote (u(h), U (h)) the solution of (2.5) with (Ω, E m ) replaced by (Ω h , E m,h ), and f = g = 0, and define the measurement map

R : (h, I) ∈ B d × R M → U (h), with B d = {h ∈ C 1 (∂Ω, R n ); h C 1 (∂Ω,R d ) < d},
d being a small enough fixed constant. The main result of [START_REF] Dardé | Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography[END_REF] is the Frechet derivability of R with respect to h ([15, Theorem 3.2]). A sampling formula is also provided ([15, Corollary 3.4]), allowing to easily compute this derivative. These results are gathered in the following Theorem.

Theorem A.2. The operator R is Fréchet differentiable at the origin with respect to the first variable. In other words, there exists a bounded bilinear operator

R : C 1 (∂Ω; R d ) × R M → R M such that lim h→0 1 h C 1 R(h, •) -R(0, •) -R h L = 0.
Furthermore, one can compute R h component by component using the following sampling formula: for I and Ĩ in R M , if one denotes (u, U ) and (ũ, Ũ ) the corresponding solutions to (2.5), then Here, H is the mean curvature, ν (resp. ν ∂Em ) is the outward normal vector to Ω (resp. to E m ), and

h ν = h • ν.
Now that we have all the necessary tools, we can prove Lemma 4.6. In the following, we focus on Θ 2 , the computations with Θ 1 being similar. We also do the computation for a single input current I ∈ R M , the extension to P input currents being straightforward. Finally, for any x ∈ ∂Ω, we define θ(x) as the unique element of [-π, π[ such that x = r(θ(x))u(θ(x)).

Let m ∈ {1, . . . , M } be fixed, and η ∈ R sufficiently small such that

Θ 1 m < Θ 2 m + η < Θ 1 m+1 ,
where, to simplify notations, we set Θ 1 M +1 = Θ 1 1 + 2π. We define Θ 2 η by Θ 2 η,k = Θ 2 k + δ km η, ∀k ∈ {1, . . . , M } . We define, for a fixed constant c > 0,

H η =    h ∈ C 1 (∂Ω; R

Lemma 2 . 1 .

 21 The bilinear form b defines a scalar product on H 1 , with corresponding norm . b equivalent to the norm . H 1 .As a consequence, (H 1 , b) is a Hilbert space.

Figure 2 .

 2 Figure 2. All possible stencils for the first-order flux discretization on the left side of the interface, with points involved in the discretization signaled by black circles.

4. 1 .

 1 Representation of the geometry. In the following experiments, the domain Ω is represented by a polar parametrization of its boundary ∂Ω. More precisely, for a given vector α = [α 0 , α 1 , . . . , α 2N ] ∈ R 2N +1 , we set ∂Ω = {r(θ)u(θ), θ ∈ [0, 2π]} , with u(θ) = [cos(θ), sin(θ)], and r(θ) = α 0 + N k=1 (α k cos(kθ) + α k+N sin(kθ)) .

  Unless otherwise specified, in the experiments below we use 16 electrodes, defined by Θ 1 k = -π+(k-1) π 8 and L = 0.35.

Figure 3 .

 3 Figure 3. The domains Ω 1 , Ω 2 and Ω 3 (in red) in the reference square [-2, 2]×[-2, 2] (in blue), with our usual configuration of 16 electrodes (in bold black)

  Figure is computed with the classical CEM. We observe that the amplitude of the numerical errors for both models as well as their rate of convergence is very close.Finally, we present some reconstructions of the inner potential, with Ω e = [-2, 2]×[-2, 2] discretized with a 700 × 700 Cartesian grid, and in different geometries, with different background conductivities and various patterns applied to the electrodes, see Figures6 and 7.

  u h -u ∞ as a function of h. Ω = Ω1. du h,1 -∂xu ∞ as a function of h. Ω = Ω1.

  u h -u ∞ as a function of h. Ω = Ω2. du h,1 -∂xu ∞ as a function of h. Ω = Ω2.

  u h -u ∞ as a function of h. Ω = Ω3.

  du h,1 -∂xu ∞ as a function of h. Ω = Ω3.

Figure 4 .

 4 Figure 4. Convergence study with manufactured solution u(x, y) = sin(x y) and U = 0 R 16 in each domain.

  (a) u h -u ∞ as a function of h.(b) du h,1 -∂xu ∞ as a function of h.

Figure 5 .Figure 6 .

 56 Figure 5. Convergence study with manufactured solution u(x, y) = exp(x 2 +y 2 ) and U = 0.5 exp(R 2 )1 R 4 . Blue points stand for the numerical errors of classical CEM, red crosses for the numerical errors of generalized CEM

Figure 7 .

 7 Figure 7. Resolution of Problem (2.5) in domain Ω 3 with input current I m = δ m 8δ m 14 . Left: inner potential. Right: conductivity.

  (d x u n , d y u n ) by d x u n = ∂u(σn,I 1 ) ∂x , . . . , ∂u(σn,I P ) ∂x ∈ L 2 (Ω) d , d y u n = ∂u(σn,I 1 ) ∂y , . . . , ∂u(σn,I P ) ∂y ∈ L 2 (Ω) d . Now, with U meas = U 1 meas , . . . , U P meas , we define (d x w n , d y w n ) by d x w n = ∂u(σn,U (σn,I 1 )-U 1 meas ) ∂x , . . . , ∂u(σn,U (σn,I P )-U P meas ) ∂x and d y w n = ∂u(σn,U (σn,I 1 )-U 1 meas ) ∂y , . . . , ∂u(σn,U (σn,I P )-U P meas ) ∂y .

Figure 8 .

 8 Figure 8. Reconstructions in Ω 1 . First column: searched conductivity. Second column: reconstructed conductivity with unperturbed data (δ = 0). Third column: reconstructed conductivity with perturbed data (δ = 2.0%).

Figure 9 .

 9 Figure 9. Reconstructions in Ω 3 . First column: searched conductivity. Second column: reconstructed conductivity with unperturbed data (δ = 0). Third column: reconstructed conductivity with perturbed data (δ = 2.0%).

  ν ∂Em )(U m -u)( Ũm -ũ)ds(x) h ν ((d -1)(U m -u)H -∂ ν u) ( Ũm -ũ)ds(x) -∂Ω h ν (σ∇u) τ (∇ũ) τ ds(x).

 2 .

 2 , . . . , M } , h(Θ k2 ) = δ km η ∀k ∈ {1, . . . , M } , h C 1 ≤ cη. The set H η is non empty: indeed, take χ ∈ C ∞ c (R; [0, 1]) such that supp(χ) ∈] -1, 1[ and χ(0) = 1 and defineϕ η (θ) = ηχ θ -Then set ω η : R → R 2 π-periodic, defined byω η : θ ∈ [Θ 1 1 , Θ 1 1 + 2 π[-→ θ + ϕ η (θ).The, if the constant c is chosen sufficiently large, the vector fieldh : x ∈ ∂Ω -→ r(ω η (θ(x))u(ω η (θ(x)) -r(θ)u(θ(x)), belongs to H η .As, for η small enough, for h η in H η , one has∂Ω hη = ∂Ω, E k,hη = E k for all k ∈ {1, . . . , M }, k = m, and E m,hη = r(θ)u(θ), θ ∈ [Θ 1 m , Θ 2 m + η] , we have U (Θ 1 , Θ 2 η , I) = R(h η , I). Now, cumbersome computations show that, for all h η in H η , for all x in ∂Ω, one has h η (x) = θ(x + h η (x)) -θ(x) ρ(θ(x))τ (x) + o( h η C 1 ).Applying Theorem A.2 with d = 2, immediately gives that for all Ĩ ∈ R M and (ũ, Ũ ) corresponding solution of (2.5), one has(R (h η )I) • Ĩ = η ρ(Θ 2 m )(U m -u(x 2 m ))( Ũm -ũ(x 2 m )) + o( h η C 1 ). As U (Θ 1 , Θ 2 η , I) -U (Θ 1 , Θ 2 , I) • Ĩ -η ρ(Θ 2 m )(U m -u(x 2 m ))( Ũm -ũ(x 2 m )) = R(h η , I) -R(0, I) -R (h η )I • Ĩ + o( h η C 1 ), Theorem A.2 directly implies Lemma 4.6. Acknowledgements. This work has been partially supported by the ANR LabEx CIMI (under grant ANR-11-LABX-0040) within the French State Programme "Investissements d'Avenir". Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine.

  1, 1] × [-1, 1], discretized by a 400 × 400 Cartesian grid, N = 4 and α = [0.8, 0.02, 0.001, 0.05, 0.001, 0.04, 0.001]. The obtained results are displayed in the following tables:

		Searched positions	Start positions	Reconstructed positions
		Θ 1	Θ 2	Θ 1	Θ 2	Θ 1	Θ 2
	E 1 -2.51327 -1.93817 -3.14159 -2.51358	-2.55152 -1.97563
	E 2 -0.94247 -0.30516 -1.57079 -0.91229	-0.90298 -0.27141
	E 3	0.62831 1.24343	0	0.57253	0.54018 1.14782
	E 4	2.19911 2.84837	1.57079 2.18999	2.08957 2.72868
		Searched positions	Starting positions	Reconstructed positions
		Θ 1	Θ 2	Θ 1	Θ 2	Θ 1	Θ 2
	E 1 -2.35619 -1.78157 -3.14159 -2.51358	-2.24361 -1.66646
	E 2 -0.78539 -0.17053 -1.57079 -0.91229	-0.70862 -0.10368
	E 3	0.78539 1.41055	0	0.57253	0.83745 1.46474
	E 4	2.19911 2.84837	1.57079 2.18999	2.23475 2.88755

Appendix A.

Conductivity reconstruction -Proof of Lemma 4.4. We prove Lemma 4.4 with P = 1, that is with a single current pattern I. The extension to several current patterns is straightforward. We therefore fix I ∈ R M and U meas ∈ R M such that U meas,1 = 0.

Let (u(σ, I), U (σ, I)) be the solution of (2.5) for a given conductivity σ ∈ L ∞ ∩ H 1 (Ω) satisfying σ ≥ c > 0 a.e. in Ω. It is known that the mapping M : σ -→ (u(σ, I), U (σ, I)), is Frechet differentiable [START_REF] Kaipio | Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography[END_REF]. More precisely, for all δσ in L ∞ ∩ H 1 (Ω) such that σ + δσ ≥ c > 0 a.e. in Ω, one has

where δu ∈ H 1 (Ω) and δU ∈ R M , δU 1 = 0, are the only solution of the following variational problem:

Equivalently, (δu, δU ) satisfies

As a consequence, we have

We now define (w, W ) as the unique solution of (2.5) with f = g = 0 and input current I = U (σ, I input ) -U meas , that is

Using (A.1) and the variational formulation of problem (2.5), we obtain that

As a consequence, for δσ in L ∞ ∩ H 1 (Ω) such that σ + δσ ≥ c > 0 in Ω, and t ∈ (0, 1), we have

Our aim is now to choose δσ such that F (σ + tδσ) ≤ F (σ) for t > 0 small enough. The previous formula shows that it is sufficient to choose δσ such that

Let us define v as the unique function of H 1 0 (Ω) such that for all ṽ in H 1 0 (Ω),


Yun-Ze Long
email: yunze.long@163.com

Duvail Jean-Luc
email: duvail@cnrs-imn.fr

Chen Zhao-Jia

Jin Ai-Zi

Electrical Conductivity and Current-Voltage Characteristics of Individual Conducting Polymer PEDOT Nanowires *

We report the current-voltage (I-V ) characteristics and electrical conductivity of individual template-synthesized poly(3,4-ethylenedioxythiophene) (PEDOT) nanowires (190 ± 6 nm in diameter and σRT = 11.2 ± 2 Ω -1 cm -1 ) over a wide temperature range from 300 to 10 K. With lowering temperature, the I -V characteristics become nonlinear around 50 K, and a clear Coulomb gap-like structure appears in the differential conductance (dI/dV ) spectra. The temperature dependence of the resistance below 70 K follows ln R ∝ T -1/2 , which can be interpreted as Efros-Shklovskii hopping conduction in the presence of a Coulomb gap. In addition, the influences of measurement methods such as the applied bias voltage magnitude, the two-probe and four-probe techniques used in the resistance measurements are also reported and discussed.

Conjugated polymer nanotubes/wires have attracted much attention due to fundamental research and applications in nanoscale devices (for a review, see Refs. [1,2] and references therein). Recent studies on conducting polymer nanostructures have revealed interesting features different from bulk polymer films and pellets. [1][2][3][4][5][6][7][8][9][10] For example, it has been shown that electrical [3][4][5][6][7] and mechanical [8] properties of the polymer nanotubes/wires strongly depend on the diameter. Moreover, Coulomb-blockade transport in helical polyacetylene nanofibres, [9] small magnetoresistance of single polymer nanotubes/wires, [10] and strong electron-electron interaction arising from nanosize effect and strong disorder have been evidenced in the polymer nanostructures. [1,6,9] Among conducting polymer nanostructures, poly (3,4-ethylenedioxythiophene) (PEDOT) nanotubes/wires have been much investigated for applications in organic light-emitting diodes, field emission, [11] organic field effect transistors, [12] and electrochromic displays et al. because of its good stability in air, relatively high conductivity, high optical transparency and high efficiency in hole injection. So far, various methods have been used to prepare PEDOT nanotubes/wires, such as template synthesis, [5,7,11,[13][14][15] template-free electrochemical polymerization, [12] nanofibre seeding, microemulsion polymerization, [16] molecular combing method, [17] etc. Some studies have explored the electrical properties of single PEDOT nanowires either by focused-ion beam deposition technique, [16] or by dispersion of nanowires on pre-patterned micro-or nanoleads. [11,17] The effect of HF and NaOH treatments on the conductivity of PEDOT-DBSA nanowires was also studied. [11] Samitsu et al. [17] have measured the I -V curves and conductivity (σ RT ∼ 0.6 and 0.09 S/cm, lnG ∝ T -1/2 ) of individual PEDOT-PSS nanowires on platinum nanoelectrodes.

Although a lot of investigations have been carried out, detailed measurements of conductivity and I -V curves as a function of temperature, and effects of the experimental conditions on measurement results have not been systematically investigated yet. In this Letter, we present the temperature dependence of I -V characteristics and the electrical conductivity of template-synthesized PEDOT nanowires after removal of the membrane and dispersion on thick SiO 2 layer. A transition from linear to nonlinear I -V curves, the emergence of a Coulomb gap-like structure, and the resistance variation as ln R ∝ T -1/2 strongly suggest that electron-electron interaction is a key issue dominating the low-temperature electronic transport processes in the polymer nanowire.

The PEDOT nanowires were prepared in templates of polycarbonate track-etched membranes. [5,7,[13][14][15] In a typical synthesis procedure, we used a gold layer evaporated on one side of the membrane as the working electrode, a platinum plate as a counter electrode, and a saturated calomel electrode (SCE) as a reference. The polymerization bath consisted of an aqueous solution containing 0.07 M sodium dodecyl sulfate, 0.1 M LiClO 4 , and 0.05 M 3,4-ethylenedioxythiophene (provided by BAYER AG and distilled before using). The electrochemical polymerization was carried out at a fixed potential of 0.8 V vs. SCE with an EGG 273 potentiostat. After the polymerization, polycarbonate (membrane template) was removed by dissolution with a flow of dichloromethane and nanowires were dispersed onto a 400-nm-thick SiO 2 wafer. The resulting PEDOT nanowires have been characterized by field-emission scanning electron microscopy (SEM), transmission electron microscopy, resonant Raman spectroscopy, in situ absorption spectroscopy, x-ray photoelectron spectroscopy and electron spin resonance. More details can be found in Refs. [13-15]. To attach Pt microleads on a dispersed nanowires, a scanning electron microscope was used to find an appropriately PEDOT nanowire on the wafer. Then, four Pt microleads with typically 0.5 µm in width and 0.4 µm in thickness were fabricated by FIB deposition (Dual-Beam 235 FIB System from FEI Company, working voltage of the system is 5 kV for the electron beam and 30 kV for the focused-ion beam, respectively, current of the focused-ion beam is 1-10 pA). Electrical connection between the Pt microleads and the sample holder was made by highly conductive silver paste and gold wires. The I -V characteristics and electrical resistance of the PEDOT nanowires were measured by using a physical property measurement system from Quantum Design and a Keithley 6487 picoammeter/voltage source over a wide temperature range from 300 down to 10 K. The I -V curves were obtained by scanning the voltage from -1.0 to 1.0 V with a step of 0.005 V at different temperatures. The dI/dV curves were numerically derived from the corresponding I -V curves.

Figure 1 shows the SEM image of a PEDOT nanowire, on which four Pt microleads are attached. The diameter of the measured nanowire is about 190± 6 nm. It is noted that the leads' resistance are small and negligible compared to the nanowire's resistance because the four-probe resistance of the nanowire is larger than 60 kΩ, whereas the resistance of the Pt leads is less than 1 kΩ (estimated using the widely recognized resistivity 5 × 10 -4 Ωcm for the Pt film under the conditions used for the FIB deposition). [18][19][20][21] The room-temperature conductivity of the measured nanowire is 11.2 ± 2 Ω -1 cm -1 . Here we note that at least two individual PEDOT nanowires (190 ± 6 nm in diameter) have been measured, and the measurement results are repeatable. Previous studies on electronic transport in individual polymeric and inorganic nanotubes/wires [1,6,9,[22][23][24] have demonstrated that Coulomb interaction in nanostructures is enhanced and must be taken into account at low temperatures due to the confinement of electrons (nanosize effect) and/or strong disorder. For example, the Coulomb gap for the conducting polyaniline tube is 20-40 times larger than that for bulk samples. [6] In order to determine whether the Coulomb interaction in the PE-DOT nanowire is strong or not, we measured the I -V characteristics over a wide temperature range (300 to 10 K), as shown in Fig. 2(a). At temperatures higher than (or equal to) 50 K, the I -V curves are linear (I -V data above 60 K not shown). With lowering temperature the I -V curves gradually become nonlinear. Figure 2(b) exhibits the differential conductance dI/dV spectra taken numerically from the corresponding I -V characteristics at 50, 40, 30, 25, and 20 K. The dI/dV spectra have a pronounced minimum around zero bias at low temperatures. According to the theory developed by Pollak, Efros and Shklovskii et al., [START_REF] Pollak | Hopping Transport in Solids[END_REF][START_REF] Efros | Electron-Electron Interaction in Disordered Systems[END_REF] a Coulomb gap occurs in the density of localized states due to long-range Coulomb interaction between the electrons in different centres of localization near the Fermi energy. The depletion of the differential conductance at low energies is widely reported and usually regarded as a signature of electron-electron interaction, [START_REF] Massey | [END_REF]28] in accordance with the gradual opening of a Coulomb gap at low temperatures. Here we note that Kaiser et al. [29] have proposed a generic expression (extending the fluctuation-assisted tunnelling and thermal excitation model) for the nonlinear I -V characteristics based on numerical calculations for quasi-one-dimensional conduction interrupted by small barriers, or for granular materials. However, our calculations indicate that the I -V data around zero bias below 40 K always deviate from the Kaiser model because Coulomb gap (or Coulomb blockade) effect (zero-bias anomaly) has not been considered by this model. [START_REF] Long | Template Synthesis of PEDOT Nanowires and Electrical Properties of Individual PEDOT Nanowires[END_REF] The above results indicate that electron-electron interaction is an important issue dominating the electronic transport processes in the polymer nanowires especially at low temperatures. According to the theory of Mott variable-range hopping (VRH), [START_REF] Pollak | Hopping Transport in Solids[END_REF] if the effect of electron-electron interaction is negligible, the temperature dependence of resistance can be expressed as R(T ) = R 0 exp(T M /T ) 1/(d+1) , where T M is the characteristic Mott temperature, and d is the dimensionality of the hopping process. However, if the electron-electron interaction is strong and taken into account, the existence of the Coulomb gap will result in an Efros-Shklovskii (ES) VRH law, namely, R(T ) = R 0 exp(T ES /T ) 1/2 , where T ES is the characteristic ES temperature, and the exponent of 1/2 is independent of the dimensionality of the disordered system. Thus it is interesting to study the temperature-dependent resistance of the single PE-DOT nanowire.

However, since the I -V characteristics are nonlinear at low temperatures (as shown in Fig. 2), one question arises before discussing the temperature dependence of resistance: how is the measured electrical resistance affected by the measurement methods such as four-probe technique, two-probe technique, and applied bias voltage magnitude? Figure 3 shows the two-probe resistance data of the middle section of the nanowire measured by applying different bias voltage V bias = 0.5, 0.1, and 0.02 V. At high temperatures the resistance data measured under different V bias are nearly consistent, whereas at low temperatures the resistance values determined under V bias = 0.1 and 0.02 V are larger than that determined under V bias = 0.5 V, which is in accordance with the dI/dV spectra shown in Fig. 2. Figure 3(b) also indicates that below 70 K, ln R is linear in T -1/2 . Here we note that Samitsu et al. [17] have already observed ln G ∝ T -1/2 in individual PEDOT-PSS nanowires with diameters under 10 nm owing to quasi-one dimensional VRH conduction. For the present case, in the view point of electrons, the PEDOT nanowire is still three-dimensional because the localization length of electrons (less than 20 nm [6] ) is much smaller than the diameter (190 nm) of the nanowire. Thus we ascribe the observed ln R ∝ T -1/2 temperature dependence to Efros-Shklovskii VRH conduction model. This result is in agreement with the above analysis. The characteristic ES temperature determined from Fig. 3(b) is T ES = 6072 K (V bias = 0.1 and 0.02 V) and 4882 K (V bias = 0.5 V). It is indicated that the values of V bias used in the two-probe measurements will affect the characteristic ES temperature T ES . In order to compare the two-probe resistance with the four-probe resistance, Fig. 4 gives the temperature dependence of the four-probe resistance measured by applying a very small current (I = 0.01-10 nA and the corresponding V = 0.0005-0.02 V). The four-probe resistance data is quite close to the two-probe resistance data determined under V bias = 0.02 V. This result indicates that the microlead-nanowire nanocontact resistance is small and negligible for the present measured PEDOT nanowire. However, if the nanocontact resistance is comparable with or even larger than the intrinsic resistance of the nanowires, the three-probe or four-probe resistance will much smaller than the two-probe resistance, and only the four-probe measurements can reveal the intrinsic electronic transport properties of the nanowires. [31] For example, the temperature dependence of two-probe resistance of a nanowire may indicate that the nanowire is semiconducting, whereas the four-probe resistance dependence of the same nanowire may demonstrate that the measured wire is metallic. [31] Thus a four-probe resistance measurement is necessary and much more reliable in transport property studies of individual nanowires because nanocontact resistances may dominate the measured resistance. In summary, the low-temperature electronic transport in template-synthesized PEDOT nanowires (190±6 nm in diameter and 11.2±2 Ω -1 cm -1 at room temperature) is investigated by two-probe and fourprobe measurements. With lowering temperature, a transition from linear to nonlinear current-voltage characteristics occurs around 40 K, clear Coulomb gap-like structure gradually emerges on the differential conductance curves, and the temperature dependence of resistance (ln R ∝ T -1/2 ) follows the Efros-Shklovskii variable-range hopping conduction model. These results demonstrate that electronelectron interaction is strong and dominates the lowtemperature electronic transport processes in the PE-DOT nanowires. The two-probe and four-probe resistance measurements also indicate that the applied bias voltage in the two-probe resistance measurements will affect the characteristic ES temperature.

Fig. 1 .

 1 Fig. 1. SEM image of a template-synthesized PEDOT nanowire attached with four Pt microleads. The diameter of the measured nanowire is about 190 nm.

Fig. 2 .

 2 Fig. 2. (a) I -V characteristics and (b) the corresponding differential conductance dI/dV spectra of the single PEDOT nanowire at 50, 40, 30, 25, and 20 K.

Fig. 3 .

 3 Fig. 3. (a) Temperature dependence of the two-probe resistance of the PEDOT nanowire measured by applying different bias voltage (V bias = 0.5, 0.1, and 0.02 V), (b) plotted as ln R vs T -1/2 . The characteristic ES temperature T ES = 4882 K(V bias = 0.5 V) and 6072 K(V bias = 0.1 V).

Fig. 4 .

 4 Fig. 4. Temperature dependence of the four-probe resistance of the PEDOT nanowire. For comparison, the two-probe resistance data measured under V bias = 0.02 V is also plotted.

* Supported by the National Natural Science Foundation of China under Grant No 10604038, and the Programme for New Century Excellent Talents in University of China under Grant No NCET-07-0472.