Jan Vogel

Salia Cherifi

Stefania Pizzini

Fabien Romanens

Julio Camarero

Frédéric Petroff

Stefan Heun

Andrea Locatelli

Layer-resolved imaging of domain wall interactions in magnetic tunnel junction-like trilayers

Keywords: numbers: 75.60.Ch, 75.70.Kw, 07.85.Tt

come

Introduction

Magnetic trilayer systems consisting of two thin ferromagnetic (FM) layers separated by a non-magnetic (NM) spacer layer exhibit a wealth of interesting physical phenomena that are increasingly used in applications. Examples are the Giant Magneto resistance effect [START_REF] Baibich | Giant magnetoresistance of (001)Fe/(001)Cr magnetic superlattices[END_REF][START_REF] Binasch | Enhanced magnetoresistance in layered magnetic structures with antiferromagnetic interlayer exchange[END_REF], the difference in resistance for parallel and antiparallel alignment of the magnetization directions of the two FM layers, which is used in magnetic read heads based on spin-valves [START_REF] Dieny | Giant magnetoresistive in soft ferromagnetic multilayers[END_REF]. The inverse effect, the influence of an electrical current on the magnetic configuration of the trilayer system, has drawn a lot of attention recently [START_REF] Slonczewski | Current-driven excitation of magnetic multilayers[END_REF][START_REF] Myers | Current-induced switching of domains in magnetic multilayer devices[END_REF][START_REF] Katine | Current-driven magnetization reversal and spin-wave excitations in co/cu/co pillars[END_REF].

In these trilayer systems, interactions between the two FM layers through the NM spacer layer can be induced either by exchange [START_REF] Bruno | Ruderman-kittel theory of oscillatory interlayer exchange coupling[END_REF], magnetostatic [START_REF] Néel | Sur un nouveau mode de couplage entre les aimantations de deux films[END_REF] or dynamic effects [START_REF] Heinrich | Dynamic exchange coupling in magnetic bilayers[END_REF]. Local variations in these interlayer interactions, caused by modulated topography [START_REF] Pennec | Switching mode dependent magnetic interlayer coupling strength in spin-valves and magnetic tunnel junctions[END_REF] or inhomogeneous magnetization, can be important for the static and dynamic magnetic properties of FM/NM/FM trilayers. In this paper, we will focus on the effect of domain walls in one or both FM layers, which lead to strongly localized magnetostatic interactions between the two layers. These interactions can lead to a local decrease of tle nucleation barrier for magnetization reversal [START_REF] Vogel | Influence of domain wall interactions on nanosecond switching in magnetic tunnel junctions[END_REF] or to the demagnetization of one of the layers [START_REF] Thomas | Domain-wall induced coupling between ferromagnetic layers[END_REF]. They might also play a role in the decrease of the critical current in trilayer systems for the propagation of domain walls induced by spin-polarized current [START_REF] Grollier | Switching a spin valve back and forth by current-induced domain wall motion[END_REF][START_REF] Laribi | Reversible and irreversible current induced domain wall motion in cofeb based spin valves stripes[END_REF]. The goal of this study was to obtain a better understanding of the interaction between magnetic layers through a nonmagnetic spacer layer in the vicinity of domain walls, using high-resolution x-ray photoemission electron microscopy combined with x-ray magnetic circular dichroism (XMCD-PEEM).

Microscopic evidence of the influence of domain wall stray fields in one layer on the static domain configuration of another layer was obtained by different groups [START_REF] Fuller | Magnetostatic interactions between thin magnetic films[END_REF][START_REF] Kuch | Layerresolved imaging of magnetic interlayer coupling by domain-wall stray fields[END_REF][START_REF] Schäfer | Domain wall induced switching of whisker-based tunnel junctions[END_REF][START_REF] Schäfer | Numerical simulation of domain walls in fe whiskers and their interaction with deposited thin films[END_REF][START_REF] Wiebel | Domain decoration in dipolar coupled ferromagnetic stacks with perpendicular anisotropy[END_REF]. Thomas et al. [START_REF] Thomas | Domain-wall induced coupling between ferromagnetic layers[END_REF] showed that stray fields of domain walls sweeping back and forth during repeated switching of one of the FM layers can demagnetise the other magnetic layer, even if it is much harder magnetically. We have recently shown [START_REF] Vogel | Influence of domain wall interactions on nanosecond switching in magnetic tunnel junctions[END_REF] that in nanosecond magnetization switching, domain walls in the hard magnetic Co layer of Fe 20 Ni 80 /Al 2 O 3 /Co trilayers can trigger nucleation in the soft FeNi layer and thereby increase the local switching speed. These observations were made taking advantage of the element selectivity, spatial and temporal resolution of time-resolved XMCD-PEEM [START_REF] Kuch | X-ray magnetic circular dichroism for the quantitative element-resolved magnetic microscopy[END_REF][START_REF] Vogel | Time-resolved magnetic domain imaging by x-ray photoemission electron microscopy[END_REF][START_REF] Bonfim | Element-selective nanosecond magnetization dynamics in magnetic heterostructures[END_REF][START_REF] Schneider | Accessing fast magnetization dynamics by xpeem: Status and perspectives[END_REF]. The high-resolution static XMCD-PEEM images of FeNi/Al 2 O 3 /Co trilayers presented in this paper nicely reveal the influence of stray fields of domain walls in the Co layer on the local magnetization of the FeNi layer. These stray fields lead to a localised antiparallel coupling between the Co and FeNi moments. For strongly charged walls, which are not parallel to the easy magnetization axis of the films, this antiparallel coupling can extend over radial distances of several micrometers from the centre of the domain wall.

Experimental details

The sample consisted of a stack of continuous layers, CoO(3nm)/Co(7nm)/Al 2 O 3 (2.6nm) /Fe 20 Ni 80 (4nm), deposited on an oxidized Si(111) substrate. It was covered with 3 nm of Al to protect the layers from oxidation. The Si(111) substrate was miscut by 8 • along the [112] direction and then heat-treated to obtain a step-bunched surface, with terraces having an average width of 20 nm separated by 4 nm high steps [START_REF] Sussiau | Bicrystalline magnetic lateral superlattices[END_REF]. For sufficiently thin films, this substrate topography is transmitted to the deposited layers [START_REF] Pennec | Switching mode dependent magnetic interlayer coupling strength in spin-valves and magnetic tunnel junctions[END_REF][START_REF] Montaigne | Effet tunnel dépendant du spin : Des simples aux doubles jonctions[END_REF]. The presence of elongated terraces leads to a shape-induced uniaxial magnetic anisotropy with the easy magnetization axis along the long axis of the terraces. The correlated roughness at the two FM/NM interfaces leads to a magnetostatic orangepeel coupling [START_REF] Néel | Sur un nouveau mode de couplage entre les aimantations de deux films[END_REF] between the two magnetic layers [START_REF] Pennec | Switching mode dependent magnetic interlayer coupling strength in spin-valves and magnetic tunnel junctions[END_REF][START_REF] Encinas | Properties of spin-valve structures deposited on step-bunched vicinal surfaces[END_REF].

In order to obtain XMCD-PEEM images with high spatial resolution (< 40nm), we made use of the SPELEEM instrument at the Nanospectroscopy beamline at the synchrotron ELETTRA (Trieste, Italy) [START_REF] Locatelli | High lateral resolution spectroscopic imaging of surfaces: The undulator beamline "nanospectroscopy[END_REF][START_REF] Locatelli | Photoemission electron microscopy with chemical sensitivity: Speleem methods and applications Surf[END_REF]. In this microscope, the x-rays hit the sample surface at a grazing incidence angle of 16 • . It is therefore particularly suited for the study of samples with in-plane magnetization. Element-selective magnetic contrast was obtained using XMCD, with the x-ray energy tuned to the maximum of the Fe L 3 absorption edge (720 eV) for the FeNi layer and of the Co L 3 -edge (793 eV) for the Co layer. The images, obtained at room temperature, represent the asymmetry (difference divided by sum) of two measurements taken with opposite photon helicity. This allows optimizing the magnetic contrast while minimizing topographic effects.

Results and discussion

The quasi-static magnetization hysteresis loops of the sample, obtained using the longitudinal magneto-optical Kerr effect, are shown in Fig. 1. The quasi-static magnetization reversal of this sample takes place by nucleation of a small number of domains and a subsequent propagation of the domain walls [START_REF] Romanens | Influence of topography and Co domain walls on the magnetization reversal of the FeNi layer in FeNi/Al2O3/Co magnetic tunnel junctions[END_REF]. The minor loops of the FeNi layer reveal a coupling strength of about 0.6 mT.

Before taking images, an AC magnetic field with decreasing amplitude was applied to the sample along the easy magnetization axis, in order to create a multidomain state in both the FeNi and Co layers. A small constant field of about -3 mT was then applied, higher than the FeNi coercivity but lower than the Co coercivity, in order to saturate the FeNi layer without changing the Co domain structure (Fig. 1). The sample was then introduced into the microscope, and several regions of interest were imaged. Figure 2 shows some typical magnetic domain images of the Co (a) and µ 0 H (mT) FeNi (b) layers. The images were obtained with the projection of the x-ray incidence direction on the sample surface parallel to the easy magnetization axis. The contrast is due to the difference in absorption for domains with their magnetization direction parallel and antiparallel to the incoming x-ray direction.

The image of the Co layer shows a domain structure typical for a magnetic layer with in-plane magnetic anisotropy, with one large white domain (top left part of the image) and one large black domain. Contrast variations inside the domains are due to inhomogeneous intensity of the x-ray beam. Two main domain wall orientations are present, one parallel to the easy magnetization axis (indicated with a black arrowhead) and the other tilted with respect to this axis of about 45 • (grey arrowhead).

In the FeNi image, clear differences are observed between the regions above the two sections of the Co domain wall with different orientations. Above the parallel DW, a faint grey line is visible, indicating that the magnetization is tilted away from the easy axis to form a 'quasi-wall' [START_REF] Vogel | Influence of domain wall interactions on nanosecond switching in magnetic tunnel junctions[END_REF][START_REF] Hubert | Magnetic Domains : The Analysis of Magnetic Microstructures[END_REF]. For the layer thickness used here, domain walls in both Co and FeNi layers are expected to be of Néel-type. The presence of the Néel wall in the Co-layer implies a local divergence of the in-plane magnetization, leading to a magnetic charge -∇•M. The associated stray field points in the direction opposite to the magnetization at the centre of the Co domain wall, as indicated in Fig. 3(left), and induces the 'quasi-wall' in the FeNi layer.

In the vicinity of the tilted Co DW, much larger white domains are visible in the FeNi layer. The magnetization directions in the Co layer around this DW are partly pointing towards each other, leading to a so-called 'head-on' DW. The stray field induced by this strongly charged DW, as indicated in Fig. 3(right), is much larger than for the parallel DW. In order to partly compensate the magnetic charges on the Co DW, an oppositely charged domain wall is formed in the FeNi layer (a tail-totail domain wall). The Co DW stray field therefore locally induces an antiparallel coupling between the two magnetic layers through the alumina spacer layer, as is clearly seen in Fig. 2: the white domains in the FeNi layer are situated on the right side of the Co DW, above the black domain in the Co layer. This DW induced antiparallel coupling apparently is much stronger than the orange-peel coupling that favors parallel alignment.

In order to allow a better observation of the domain walls, we rotated the sample by 90 • . The corresponding images, thus taken with the x-rays coming in along the hard magnetization axis, are shown in Fig. 4. The azimuthal rotation was performed outside the microscope, since an in-situ rotation was not possible. This makes it almost impossible to find back exactly the same region of the sample in this kind of continuous films. The region in Fig. 4 shows therefore a domain structure that is similar to the one in Fig. 2, but not exactly the same. In this case, the observed black and white regions correspond to the domain walls or quasi-walls, while the domains show intermediate grey contrast. In the left (Co) image, a zigzag domain wall is observed as in Fig. 2, and the magnetization direction at the centre of the DW is indicated by the arrows. The expected magnetization direction inside the domains is also given. Two Bloch lines (vortices), characterised by a reversal of the direction of magnetization at the centre of the DW, can be clearly observed. They are indicated by black arrowheads in Fig. 4a).

Figure 4 b) shows that at positions where a domain wall is present in the Co layer, domain walls or quasi-walls are present also in the FeNi layer. The contrast at these positions and thus the magnetization direction at the centre of the wall is opposite in the Co and FeNi layers. This is due to the direction of the domain wall stray field, and allows a partial compensation of the domain wall charges in both layers. Micromagnetic simulations confirm this picture [START_REF] Vogel | Influence of domain wall interactions on nanosecond switching in magnetic tunnel junctions[END_REF]. Additional domain walls in the

Parallel

Head-on FeNi layer, which are not associated to the Co DW, are due to the extended domains generated in the FeNi layer by the stray field of the tilted Co DW (see Fig. 2b).

Outside the region where domain walls are present, contrast modulations are visible. These modulations are not due to noise, since they are reproducible between different series of images taken at the same position. This indicates that the local magnetization direction can deviate from the global easy magnetization axis. This is a remainder of the so-called 'ripple' structure often observed for uniaxial polycrystalline thin films when a magnetic field is applied perpendicular to the easy magnetization axis. The phenomenon can be associated to the crystalline anisotropy of individual grains [START_REF] Harte | Theory of magnetization ripple in ferromagnetic films[END_REF]. In our case, it is probably associated to a modulation of the long axis of the substrate terraces with respect to the easy [110] axis. The long axis of each terrace determines the local shape anisotropy, leading to a distribution of magnetization directions around the average, macroscopic easy magnetization axis.

In Fig. 5, we show an ellipse-shaped domain structure in the Co layer, obtained after complete demagnetization of the sample. In this case, the magnetization directions in both layers are expected to be globally parallel far from the domain wall regions, due to the orange-peel coupling. The expected magnetization directions are indicated by arrows in Fig. 5. In the Co domain wall, several Bloch lines are clearly visible. At these Bloch lines, the magnetization direction at the centre of the DW reverses direction, and therefore also the emitted stray field. In most cases, the corresponding magnetization direction at the centre of the FeNi DW also changes sign. This can be seen for instance at the position indicated by the white arrowhead in Fig. 5, where the small black segment in the Co DW is accompanied by a small white segment in the FeNi DW. At the centre of the Bloch lines, the magnetization direction is perpendicular to the plane of the film [START_REF] Wachowiak | Direct observation of internal spin structure of magnetic vortex cores[END_REF], leading to the presence of additional stray field components. This leads to a more complex magnetization configuration in the FeNi film close to the Bloch lines.

The strongly charged domain walls at the extremities of the ellipse-shaped Co domain induce small domains in the FeNi layer at both sides of the Co DW. As mentioned above, the magnetic charges in the Co layer can be partly compensated by creating charges of opposite sign in the FeNi layer. This leads to a strong antiparallel coupling between Co and FeNi layers close to the charged Co domain walls. This influence of the domain wall charges is felt up to lateral distances of several micrometers from the domain wall centre.

At several positions, a grey/black/white/grey contrast can be seen around domain walls in the FeNi layer. This is what is expected for a 360 • domain wall, where the magnetization vector makes a complete 360 • turn going from one side of the domain wall to the other. However, also other configurations can give rise to this contrast. As an example consider the region in the circle of the FeNi image of Fig. 5. On the lower left, a black/white contrast is present. Going to the upper right, the black and white parts split. The white part is not associated to a DW in the Co layer, indicating that it has to be a normal 180 • Néel wall. The black part separates two domains with the same magnetization direction and is therefore a 'quasi-wall', induced by the stray field of the Co domain wall. This means that, crossing the black/white region from lower right to upper left, the magnetization vector turns over an angle of about 270 • , before turning back -90 • to end up at an angle of 180 • with respect to the initial magnetization direction. The contrast is thus induced by the combination of a Néel wall and a 'quasi-wall' with the same chirality. A zoom of the magnetic configuration is schematically shown in the circle.

These images can be used to obtain an estimate of the domain wall width in both layers. The definition we use for determining the domain wall width is based on the total wall flux and is given by W F = ∞ -∞ cos ϕ(x)dx [START_REF] Hubert | Magnetic Domains : The Analysis of Magnetic Microstructures[END_REF], where ϕ(x) is the angle between the local magnetization direction and the hard magnetization axis. This definition is the most convenient in our case, since the XMCD-intensity in the images of Fig. 4 is proportional to the projection of the local magnetization direction on the xray incidence direction and thus to cos ϕ. The wall width can be obtained by taking an image profile along a direction perpendicular to the DW. Using the highest resolution images (field of view 10 µm), we find widths of 130 ± 20 nm for the Co domain walls and 170 ± 15 nm for the FeNi domain walls. However, the domain wall widths seem to be quite irregular, probably due to the specific, irregular topography of terraces and steps that can locally lead to confinement or, on the contrary, largening of the domain walls [START_REF] Pennec | Switching mode dependent magnetic interlayer coupling strength in spin-valves and magnetic tunnel junctions[END_REF]. The DW width of Co is smaller than for FeNi. According to Ref. [START_REF] Hubert | Magnetic Domains : The Analysis of Magnetic Microstructures[END_REF], Chap. 5.5.7, the Néel wall width in each layer of a magnetic/non-magnetic/magnetic trilayer system is determined by the magnetic anisotropy K u and the stray field energy (with parameter K D = J 2 S /µ 0 , where J S is the spontaneous magnetization). While the magnetic anisotropy tends to decrease the domain wall width, the stray field energy tends to increase it. In our case, the magnetic anisotropy is larger for the Co layer (K Co U ≈ 11200J/m 3 against K F eN i U ≈ 1720J/m 3 (Ref. [START_REF] Vogel | Influence of domain wall interactions on nanosecond switching in magnetic tunnel junctions[END_REF]) but also the exchange and the stray field energy (J Co S = 1.76T, J F eN i S = 1T ). The wall width of the Co domain walls is not much smaller than the one of the FeNi layer, but the stray field energy is about three times higher. Moreover, it is easier to deviate the FeNi magnetization from its easy axis due to the smaller anistropy. The influence of the Co DW stray field on the FeNi layer is therefore much larger than the reciprocal effect, which is however still visible in some images (see Fig. 4a, for example).

Co

FeNi a) b) Time-resolved XMCD-PEEM measurements on this same sample already revealed the strong influence of the Co domain walls on the magnetization reversal of the FeNi layer (Fig. 7 of Ref. [START_REF] Romanens | Influence of topography and Co domain walls on the magnetization reversal of the FeNi layer in FeNi/Al2O3/Co magnetic tunnel junctions[END_REF]). In that case, a non-saturated region was observed in the FeNi layer just above the Co domain wall even for fields of 6 mT. This region seemed to be rather large (about 1 µm), larger than we observe in these static images, even taking into account the difference in spatial resolution (about 50 nm here against 300 nm in Ref. [START_REF] Romanens | Influence of topography and Co domain walls on the magnetization reversal of the FeNi layer in FeNi/Al2O3/Co magnetic tunnel junctions[END_REF]). In order to investigate the influence of repeated magnetic pulses on the domain walls, we recorded images of the Co and FeNi layers after application of some thousands of bipolar magnetic pulses of about 5 mT, using the same coil and pulsed current supply as in Refs. [START_REF] Romanens | Influence of topography and Co domain walls on the magnetization reversal of the FeNi layer in FeNi/Al2O3/Co magnetic tunnel junctions[END_REF][START_REF] Vogel | Time and layer resolved magnetic domain imaging of FeNi/Cu/Co trilayers using x-ray photoelectron emission microscopy[END_REF]. The corresponding images are shown in Fig. 6.

The images were recorded with the x-rays coming in at an azimuthal angle of 80 • with respect to the easy magnetization axis. In this way it is possible to reveal the magnetization direction in both the domains and the domain walls. The FeNi image clearly shows signs of domain accumulation around the Co domain wall over distances well above 1 µm. Compression of these domains during the application of the field pulses might lead to 360 • domain walls and be at the origin of the larger grey region observed in the images of Ref. [START_REF] Romanens | Influence of topography and Co domain walls on the magnetization reversal of the FeNi layer in FeNi/Al2O3/Co magnetic tunnel junctions[END_REF].

Conclusion

We have performed layer-selective imaging of magnetic domains in FeNi/ Al 2 O 3 /Co trilayers with high spatial resolution. These images reveal the strong interaction between domain walls in the two magnetic layers. Co domain walls parallel to the easy magnetization axis induce 'quasi-walls' in the FeNi layer, with the magnetization at the centres of the two walls in opposite directions. For domain walls that are tilted with respect to the easy magnetization axis magnetostatic interactions favor an anti-parallel alignment of the magnetization directions in both magnetic layers over a distance up to several microns from the domain wall. Upon repetitive application of magnetic field pulses, these interactions can lead to an accumulation of domains and domain walls in the FeNi layer around the Co DW positions.

Figure 1 .

 1 Figure1. Quasi-static hysteresis loops of the magnetic tunnel junction-like trilayer obtained by longitudinal Kerr effect. Loops obtained with the field applied along the easy (black line) and hard (grey line) magnetization axes are shown, as well as minor loops for the FeNi layer taken along the easy axis.

Figure 2 .

 2 Figure 2. Magnetic domain images of the Co (a) and FeNi (b) layers, after demagnetising the sample and consecutively re-saturating the FeNi layer. The directions of local magnetization and of incoming photons are indicated with arrows. The black and grey arrowheads next to the Co image indicate, respectively, the Co domain wall oriented parallel to the easy magnetization axis and the one tilted of about 45 • with respect to this axis

Figure 3 .

 3 Figure 3. Schematic representation of domain wall configurations for 'parallel' (left) and 'head-on' (right) Néel-type domain walls and the associated domain wall stray fields. The black arrows indicate the local magnetization direction in the film, the grey arrows the associated stray field directions.

Figure 4 .Figure 5 .

 45 Figure 4. Magnetic domain images of the Co (a) and FeNi (b) layers, taken with the incoming x-rays perpendicular to the easy magnetization axis. The black and white regions correspond to domain walls, the grey regions to domains. The black arrowheads in the Co image indicate the positions of two Bloch lines.

Figure 6 .

 6 Figure 6. Magnetic domain images of the Co (a) and FeNi (b) layers after application of several bipolar pulses with an amplitude of 5 mT. The images were taken with the x-rays coming in at an azimuthal angle of 10 • with respect to the hard magnetization axis.

Acknowledgments

We acknowledge partial funding by the European Union under contract No. HPRN-CT-2000-00134 and by a Programme d'Actions Intégrées 'Picasso' (J.V. and J.C.), through Grant No. HF2003-0173. We thank A. Vaurès for her invaluable help in sample preparation.


Habib Ammari
email: habib.ammari@polytechnique.fr

Hyeonbae Kang
email: hbkang@snu.ac.kr

Eunjoo Kim

Kaouthar Louati
email: louati@cmapx.polytechnique.fr

Michael S Vogelius
email: gelius@math.rutgers.edu

A MUSIC-type Algorithm for Detecting Internal Corrosion from Electrostatic Boundary Measurements *

Keywords: Mathematics subject classification (MSC2000): 35R30 corrosion, asymptotic representation formula, reconstruction, layer potentials

We establish an asymptotic representation formula for the steady state current perturbations caused by internal corrosive boundary parts of small surface measure. Based on this formula we design a non-iterative method of MUSIC (multiple signal classification) type for localizing the corrosive parts from voltage-to-current observations. We perform numerical experiments to test the viability of the algorithm and the results clearly demonstrate that the algorithm works well even in the presence of relatively high noise ratios.

Introduction

Hidden corrosion, for instance in pipes, is a serious problem that is responsible for significant economic losses and represents a dangerous threat to safety. In the field of nondestructive evaluation, new and improved techniques are therefore constantly being sought to facilitate the detection of such hidden corrosion.

Corrosion occurs in many different forms, and several mathematical models for electrostatic fields in the presence of corrosion have been studied in the literature (see, for example, Kaup and Santosa [START_REF] Kaup | Nondestructive evaluation of corrosion damage using electrostatic measurements[END_REF], Kaup et al. [START_REF] Kaup | A method for imaging corrosion damage in thin plates from electrostatic data[END_REF], Vogelius and Xu [START_REF] Vogelius | A nonlinear elliptic boundary value problem related to corrosion modelling[END_REF], Inglese [START_REF] Inglese | An inverse problem in corrosion detection[END_REF], Luong and Santosa [START_REF] Luong | Quantitative imaging of corrosion inplates by eddy current methods[END_REF], Banks et al. [START_REF] Banks | Real time computational algorithms for eddy-current-based damage detection[END_REF] and references therein).

In this paper, we adopt the potential model: Laplace's equation holds in the metal pipe, and the effect of corrosion is described by means of a linear Robin boundary condition. A model like this is motivated by a number of observations. On the one hand, corrosion tends to roughen a surface: this effect can be modelled by the introduction of a thin coating characterized by rapid oscillations. In the limit as the thickness of the coating goes to zero and the rapidity of the oscillations diverges, the emergence of a Robin boundary conditions has been shown by Buttazo and Kohn [START_REF] Buttazzo | Reinforcement by a thin layer with oscillating thickness[END_REF]. On the other hand, the study of active electrochemical corrosion processes can be based on "first principles", which conjecture a certain relation between the number of "free electrons" and the energy level (typically an exponential Boltzmann-type distribution). Potential models with exponential boundary conditions that arise from such considerations are often associated with the names of Butler and Volmer. Vogelius and Xu [START_REF] Vogelius | A nonlinear elliptic boundary value problem related to corrosion modelling[END_REF] study a potential model of this kind. If we linearize the boundary condition of this model, we get the linear Robin boundary condition studied here.

To set up the problem mathematically, consider a simply connected, bounded C 2,α domain U in R 2 for some 0 < α ≤ 1, and a simply connected C 2,α domain D, compactly contained in U . Ω = U \ D represents the specimen to be inspected, e.g., the cross section of the pipe. We define Γ e = ∂U and Γ i = ∂D, so that ∂Ω = Γ i ∪ Γ e . Suppose that the inaccessible boundary Γ i contains some corrosive, connected parts I s , s = 1, . . . , m. The open curves I s are well-separated, i.e., there is a constant c 0 > 0 such that dist (I i , I j ) > c 0 for i = j, and the corrosion coefficient, γ s ∈ C 1 (I s ), of each I s , s = 1, . . . , m, satisfies 0 ≤ γ s ≤ C 0 , for some constant C 0 , and is not identically zero. Let

γ(x) = m s=1 γ s χ s (x), x ∈ Γ i , (1.1)
where χ s denotes the characteristic function of I s . A typical shape of Ω would be an annulus. We assume that the curve length of each I s is small, to be specific,

d 0 ≤ |I s | Γi ≤ D 0 , s = 1, . . . , m , (1.2)
where is a small parameter, representing the common order of magnitude of I s , and d 0 , D 0 are positive constants. Each curve I s is centered around a fixed point z s :

I s = {z ∈ Γ i : d Γi (z, z s ) < |I s |/2} . (1.3)
Here | • | Γi and d Γi (•, •) refer to the one-dimensional curve length on Γ i . The internal voltage potential u , generated by a voltage f applied to Γ e , satisfies

           ∆u = 0 in Ω , - ∂u ∂ν + γu = 0 on Γ i , u = f on Γ e , (1.4)
where ν is the inward normal to Ω on Γ i (i.e., the outward normal to D).

A realistic inverse problem in corrosion detection consists of the determination of the corrosion damage on the inaccessible boundary part, Γ i , when the available data are a finite number of voltage-to-current pairs on the accessible part Γ e . The difficulties of this inverse problem stem from its inherent ill-posedness and nonlinearity. Most of the techniques for detecting the corrosion are based on iterative algorithms: least-square-fit algorithms and Newton-type iteration schemes. These methods typically incur tremendous computational costs to get close to the true corrosion damage, but more troubling: they frequently get "stuck" near local minima.

The purpose of this work is to design a direct (non-iterative) technique for electrostatic corrosion detection. This algorithm is of MUSIC-type (multiple signal classification) and is based on an accurate asymptotic representation formula for the steady state boundary currents. In many ways it is closely related to the algorithm developed in [START_REF] Brühl | A direct impedance tomography algorithm for locating small inhomogeneities[END_REF] for the purpose of detection of small internal defects.

The Cauchy data continuation technique, as studied in the paper by Yang, Choulli, and Cheng [START_REF] Yang | An iterative BEM for the inverse problem of detecting corrosion in a pipe[END_REF] for thermal detection of corrosion, is also a direct method. For reasons of instability of the solution of the Cauchy problem it is less accurate than the method described here. Furthermore it makes use of the formula γ = (1/u ) ∂u /∂ν which, in the electrostatic case is somewhat "unfortunate", due to potential sign changes (and corresponding zeroes) of u . A significant ("signal-to noise") advantage of the present method is that it directly uses the boundary current perturbation data ∂u ∂ν -∂u0 ∂ν | Γe (as opposed to ∂u ∂ν | Γe ). Here u 0 is the voltage potential in the absence of corrosion, i.e., the solution to

           ∆u 0 = 0 in Ω , - ∂u 0 ∂ν = 0 on Γ i , u 0 = f on Γ e .
(1.5)

The present paper is organized as follows: In Section 2 we review some basic facts about layer potentials. In the following section we establish a representation formula for the unique solution to (1.4). This formula generalizes the formula proved by Kang and Seo in [START_REF] Kang | Layer potential technique for the inverse conductivity problem[END_REF][START_REF]Recent progress in the inverse conductivity problem with single measurement, in Inverse Problems and Related Fields[END_REF]. The main goal of Section 4 is to rigorously derive an asymptotic expansion for ∂u ∂ν -∂u0 ∂ν on Γ e . Section 5 is then devoted to the design of the numerical algorithm, and several computational experiments with it.

Layer Potentials

In this section we review some well-known properties of layer potentials for the Laplacian, and establish some useful identities. Let Φ(x) be the fundamental solution associated with the Laplacian ∆,

Φ(x) = 1 2π ln |x| , x = 0 .
Let D be a simply connected, bounded C 2,α domain in R 2 , and let Γ :

= ∂D. Let H 1 (D) denote the set of functions f ∈ L 2 (D) such that ∇f ∈ L 2 (D). Furthermore, we define H 2 (D) as the space of functions u ∈ H 1 (D) such that ∂ 2 u ∈ L 2 (D)
, and the space

H 3/2 (D) as the interpolation space [H 1 (D), H 2 (D)] 1/2 . Let H 1/2 (∂D) be the set of traces of functions in H 1 (D) and H -1/2 (∂D) = (H 1/2 (∂D)) * . Finally, let H 1 (∂D) denote the set of functions f ∈ L 2 (∂D) such that ∂f /∂τ ∈ L 2 (∂D)
, where ∂/∂τ is the tangential derivative.

The single-layer and double-layer potentials of a function ϕ ∈ L 2 (Γ), denoted S Γ [ϕ] and D Γ [ϕ], respectively, are defined by

S Γ [ϕ](x) := Γ Φ(x -y)ϕ(y) dσ y , x ∈ R 2 ,
(2.1)

D Γ [ϕ](x) := Γ ∂ ∂ν y Φ(x -y)ϕ(y) dσ y , x ∈ R 2 \ Γ . (2.2)
For a function u, defined in R 2 \ Γ, we denote

∂u ∂ν ± (x) := lim t→0 + ∇u(x ± tν x ), ν x , x ∈ Γ ,
if the limit exists. Here ν x is the outward unit normal to Γ = ∂D at x, and , denotes the scalar product in R 2 . u| ± are defined likewise.

It is well-known (see for example [START_REF] Folland | Introduction to Partial Differential Equations[END_REF]) that for ϕ ∈ L 2 (Γ)

∂(S Γ [ϕ]) ∂ν ± (x) = ± 1 2 I + K * Γ [ϕ](x) , x ∈ Γ , (2.3)
D Γ [ϕ] ± (x) = ∓ 1 2 I + K Γ [ϕ](x) , x ∈ Γ , (2.4)
where K Γ is defined by

K Γ [ϕ](x) = 1 2π Γ y -x, ν y |x -y| 2 ϕ(y) dσ y ,
and

K * Γ is the L 2 -adjoint of K Γ , i.e., K * Γ [ϕ](x) = 1 2π Γ x -y, ν x |x -y| 2 ϕ(y) dσ y .
Since Γ is a C 2,α curve, there is a constant C such that

| x -y, ν x | |x -y| 2 ≤ C , x, y ∈ Γ . (2.5)
Thus it immediately follows that K Γ is bounded on L p (Γ) for 1 ≤ p ≤ ∞. Moreover, the operator ϕ → ∇ T K Γ [ϕ], where ∇ T denotes the tangential derivative on Γ, is a Calderón-Zygmund operator on Γ which is known to be bounded on L p (Γ) (1 < p < ∞), C β (Γ) and C 1,β (Γ) (0 < β < 1) (see [START_REF] Stein | Singular Integrals and Differentiability Properties of Functions[END_REF]). Therefore K Γ and

K * Γ map continuously L p (Γ) into H 1,p (Γ), C β (Γ) into C 1,β (Γ) (0 < β < 1), and C 1,β (Γ) into C 2,β (Γ) (0 < β ≤ α) since Γ is C 2,α .
Here H 1,p (Γ) is the collection of functions in L p (Γ) whose derivatives also belong to L p (Γ).

Lemma 2.1 The operators 1 2 I + K Γ and 1 2 I + K * Γ are invertible on L p (Γ), 1 < p, as well as on C k (Γ) for k = 1, 2.

Proofs of Lemma 2.1 on the spaces L 2 (Γ) and C(Γ) can be found in [START_REF] Folland | Introduction to Partial Differential Equations[END_REF]. The invertibility on the other spaces can seen as follows: consider the equation

1 2 I + K Γ [u] = f or 1 2 I + K * Γ [u] = f.
If f ∈ L p (Γ) for p > 2 and u ∈ L 2 (Γ), then K Γ u ∈ L p (Γ) and hence u ∈ L p (Γ). This shows that 1 2 I + K Γ and 1 2 I + K * Γ are invertible on L p (Γ) for p > 2. Invertibility on L p (Γ) for 1 < p < 2 can be seen by duality. Similarly one can prove invertibility on C k (Γ) for k = 1, 2.

In the case where Γ is a circle of radius r

ν x , x -y |x -y| 2 = 1 2r , x, y ∈ Γ ,
and hence

K Γ [ϕ](x) = K * Γ [ϕ](x) = 1 4πr Γ ϕ(y) dσ y . (2.6)
It now follows from (2.3) and (2.6) that

S Γ [1](x) =    r ln r if |x| ≤ r , r ln |x| if |x| ≥ r , (2.7)
in the case where Γ is a circle of radius r. Finally, we recall some boundedness properties of single layer potentials for later use. If p > 1, there is a constant C p such that

S Γ [ϕ] L ∞ (Γ) ≤ C p ϕ L p (Γ) (2.8)
for all ϕ ∈ L p (Γ). In fact, (2.8) can be proved using Hölder's inequality. It is also well-known that

S Γ [ϕ] H 1 (D) ≤ C ϕ L 2 (Γ) (2.9)
for any ϕ ∈ L 2 (Γ).

Representation Formula

Recall that Ω = U \ D, where U and D are simply connected, bounded C 2 domains. Let Γ e = ∂U and Γ i = ∂D. For f ∈ H 1 (Γ e ), let u 0 be the solution in the absence of corrosion, i.e., the solution to the problem (1.5). We seek to represent u 0 in the following form

u 0 = D Γe [ϕ 0 ] + S Γi [ψ 0 ] in Ω
for some functions ϕ 0 ∈ H 1 (Γ e ) and ψ 0 ∈ L 2 (Γ i ). The boundary conditions in (1.5) are, due to the relations (2.4) and (2.3), equivalent to    

1 2 I + K Γe S Γi - ∂ ∂ν i D Γe - 1 2 I -K * Γi     ϕ 0 ψ 0 = f 0 ∈ H 1 (Γ e ) × L 2 (Γ i ) , (3.1)
where ν e and ν i denote outward normal ν to Ω on Γ e and Γ i (the subscripts e and i emphasize that they are defined on Γ e and Γ i , respectively). For example,

∂ ∂νi D Γe [ϕ] = ν i • ∇D Γe [ϕ].
Lemma 3.1 For p > 1, let X p := L p (Γ e ) × L p (Γ i ) and

A 0 :=     1 2 I + K Γe S Γi - ∂ ∂ν i D Γe - 1 2 I -K * Γi     . (3.2)
Then A 0 is invertible on X p , as well as on H 1 (Γ e ) × L 2 (Γ i ).

Proof.

Since there is a distance between Γ e and Γ i , the operator A 0 is a compact perturbation of   

1 2 I + K Γe 0 0 - 1 2 I -K * Γi    ,
which is known to be invertible on X p (Lemma 2.1). Therefore, it suffices, by applying the Fredholm alternative, to show that the operator A 0 is injective. Suppose that (ϕ, ψ) ∈ X p satisfies

A 0 ϕ ψ = 0 . Since ( 1 2 I + K Γe )[ϕ] = -S Γi [ψ] on Γ e and S Γi [ψ] is C 2 on Γ e , we obtain that ϕ ∈ C 2 (Γ e ). Likewise it follows that ψ ∈ C 2 (Γ i ). Therefore, the function u defined by u = D Γe [ϕ]+S Γi [ψ]
in Ω is a classical solution to (1.5), with f = 0. Since this solution is unique, we get

D Γe [ϕ] + S Γi [ψ] = 0 in Ω , and hence D Γe [ϕ] + S Γi [ψ] = 0 in U .
The "continuity" of the normal derivative of a double-layer potential, in combination with (2.3) now gives

ψ = ∂ ∂ν S Γi [ψ] + - ∂ ∂ν S Γi [ψ] - = ∂ ∂ν (D Γe [ϕ] + S Γi [ψ]) + - ∂ ∂ν (D Γe [ϕ] + S Γi [ψ]) - = 0 on Γ i .
In other words ψ = 0, and D Γe [ϕ] = 0 in U . It follows from (2.4) that 1 2 I + K Γe [ϕ] = 0, and the invertibility of

1 2 I + K Γe now insures that ϕ = 0. Invertibility of A 0 on H 1 (Γ e ) × L 2 (Γ i ) follows easily by the fact that if (ϕ, ψ) ∈ L 2 (Γ e ) × L 2 (Γ i ) and A 0 ϕ ψ ∈ H 1 (Γ e ) × L 2 (Γ i ), then K Γe [ϕ] ∈ H 1 (Γ e )
, and hence ϕ ∈ H 1 (Γ e ). This completes the proof.

As an immediate consequence of Lemma 3.1, we obtain the following theorem.

Theorem 3.2

The solution u 0 to (1.5) can be represented as

u 0 = D Γe [ϕ 0 ] + S Γi [ψ 0 ] in Ω ,
where (ϕ 0 , ψ 0 ) is the unique solution to (3.1).

Using the exact same techniques we may also derive a representation formula for the unique solution to the boundary value problem

           ∆v 0 = 0 in Ω , - ∂v 0 ∂ν = g on Γ i , v 0 = 0 on Γ e , (3.3)
where g ∈ L 2 (Γ i ). The formula simply reads

v 0 = D Γe [ φ0 ] + S Γi [ ψ0 ] in Ω ,
where ( φ0 , ψ0 ) is the solution to

A 0 φ0 ψ0 = 0 g .
We proceed to derive a representation formula for u , the solution to (1.4), valid for sufficiently small. In this case the relevant integral operator is

A γ :=     1 2 I + K Γe S Γi - ∂ ∂ν i D Γe + γD Γe - 1 2 I -K * Γi + γS Γi     .
Observe that

A γ = A 0 + M γ B , (3.4)
where M γ represents multiplication by γ, and

B = 0 0 D Γe S Γi . Since γ L p (Γi) ≤ C 1/p ,
we have

M γ B ϕ ψ Xp ≤ γ L p (Γi) D Γe [ϕ] + S Γi [ψ] L ∞ (Γi) ≤ C p 1/p ϕ ψ Xp , for any p > 1. Here (ϕ, ψ) Xp := ϕ L p (Γe) + ψ L p (Γi)
and C p is a constant depending only on p. In other words, we have a bound for the operator norm,

M γ BA -1 0 p , of M γ BA -1 0 on X p : M γ BA -1 0 p ≤ C p 1/p . (3.5)
Likewise, one can show that

M γ BA -1 0 ϕ ψ H 1 (Γe)×L 2 (Γi) ≤ C 2 1/2 ϕ ψ H 1 (Γe)×L 2 (Γi) (3.6)
Since A γ = (I + M γ BA -1 0 )A 0 , we arrive at the following theorem. Theorem 3.3 There exists a positive number 0 such that for any 0 < < 0 A γ is invertible on X p and on H 1 (Γ e ) × L 2 (Γ i ). Moreover, for any f ∈ H 1 (Γ e ), the solution u to (1.4) can be represented as

u = D Γe [ϕ ] + S Γi [ψ ] in Ω , where (ϕ , ψ ) ∈ H 1 (Γ e ) × L 2 (Γ i ) is the unique solution to A γ ϕ ψ = f 0 . (3.7)
The last result in this section gives a complete expansion for u in terms of iterates easily constructed from u 0 . We can only prove convergence of this expansion for sufficiently small. First note that

A -1 γ = A -1 0 (I + M γ BA -1 0 ) -1 , (3.8)
and

M γ B ϕ ψ = 0 γ (D Γe [ϕ] + S Γi [ψ]) . (3.9)
By expanding (3.8) in a Neumann series, we get

A -1 γ = A -1 0 + +∞ n=1 (-1) n A -1 0 M γ BA -1 0 n . Now let ϕ n ψ n := A -1 0 M γ BA -1 0 n f 0 , n = 0, 1, 2, . . . , (3.10)
then ϕ ψ = A -1 γ f 0 = +∞ n=0 (-1) n ϕ n ψ n . (3.11)
Moreover, the following recursive relation holds:

ϕ n+1 ψ n+1 = A -1 0 M γ B ϕ n ψ n , n = 0, 1, 2, . . . . (3.12)
It then follows from (3.6) and (3.10) that

ϕ n H 1 (Γe) + ψ n L 2 (Γi) ≤ C 0 C n n/2 f H 1 (Γe) , n = 0, 1, 2, . . . , (3.13)
for some constant C independent of . Here C 0 is a bound for A -1 0 .

If we define

u n := D Γe [ϕ n ] + S Γi [ψ n ] in Ω , n = 0, 1, 2, . . . , (3.14)
which as far as u 0 is concerned is consistent with being the solution to (1.5), then the identity (3.11) immediately yields

u = +∞ n=0 (-1) n u n in Ω .
It is not difficult to see that this series converges in H 3/2 (Ω) (for sufficiently small). In fact, thanks to elliptic estimates, we have

u n H 3/2 (Ω) ≤ C ( D Γe [ϕ n ]| -H 1 (Γe) + S Γi [ψ n ] H 1 (Γi) ) ≤ C ( ( 1 2 I + K Γe )[ϕ n ] H 1 (Γe) + S Γi [ψ n ] H 1 (Γi) ) ≤ C 1 ( ϕ n H 1 (Γe) + ψ n L 2 (Γi) )
for some constant C and C 1 , and hence, by (3.13),

u n H 3/2 (Ω) ≤ C 0 C 1 C n n/2 f H 1 (Γe) . (3.15)
Thus the series for u converges strongly in H 3/2 (Ω) as long as is sufficiently small. By (3.9) and (3.12),

ϕ n ψ n = A -1 0 0 γ (D Γe [ϕ n-1 ] + S Γi [ψ n-1 ]) = A -1 0 0 γu n-1 | Γi , n ≥ 1 ,
and so u n , n ≥ 1, is the solution to            ∆u n = 0 in Ω , - ∂u n ∂ν = γu n-1 on Γ i , u n = 0 on Γ e .
(3.16)

Here we have used the observation immediately following Theorem 3.2. In summary we have proven:

Theorem 3.4 There exists a positive number 0 such that, for any 0 < < 0 , and any f ∈ H 1 (Γ e ), the solution u to (1.4) has the expansion

u = +∞ n=0 (-1) n u n in H 3/2 (Ω) , (3.17)
where u 0 is the solution to (1.5) and u n , n ≥ 1, are the solutions to (3.16).

Asymptotic Expansion

We now derive a formula for the principal term of ∂ ∂ν (u -u 0 ) (on Γ e ) as → 0. This formula forms the basis for our corrosion detection algorithm. In order to derive our formula we need an elliptic regularity estimate for the solution to the boundary value problem

           ∆v = 0 in Ω , - ∂v ∂ν = γg on Γ i , v = 0 on Γ e . (4.1)
In particular we need that

∂v ∂ν C 1 (Γe) ≤ C p γg L p (Γi) , (4.2)
for any p > 1. This regularity estimate is a consequence of well-known Schauder estimates, and could just be stated without proof -however, for completeness we show how it also immediately follows from our integral representation of the solution. We have seen that the solution v to (4.1) can be represented as

v = D Γe [ϕ] + S Γi [ψ] in Ω ,
where (ϕ, ψ) is the solution to

A 0 ϕ ψ = 0 γg , i.e. ,          1 2 I + K Γe [ϕ] + S Γi [ψ] = 0 on Γ e , - ∂ ∂ν D Γe [ϕ] - 1 2 I + K * Γi [ψ] = γg on Γ i . (4.3)
Since A 0 is invertible on X p , we have

ψ L p (Γi) ≤ C p γg L p (Γi) . (4.4)
Let 2δ = dist(Γ i , Γ e ), and define Ω δ := {x ∈ Ω : dist(x, Γ e ) ≤ δ}. Then

S Γi [ψ] C 2,α (Ω δ ) ≤ C ψ L p (Γi) ≤ C p γg L p (Γi) ,
where constants depend on δ. It now follows from the first equation in (4.3) and Lemma 2.1 that

ϕ C 2,α (Γe) = 1 2 I + K Γe -1 [(S Γi [ψ])| Γi ] C 2,α (Γe) ≤ C S Γi [ψ] C 2,α (Γe) ≤ C p γg L p (Γi) . (4.5) Note that if ϕ ∈ C 2,α (Γ e ), then D Γe [ϕ] ∈ C 2,α (U ). In fact, if ϕ ∈ C 2,α (Γ e ), then K Γe [ϕ] ∈ C 2,α (Γ e ). Since u := D Γe [ϕ] is the unique solution to the Dirichlet problem ∆u = 0 in U and u| Γe = 1 2 ϕ + K Γe [ϕ], we have D Γe [ϕ] ∈ C 2,α (U ). Therefore, we have in particular ∂v ∂ν C 1 (Γe) ≤ D Γe [ϕ] + S Γi [ψ] C 2 (Ω δ ) ≤ C p γg L p (Γi) ,
which is exactly the estimate (4.2). By a combination of (4.2) with the estimate

γg L p (Γi) ≤ C 1/p g L ∞ (Γi) ,
we obtain the following lemma.

Lemma 4.1 Let v be the solution to (4.1). Then, for each p > 1, there is a constant C p , independent of g, such that

∂v ∂ν C 1 (Γe) ≤ C p 1/p g L ∞ (Γi) . (4.6)
Based on (3.16), (3.17), and (4.6)

+∞ n=k ∂u n ∂ν C 1 (Γe) ≤ C p 1/p +∞ n=k u n-1 L ∞ (Γi) , (4.7)
for any k ≥ 1. From (3.14) we have

u n-1 L ∞ (Γi) = D Γe [ϕ n-1 ] + S Γi [ψ n-1 ] L ∞ (Γi) ≤ C p ϕ n-1 ψ n-1 Xp ,
which in combination with (3.5) and (3.10) yields

u n-1 L ∞ (Γi) ≤ C p C 1/p n-1 f L p (Γe) .
Insertion into (4.7) now gives

+∞ n=k ∂u n ∂ν C 1 (Γe) ≤ C p 1/p +∞ n=k C 1/p n-1 f L p (Γe) ≤ C k,p k/p 1 -C 1/p f L p (Γe) ,
for any k ≥ 1, provided C 1/p < 1. The preceeding analysis immediately leads to the following theorem.

Theorem 4.2 For all p > 1, there is a constant C p such that

∂u ∂ν - ∂u 0 ∂ν C 1 (Γe) ≤ C p 1/p f L p (Γe) . (4.8)
We now derive the desired formula for the principal term of ∂u ∂ν -∂u0 ∂ν . Using the notation from before we easily obtain that

∂u ∂ν (x) - ∂u 0 ∂ν (x) = - ∂u 1 ∂ν (x) + O( 2/p ) ,
uniformly on Γ e . We thus need to derive a formula for the principal term of ∂u1 ∂ν . For that purpose, let G(x, y) denote the Green's function for the problem (3.3), i.e., for each x ∈ Ω, G(x, y) is the solution to

           ∆ y G(x, y) = -δ x (y) in Ω , ∂ ∂ν y G(x, y) = 0, y ∈ Γ i , G(x, y) = 0, y ∈ Γ e . (4.9)
By continuity G(•, •) may be extended to all of Ω × Ω \ {(x, y) : x = y}, and it satisfies G(x, y) = G(y, x). In terms of G the solution v to (4.1) can be written

v(x) = Γi G(x, y)γ(y)g(y)dσ y , x ∈ Ω .
Remember: ν is the inward normal to Ω on Γ i . With g = u 0 this yields

u 1 (x) = Γi G(x, y)γ(y)u 0 (y)dσ y , x ∈ Ω ,
and, by taking the normal derivative on Γ e we now obtain

∂u 1 ∂ν (x) = Γi ∂ ∂ν x G(x, y)γ(y)u 0 (y)dσ y , x ∈ Γ e .
Note that ∂ ∂νx G(x, •) is a C 2 function on Γ i for any x ∈ Γ e and that u 0 is also C 2 on Γ i since Γ i is away from Γ e . Thus we have

∂ ∂ν x G(x, y)γ(y)u 0 (y) = m s=1 ∂ ∂ν x G(x, z s )u 0 (z s ) + O( ) γ s (y)χ s (y) ,
where χ s is the characteristic function of I s . We thus finally get

∂u 1 ∂ν (x) = m s=1 ∂ ∂ν x G(x, z s )u 0 (z s ) Is γ s dσ + O( 2 ) . Define γ s := Is γ s dσ .
Then we have arrived at the following formula for the principal term of

∂ ∂ν u -∂ ∂ν u 0 . Theorem 4.3 For x ∈ Γ e : ∂u ∂ν (x) - ∂u 0 ∂ν (x) = - m s=1 γ s u 0 (z s ) ∂ ∂ν x G(x, z s ) + O( 1+β ) (4.10)
for any 0 < β < 1.

MUSIC type algorithm for reconstruction

We now apply the asymptotic formula (4.10) to design a MUSIC (MUltiple SIgnal Classification) algorithm for locating small internal corrosive parts from boundary measurements. MUSIC is generally used in signal processing problems as a method for estimating the individual frequencies from multiple-harmonic signals [START_REF] Therrien | Discrete Random Signals and Statistical Signal Processing[END_REF]. In the present context the MUSIC algorithm uses information about the spectral structure of the (exterior Dirichletto-Neumann) boundary map for the Laplace operator with the "internal" corrosion (Robin) boundary condition. The eigenvectors corresponding to the most significant eigenvalues span a "signal subspace" in the sense that they contain nearly all the information about the corrosive parts that can be extracted from the boundary map. The remaining eigenvectors span a sort of "noise subspace". In order to be of any practical interest our algorithm uses spectral information about a particular discrete version of the boundary map.

It is worth mentioning that the present algorithm is also related to the linear sampling method of Colton and Kirsch [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF]. We refer to Cheney [START_REF] Cheney | The linear sampling method and the MUSIC algorithm[END_REF], and Kirsch [START_REF] Kirsch | The MUSIC algorithm and the factorisation method in inverse scattering theory for inhomogeneous media[END_REF], for detailed discussions of the connection between MUSIC-type algorithms and the linear sampling method.

Define the (exterior Dirichlet-to-Neumann) map

Λ γ from H 1/2 (Γ e ) into H -1/2 (Γ e ) by Λ γ (f ) = ∂u ∂ν Γe , (5.1)
where u is the solution to (1.4). Let Λ 0 be the exterior Dirichlet-to-Neumann map for the case when no corrosion is present. We proceed to show how partial spectral information about Λ γ -Λ 0 may effectively be used to determine the corroded internal boundary parts. The estimate (4.8) immediately gives that

(Λ γ -Λ 0 )(f ) C 1 (Γe) ≤ C p 1/p f L 2 (Γe) , (5.2)
for any 1 < p ≤ 2. In particular this shows that the operator Λ γ -Λ 0 , originally defined on H 1/2 (Γ e ), can be extended as an operator on all of L 2 (Γ e ). This operator is compact when regarded as an operator from L 2 (Γ e ) into itself. We also note that for g ∈ H 1/2 (Γ e )

Γe (Λ γ -Λ 0 )(f ) g dσ = Ω [∇u • ∇v -∇u 0 • ∇v 0 ] dx + Γi γu v dσ (5.3)
where v and v 0 are solutions to (1.4) and (1.5) with f replaced by g, respectively. Note that

Ω ∇u • ∇u 0 dx = Γe f ∂u 0 ∂ν dσ = Ω ∇u 0 • ∇u 0 dx.
Thus we obtain by substituting f in place of g in (5.3)

Γe (Λ γ -Λ 0 )f f dσ = Ω |∇(u -u 0 )| 2 dx + Γi γ|u | 2 dσ .
In summary we have established the following lemma.

Lemma 5.1 Λ γ -Λ 0 is self-adjoint, positive, and compact on L 2 (Γ e ).

The identity (4.10) shows that

(Λ γ -Λ 0 )(f )(x) = - m s=1 γ s u 0 (z s ) ∂ ∂ν x G(x, z s ) + O( 1+β ), x ∈ Γ e , (5.4)
for any 0 < β < 1. It is not difficult to see that the remainder term O( 1+β ) is actually bounded by C β 1+β f L 2 (Γe) with C β independent of f . Define the operator T on H 1/2 (Γ e ) by

(T f )(x) = - m s=1 γ s u 0 (z s ) ∂ ∂ν x G(x, z s ), x ∈ Γ e . (5.5)
Since u 0 depends linearly on f , T is a linear operator and can be extended as a bounded operator on L 2 (Γ e ). Furthermore the estimate we have just verified may easily be sharpened to show that Λ γ -

Λ 0 = T + O( 1+β ) , (5.6)
in the operator norm on L 2 (Γ e ).

Lemma 5.2

The operator T is a compact, self-adjoint, positive semi-definite operator on L 2 (Γ e ).

If γ s = c s for some fixed constants c s , then this follows immediately by multiplication of (5.4) by -1 , consideration of the → 0 limit, and use of Lemma 5.1. Here it is essential that the remainder term is bounded by C β 1+β f L 2 (Γe) , i.e., that we have norm convergence. More generally the lemma can also be verified by noting that T is a finite-dimensional operator, and that u 0 is given by

u 0 (y) = - Γe ∂ ∂ν x G(x, y)f (x) dσ(x), y ∈ Ω , (5.7) so that Γe T (f )h = m s=1 γ s u 0 (z s )v 0 (z s ) ,
where v 0 is the solution to (1.5), with f replaced by h. We introduce the linear operator G : L 2 (Γ e ) → R m Gf = (u 0 (z 1 ), . . . , u 0 (z m )) ,

where u 0 , as before, is the solution to (1.5). If R m is endowed with the standard Euclidean inner product then

Gf, a = m s=1 a s u 0 (z s ) = - Γe m s=1 a s ∂ ∂ν x G(x, z s ) f (x) dσ(x) ,
for arbitrary a = (a 1 , . . . , a m ) ∈ R m . Therefore, the adjoint operator G * : R m → L 2 (Γ e ) is given by

G * a = - m s=1 a s ∂ ∂ν x G(•, z s ) .
(5.9)

Exactly as in Brühl et al. [START_REF] Brühl | A direct impedance tomography algorithm for locating small inhomogeneities[END_REF] there is a simple factorization of the operator T , and a corresponding characterization of its range.

Lemma 5.3 (i) G * is injective;
(ii) G is surjective;

(iii) T = G * MG, where Ma = γ 1 a 1 , . . . , γ m a m , a = (a 1 , . . . , a m ) ∈ R m ;

(iv) Range(T ) = span ∂ ∂ν x G(x, z s ); s = 1, . . . , m .

Proof.

Suppose G * a = 0, and let w denote the function w(x) = m s=1 a s G(x, z s ). Then w = ∂w ∂ν = 0 on Γ e , and therefore by unique continuation w(x) = m s=1 a s G(x, z s ) = 0 in all of Ω. It immediately follows that a s = 0 for s = 1, . . . , m. Assertion (ii) follows from (i) and the well-known relation between the range and the null space of adjoint finite-dimensional operators: Range(G) = Ker(G * ) ⊥ . Using (5.5), (5.8), and (5.9), it is easy to see that (iii) holds. Now, according to (iii) we have Range(T ) = Range(G * MG) = Range(G * ), since M and G are surjective. This yields (iv), and the proof is complete.

The following result forms the basis for our algorithm to identify the midpoint locations, z s , of the small corroded boundary parts.

∂ ∂ν x G(•, z) Γe ∈ Range (T ) .
(5.10)

Proof. Assume ∂ ∂νx G(•, z) Γe ∈ Range (T ).
As a consequence of (iv) of Lemma 5.3 there exist coefficients a s , 1 ≤ s ≤ m such that

∂ ∂ν x G(x, z) = m s=1 a s ∂ ∂ν x G(x, z s ) , for x ∈ Γ e .
(5.11)

Since G(x, z) = 0 for any x ∈ Γ e and z ∈ Γ i , we also have

G(x, z) = m s=1 a s G(x, z s ) , for x ∈ Γ e , (5.12)
and by unique continuation it follows that G(x, z) = m s=1 a s G(x, z s ) for all x ∈ Ω. Due to the singularity of G(x, z) at x = z this is only possible if z ∈ {z s : s = 1, . . . , m}, and so we have established the sufficiency of the condition (5.10). The necessity follows immediately from (iv) of Lemma 5.3.

The finite-dimensional, self-adjoint operator T has a spectral decomposition

T = m p=1 λ p v p ⊗ v p , ||v p || L 2 (Γe) = 1 ,
where λ p , 1 ≤ p ≤ m, are the (positive) eigenvalues of T , and v p is the corresponding eigenfunction. Note that the rank of T , m, is exactly the number of corroded parts. We where the angle θ(z)

assume that λ 1 ≥ λ 2 ≥ . . . ≥ λ m > 0. Let P k : L 2 0 (Γ e ) → span {v 1 , . . . , v k }, k = 1, . . . ,
∈ [0, π/2) is defined by cot θ(z) = P m ∂ ∂νx G(•, z) Γe L 2 (Γe) (I -P m ) ∂ ∂νx G(•, z) Γe L 2 (Γe)
.

(5.14)

Since Λ γ -Λ 0 is a self-adjoint, positive, and compact operator on L 2 (Γ e ), it admits a spectral decomposition

Λ γ -Λ 0 = +∞ p=1 λ p v p ⊗ v p , ||v p || L 2 (Γe) = 1 , (5.15) with λ 1 ≥ λ 2 ≥ . . . ≥ λ m ≥ . . . > 0. Let P k : L 2 (Γ e ) → span {v 1 , . . . , v k }, k = 1, 2, . . . , be the orthogonal projector P k = k p=1 v p ⊗ v p .
The data for our identification algorithm consist of (an appropriate approximation to) the spectral decomposition of Λ γ -Λ 0 . Motivated by (5.10), we seek to find those points

z ∈ Γ i such that ∂ ∂ν x G(•, z) Γe ∈ Range (Λ γ -Λ 0 ) . (5.16)
This condition is fulfilled approximately by z ∈ Γ i if the following quantity is sufficiently large:

cot θ k (z) := P k ∂ ∂νx G(•, z) Γe L 2 (Γe) (I -P k ) ∂ ∂νx G(•, z) Γe L 2 (Γe)
.

(5.17)

Since we don't a priory know m, we plot the angles θ k (z) (for k = 1, 2, . . .) as a function of z. When the plots "stabilize", i.e., don't significantly change by increasing k, we have a good estimate for m, and we interpret locations where small angles are attained as good candidates for the z s .

Once the locations z s are approximately found, we can also estimate the integrated corrosion coefficients γ s . Our procedure for doing this is the following. Define w s , s = 1, . . . , m, to be the solution to

             ∆w s = 0 in Ω , - ∂w s ∂ν = 0 on Γ i , w s (x) = ∂ ∂ν x G(x, z s ) x ∈ Γ e . (5.18)
It then follows from (5.4) and (5.15) that

- m s=1 γ s w s (z s ) ∂ ∂ν x G(•, z s ) ≈ (Λ γ -Λ 0 ) ∂ ∂ν x G(•, z s ) ≈ m p=1 λ p v p , ∂ ∂ν x G(•, z s ) v p .
By integrating both sides of the above formula against v s , we obtain

- m s=1 w s (z s ) v s , ∂ ∂ν x G(•, z s ) γ s ≈ λ s v s , ∂ ∂ν x G(•, z s ) , (5.19)
for s = 1, . . . , m. Therefore the values of γ s , s = 1, . . . , m, can be calculated by solving a linear system. For the simplest case, m = 1, the formula reads

γ 1 ≈ - λ 1 w 1 (z 1 )
.

Numerical Results

In this section we present some numerical experiments with the MUSIC-type algorithm we have designed to find the internally corroded parts,

I s ⊂ R 2 , s = 1, 2, • • • , m. The domain Ω ⊂ R 2 is
taken to be an annulus, i.e., Ω = U \ D, where U and D are disks centered at (0, 0) with radii r e and r i , respectively. The measurements used for the identification represent a discrete analogue of the Dirichlet-to-Neumann map Λ γ -Λ 0 on Γ e . In order to find a good approximation from the "identification point of view" the following observation will be useful.

Lemma 5.5 The functions ∂G ∂νx (•, y), y ∈ Γ i span a dense subspace of L 2 (Γ e ).

Proof. The solution to the boundary value problem (1.5), for a given f ∈ L 2 (Γ e ), is given by

u 0 (y) = - Γe ∂ ∂ν x G(x, y)f (x) dσ(x) .
If f ∈ L 2 (Γ e ) were orthogonal to { ∂G ∂νx (•, y) : y ∈ Γ i } it thus follows that ∂u 0 ∂ν = u 0 = 0 on Γ i .

By unique continuation this implies that u 0 = 0 in all of Ω, and so consequently f = 0.

With this observation in mind we approximate Λ γ -Λ 0 by using Dirichlet data of the form f (•) = ∂G ∂νx (•, y), y ∈ Γ i , and solving the integral equations (3.1) and (3.7). Note that G is a quite particular fundamental solution, namely the solution to (4.9). Practically it is very important that these normal derivatives can be calculated efficiently. To understand how, let G be the Green's function for U , that is, for each x ∈ U , let G(x, y) be the solution to

   ∆ y G(x, y) = -δ x (y) in U, G(x, y) = 0 y ∈ Γ e .
(5.20)

Since U is a disk, G(x, y) has a simple explicit expression. In fact, G is given by

G(x, y) = - 1 2π ln |x -y| + 1 2π ln r e |x| x - |x| r e y . (5.21)
For any y ∈ Γ i we may now compute ∂G ∂νx (•, y) on Γ e by solving the integral equation (5.22) of the following lemma. Lemma 5.6 For any y ∈ Γ i , let G y (x) := G(x, y), and G y (x) := G(x, y) for x ∈ Γ e . Then

I + 2K * Γe -4 ∂S Γi ∂ν e ∂S Γe ∂ν i ∂G y ∂ν e (x) = 2 ∂ G y ∂ν e (x), x ∈ Γ e , (5.22)
where ν e and ν i denote outward normal vectors to U on Γ e , and to D on Γ i , respectively.

Proof. Fix y ∈ Ω and let u(x) = (G y -G y )(x) for x ∈ Ω. Because of the symmetry of G and G, u satisfies ∆u = 0 in Ω, u = 0 on Γ e , and ∂u ∂νi = -∂ e Gy ∂νi on Γ i . Thus by the Green theorem, we have

u(x) = -S Γe ∂u ∂ν e (x) + D Γi [u](x) + S Γi ∂u ∂ν i (x), x ∈ Ω .
However, in view of (2.3), we have

D Γi [u](x) = S Γi 1 2 I + K * Γi -1 ∂ ∂ν i D Γi [u] Γi (x)
for all x ∈ R 2 \ D. Therefore there is a density function ψ with Γi ψ = 0 such that

(G y -G y )(x) = -S Γe ∂G y ∂ν e - ∂ G y ∂ν e (x) + S Γi [ψ](x) , x ∈ Ω . (5.23)
Similarly to (3.1) we now obtain from (5.23)

1 2 I + K * Γi [ψ] = - ∂ G y ∂ν i + ∂S Γe ∂ν i ∂G y ∂ν e - ∂ G y ∂ν e on Γ i . Since Γ i is a circle, K * Γi [ϕ] = 0 when Γi ϕ = 0 (cf. (2.6)) and hence ψ = -2 ∂ G y ∂ν i + 2 ∂S Γe ∂ν i ∂G y ∂ν e - ∂ G y ∂ν e . Since G y (x) + S Γe [ ∂ e Gy ∂νe ](x) = -Φ(x,
(x) = - 1 2 Φ x (y) + 1 4πr i Γi Φ x (y)dσ(y) = - 1 4π [ln |x -y| -ln |x|] ,
which is the desired formula. Note that the last equality holds because of (2.7).

A similar statement holds for any derivative with respect to x, |x| > r i . In view of (5.21), it is now a matter of straight-forward computation to check that in the limit as y approaches

Γ i from inside Ω lim K * Γe ∂ G y ∂ν e + 2 ∂S Γi ∂ν e ∂Φ y ∂ν i = 1 2 ∂ G y ∂ν e on Γ e ,
and we obtain (5.25). This completes the proof.

For computation, we discretize Γ e and Γ i by choosing points The discrete data we use for our numerical algorithm consist of the 256 × 256 matrix ((Λ γ -Λ 0 )( ∂G ∂νx (•, y n ))(x m )) 256 m,n=1 , or rather its SVD (singular value decomposition). For comparison we also calculate the 256 × 256 matrix (T ( ∂G ∂ν (•, y n ))(x m )) 256 m,n=1 and its SVD. Notice than none of these two matrices are symmetric.

Figure 1 displays the singular values of our discrete versions of Λ γ -Λ 0 and T when there are two internally corroded parts. The singular values of our approximation to T exhibit a significant drop, to machine epsilon, after the first two (consistent with the fact that the range of T is two-dimensional) there is also a drop in the singular values of our approximation to Λ γ -Λ 0 , but it is of course not nearly as significant due to the discrepancy between Λ γ -Λ 0 and T , and due to the very approximate nature of our discrete data. It is nonetheless still possible to estimate the location of the corroded parts very accurately from the SVD of our discrete version of Λ γ -Λ 0 .

Using the first k (left-)singular vectors of our data, we plot approximate values of θ k (z) = cot -1 L k (z) where

L k (z) := P k ( ∂G(•,z) ∂ν ) L 2 (Γe) (I -P k )( ∂G(•,z) ∂ν ) L 2 (Γe)
.

(5.27)

We start this process with just one singular vector (k = 1) and continue with increasing k until the plot of θ k (z) stabilizes. Values near zero correspond to locations are that likely to be close to the I s (or z s ).

Example 1. In this example, the outer radius r e = 1 and the inner one r i = 0.8 and there are two corroded parts. The top part of Figure 1 shows the locations of the corrosion and the approximate SVD of Λ γ -Λ 0 and T . The bottom of Figure 1 clearly shows that there are two corroded parts, since the minima of θ k stabilize already at k = 2. The minima of θ 2 are the locations of the corroded parts. Since we use the same 256 points for the calculation and plotting of θ k and for the potential location of the points z s , the fact that |z s -z c s | = 0 really asserts that z s is accurately estimated up to the spacing between two adjacent mesh points, which is approximately 0.0172. By solving (5.19), we compute the integrated corrosion coefficients: γ c 1 = 0.3846, γ c 2 = 0.1660. The actual, integrated coefficients are γ 1 = 0.3927, γ 2 = 0.1767.

Example 2. In this example we consider a case with five corroded parts. The actual data of the configuration is summarized in the top of Table 1. The corroded parts are numbered counter-clock-wise starting with the smallest positive angle. Note that the first three corroded parts have low corrosion coefficients while the remaining two have relatively high ones. The computational results with 1% noise are summarized in Table 1 and Figure 2. It is interesting to note that the corroded parts 1 and 2 which are close to each other and have low corrosion coefficients, are detected as a single one. The actual midpoint locations of corroded parts 1 and 2 are (0.0686, 0.6966) and (-0.1197, 0.6897), and the single detected location is (-0.1027, 0.6924). On the other hand, the corroded parts [START_REF] Brühl | A direct impedance tomography algorithm for locating small inhomogeneities[END_REF] The results for steps number of sing. vectors(k) In the table, m c is the number of estimated corroded parts, z s and z c s are the actual midpoint locations and the estimated midpoint locations, and γ s , γ c s are the actual and the estimated, integrated, corrosion coefficients, for a particular k. Example 3. Figure 3 shows the computational results with various degree of noise. The information about the actual locations and corrosion coefficients is summarized in Table 2. Observe that the first two corroded parts have low corrosion coefficients. They are detected as a single corroded part with 5% noise. The other two corroded parts, which have high corrosion coefficients, are detected very well even with a high ratio of noise.

m c z c s |z s -z c s | γ c s | γ s -γ c s | 1

Conclusion

We have designed a non-iterative algorithm of MUSIC-type for detecting small, internal, corroded boundary parts based on external, electrostatic boundary measurements. Our The numerical tests clearly demonstrate that the algorithm works well even in the presence of relatively high noise ratios, provided the corrosion coefficients are sufficiently large. The latter restriction is to be expected, since only those corroded parts with reasonably sized corrosion coefficients can cause sufficient changes in currents. Furthermore, it is practically most important to detect the corroded parts with the highest coefficients, since they are likely to "develop" more rapidly, and therefore potentially cause the most serious damage (to the pipe). It should be mentioned that our algorithm cannot identify separately the size of the corroded parts (without assuming that we know the corrosion coefficients). We can only reconstruct the integral γ s = Is γ dσ for s = 1, . . . , m. It is likely that a higher-order asymptotic expansion of the boundary perturbations would yield formulas that permit such identification.

Theorem 5 . 4 A

 54 point z ∈ Γ i belongs to the set {z s : s = 1, . . . , m} if and only if

  m, be the orthogonal projector P k = k p=1 v p ⊗ v p . It then follows from Theorem 5.4 that z ∈ {z s : s = 1, . . . , m} iff (I -P m ) ∂ ∂ν x G(•, z) Γe = 0 , or equivalently, z ∈ {z s : s = 1, . . . , m} iff cot θ(z) = +∞, (5.13)

  {r e (cos θ n , sin θ n )|θ n = 2π(n -1)/N, n = 0, 1, . . . , N -1} and {r i (cos θ n , sin θ n )|θ n = 2π(n -1)/N, n = 0, 1, . . . , N -1}, with N = 256. Put x n := r e (cos θ n , sin θ n ) and y n := r i (cos θ n , sin θ n ) for n = 0, 1, . . . , 255.

Figure 1 :

 1 Figure 1: Pipe with two internal corroded parts, and "approximate" singular values of Λ γ -Λ 0 and T in upper figures. The bottom figure shows θ k for three different values of k.In the table, m c is the number of estimated corroded parts, z s and z c s are the actual midpoint locations and the estimated midpoint locations, and γ s , γ c s are the actual and the estimated, integrated, corrosion coefficients, for a particular k.

Figure 2 :Table 2 :

 22 Figure 2: The computational results with 1% noise. The last plot shows the values of (Λ γ -Λ 0 )( ∂G ∂ν ) on Γ e × Γ i .

Figure 3 :

 3 Figure 3: Computational results with varying degree of noise

Table 1 :

 1 and 5, which have high corrosion coefficients, are clearly detected already with k = 3. Summary of computational results with five corroded parts, and 1% noise. Two corroded parts with low corrosion coefficient are detected as a single one.

		m γ s	actual data z s		γ s
		5 0.1	(0.0686, 0.6966)	0.0120
			0.5 (-0.1197, 0.6897) 0.0430
			0.3 (-0.4307, 0.5518) 0.0206
			1.0 (0.2519, -0.6531) 0.1203
			1.2 (0.5723, -0.4031) 0.1031
	m c		z c s	detected data γ c s |z s -z c s | | γ s -γ c s |
	4	(-0.1027, 0.6924) 0.0514	0.0172	0.0084
		(-0.4170, 0.5622) 0.0212	0.0172	0.0006
		(0.2519, -0.6531) 0.1181	0		0.0021
		(0.5723, -0.4031) 0.1026	0		0.0005
				result for each step
	number of sing. vectors(k)	m c		z c s	|z s -z c s |	γ c s	| γ s -γ c s |
	1	1	(0.3889, -0.5820)			0.1554
	2	2	(-0.1366, + 0.6865)			0.0579
			(0.3889, -0.5820)			0.1554
	3	3	(-0.1366, + 0.6865)			0.0579
			(0.2519, -0.6531)			0.1181
			(0.5723, -0.4031)			0.1026
	4		(-0.0686, + 0.6966)	0.0515	0.0486	0.0056
			(-0.4031, + 0.5723)	0.0344	0.0237	0.0030
			(0.2519, -0.6531)		0	0.1181	0.0021
			(0.5723, -0.4031)		0	0.1026	0.0005
	5	4	(-0.1027, + 0.6924)	0.0172	0.0514	0.0084
			(-0.4170, + 0.5622)	0.0172	0.0212	0.0006
			(0.2519, -0.6531)		0	0.1181	0.0021
			(0.5723, -0.4031)		0	0.1026	0.0005
	6	4	( -0.1027, + 0.6924)	0.0172	0.0514	0.0084
			(-0.4170, + 0.5622)	0.0172	0.0212	0.0006
			(0.2519,-0.6531)		0	0.1181	0.0021
			(0.5723, -0.4031)		0	0.1026	0.0005

Acknowledgments. We would like to thank the referees for their careful reading of the manuscript and their many suggestions. We feel the implementation of some of these suggestions has made our presentation more accessible.

* H.A. is partially supported by the Brain Pool Korea Program at Seoul National University, H.K. is partially supported by KOSEF grant R01-2006-000-10002-0, E.K. is supported by BK21 Math. Division at Seoul National University, and M.S.V. is partially supported by NSF grant DMS-0604999.


David Sferruzza
email: david.sferruzza@ls2n.fr

Top-Down Model-Driven Engineering of Web Services from Extended OpenAPI Models

Keywords: Domain specific languages, Software prototyping, Web Services, Model-Driven Engineering, Documentation, Code Generation, OpenAPI 3.0

Web services engineering is a crucial subject, because web services are often built to be used by other programs; thus they should have a good documentation targeting developers. Furthermore, when building a digital product, engineers need to build several programs that interact with a central instance of web services. OpenAPI, a popular industry standard, makes possible to document web services in order to quickly make a prototype of the product. It allows a top-down process where developers iterate to build an OpenAPI model that describes the web services they want, and then implement both the web services and the programs that will consume them. However, when doing such rapid prototyping, developers tend to either skip this design phase and implement web services right away, or stop updating the OpenAPI model when the product is released; in both cases they cannot take advantage of having an OpenAPI model aligned with the implementation. We show how OpenAPI can be extended to add implementation details inside models. These extensions link services to assemblies of components that describe computations. Hence a top-down development process that keeps model and implementation aligned. Moreover, this makes possible for developers to benefit from more support features while keeping the same flexibility.

CCS CONCEPTS

• Software and its engineering → Model-driven software engineering;

INTRODUCTION

Context. Web services are key elements of many digital products. Indeed, these products often require to be available on several platforms (web, Android, iOS, . . . ) thus can benefit from separating concerns: web services centralizing data and business logic have different life cycles than the several user interfaces that rely on them. This architecture is quite common; for example, after several years of coaching startups at Startup Palace 1 , it appears that a lot of them are using it or should be. One of the difficulties of this approach is to design and maintain a consistent interface between the web services and the other programs, which are sometimes developed by other teams or even unrelated people. Fortunately, many tools and languages exist to help designing, writing and maintaining such interfaces. OpenAPI [START_REF]OpenAPI Specification[END_REF] is one of them: it defines a specification to describe the HTTP APIs of web services in a language-agnostic way. It is quite famous in the industry and has a rich ecosystem of tools 2 .

Motivation. Building or fast-prototyping digital products can be achieved using a top-down approach. First, they design an OpenAPI model and make sure it is suitable as an interface between the different actors. Then, they can rely on this model to implement them independantly using programming languages. Some tools like Swagger Code Generator [START_REF]Swagger Code Generator[END_REF] provide support to developers in this process: they generate code based on the OpenAPI model, targeting one of the numerous supported languages and frameworks. However, because OpenAPI is language-agnostic, the generated code does not contain implementation details; developers still have to complete it in a way that respects the contract defined by the OpenAPI model, so that they obtain a working implementation. Focusing on web services, this means that, apart from this one-time generation, the OpenAPI model and the actual implementation will need to be maintained in a separate manner, which in many cases results in the implementation being maintained and the OpenAPI model being outdated. This becomes more problematic when the digital product evolves a lot (i.e. it is a MVP3 ) because developers cannot benefit again from the code generation step when doing another top-down cycle as it would override precedently customized code.

Contributions. To overcome this issue and give even more support to developers, we extend the OpenAPI 3.0 Specification to make possible to add some implementation details in models, and provide a tool to generate working web services from extended models.

Implementation details added to OpenAPI models consist of assemblies of components that represent computation units. These components can be defined in two ways: (i) by a contract and an implementation in a programming language, or (ii) in term of already defined components. Every service is linked to a component whose purpose is to generate an HTTP response.

Our tool, SWSG [START_REF] Sferruzza | Safe Web Services Generator[END_REF], provides support to automatizes the topdown development process. It can read these extended models, check several properties to ensure their consistency, and generate code of working web services. On top of solving the alignment problem, this gives even more support to developers.

The article is structured as follows. Section 2 describes the component system and how it is integrated in OpenAPI 3.0. Section 3 introduces SWSG and shows how it checks models and generates code. Section 4 presents related work. Finally, Section 5 concludes the article with some lessons and future work.

EXTENDING OPENAPI

The OpenAPI Specification defines a standard to express interfaces to HTTP APIs in a language-agnostic way. It aims at allowing "both humans and computers to discover and understand the capabilities of the service without access to source code, documentation, or through network traffic inspection" [START_REF]OpenAPI Specification[END_REF]. As in Model-Driven Engineering (MDE), the point of having such meta-models is to have tools that can rely on them in order to safely manipulate models and offer support to developers. Indeed, the growing ecosystem of OpenAPI contains numerous tools that can be used in a top-down process, which often means that an OpenAPI is given as input to a tool that will refine it and produce another artifact. For exemple, both Swagger Editor [START_REF]Swagger Editor[END_REF] and Swagger UI [START_REF]Swagger UI[END_REF] provide an interactive graphical user interface from a model, and Dredd [1] generates functional tests from a model.

By design, OpenAPI does not contain implementation details that describe how the web services are to be implemented. As stated in the introduction, this can lead to a misalignment between the OpenAPI model and its implementation. To overcome this limitation while keeping the benefits of using OpenAPI models in a top-down process, we extend OpenAPI to add a new layer of information that describes the implementation of the web services at a high-level of abstraction. This layer consists of components definitions and instanciations such as those we introduced in a previous work [START_REF] Sferruzza | A Model-Driven Method for Fast Building Consistent Web Services in Practice[END_REF].

Section 2.1 describes this new layer, and Section 2.2 shows how it is integrated in OpenAPI 3.0.

The Component System

Components are units of processes and computations that occur inside web services and whose purpose is to produce HTTP responses. Their execution happens in an isolated context, that can contain variables. They can mutate this context by adding and removing variables, or return an HTTP response. Components can be of two kinds: either atomic or composite.

Atomic components are defined by a name and four sets of variables: parameters, preconditions, additions and removals. Parameters are variables whose values must be provided when instanciating the component; it allows to design generic components that are easier to reuse. The other three sets of variables form the contract of the component: they define what variables it needs to access from the context (preconditions), what variables it will add to the context (additions) and what variables it will remove (removals). The point of this contract is to support a lightweight kind of static verification we call structural consistency. Atomic components must go along with an implementation in a programming language, which makes them very flexible.

Composite components are defined by a name, a set of parameters and a list of instances of components. Parameters work in the same way as for atomic components. Instances of components reference a component by its name and provide bindings and aliases, that are specific to the instance. Bindings associate a parameter of the component to a value, and aliases allow to rename a variable from the component's contract in a local way. Composite components do not need to be implemented using a programming language because they are defined in term of other components. Indeed, their behavior consists of executing their subcomponents sequentially, each getting the context output by the previous. However, this sequential execution is interrupted if any of the components outputs an HTTP response instead of a new context.

This component system was voluntarily designed to be simple and offer a good trade-off between expressiveness and ease of support. That is why, for example, sequential execution is the only way to combine components.

Extensions to OpenAPI 3.0

To make it possible to use this component system and describe implementations of web services from inside an OpenAPI model, the OpenAPI 3.0 Specification must be extended. Two kinds of informations must be added: definitions of components and, for each service, an instance of component. Hopefully, OpenAPI provides an extension mechanism that preserves tools compatibility: most of the schemas of the specification can be enhanced with more attributes, if their names start with x-. The full specification of our extensions is available in [START_REF] Sferruzza | Specification of SWSG extensions for OpenAPI[END_REF].

The Listing 1 shows an extract of the Petstore example [START_REF]The Petstore Example[END_REF], an official example of OpenAPI 3.0 model, that has been enhanced with our extensions. The components object contains two new properties: x-swsg-cc and x-swsg-ac. The first is a set of composite component definition; here we define a component called FindPet in terms of two other components named GetPetById and RenderPet. The second is a set of atomic component definition that contains the definitions of GetPetById and RenderPet. The Listing 2 shows how the FindPet component is instanciated from a service, in the x-swsg-ci attribute.

GENERATING WEB SERVICES

To use our extensions to OpenAPI 3.0, we propose the following top-down development process:

(1) developers write a standard OpenAPI model that fits their needs; The step (2) is made possible by our extensions and provide a transitional step between the modelling (1) and the implementation (3) in the refinement process. It simplifies the implementation (3) by reducing it to several small functions instead of a whole program.

The step (4) makes use of our tool, SWSG [START_REF] Sferruzza | Safe Web Services Generator[END_REF], that follows the process shown in Figure 1. We present this tool in the following sections: Section 3.1 focuses on the model consistency verification and Section 3.2 on the code generation.

Consistency Verification

The third step of the process of SWSG (Figure 1) is a static verification of the structural consistency of the model. It consists of a set of formal rules that check the semantics of models to determine whether it allows for a safe code generation. This set of rules is derived from the one we presented in [START_REF] Sferruzza | A Model-Driven Method for Fast Building Consistent Web Services in Practice[END_REF][START_REF] Bernardi | Automated development of constraint-driven web applications[END_REF]. Its purpose is to help developers to spot inconsistencies in models as earlier as possible.

For example, the model presented in Listings 1 and 2 is inconsistent: SWSG indicates a PreconditionError. Details of the error tell us that the GetPetById component is not given a variable of type String called id, as it requires. Indeed, the service definition shows that the id variable comes from a service parameter that has the Integer type. This kind of mistakes is quite common in software development, but having the right tools can mitigate its consequences. Here, SWSG detected the problem before we even run the code in production, so we just lost a small amount of time; we could have lost a lot more time if this code had been deployed (just for a typo).

When this error is fixed in the model, SWSG can proceed to the code generation step.

Code Generation

The process depicted by Figure 1 is generic: it does not rely on a specific language or technology. Yet the language and technologies used to implement atomic components must be identical or compatible with those of the code generation target. Because we experiment in Startup Palace's context, the current implementation of SWSG targets the PHP programming language 4 with the Laravel web framework 5 , which is a common tool stack.

The generation in itself is done with a model-to-text transformation that uses a template engine called Twirl6 . Apart from some static classes and from the implementations of atomic components (that are not altered), SWSG generated two kinds of source files: implementations of composite components and a route file. The Listing 3 shows the code generated for the FindPet component (defined in Listing 1).

We took the approach of generating code that should never be manually edited. The generated code does not override any existing Listing 3: Generated Code for the FindPet Component files in a Laravel's architecture and can be easily hooked to an existing web application through configuration.

In case of evolution of the needs (which is likely when working with MVPs), developers have to go through the top-down development process again, but can reuse most of the existing items (parts of the extended OpenAPI model, implementation of atomic components). This approach preserves alignment between the model and the web services, and fosters reuse of components.

RELATED WORK

The use of MDE for development and automatic generation of web services or web applications is not a new topic [START_REF] Bernardi | M3D: a tool for the Model Driven Development of Web Applications[END_REF][START_REF] Bernardi | Automated development of constraint-driven web applications[END_REF][START_REF] Scheidgen | Metamodeling vs Metaprogramming: A Case Study on Developing Client Libraries for REST APIs[END_REF]. Indeed, this work is built on top of the approach of SWSG [START_REF] Sferruzza | A Model-Driven Method for Fast Building Consistent Web Services in Practice[END_REF] and Reifier [START_REF] Rocheteau | Reifier: Model-Driven Engineering of Component-Based and Service-Oriented JEE Applications[END_REF].

Through SWSG, our process shares the meta-modelling approach with tools such as M3D (introduced in [START_REF] Bernardi | M3D: a tool for the Model Driven Development of Web Applications[END_REF] and extended in [START_REF] Bernardi | Automated development of constraint-driven web applications[END_REF]) that also focus on building web services using MDE. One of the main differences between SWSG and M3D is that SWSG was developed with a focus on design-time support. Even if SWSG is definitely related to existing standards such as BPEL [START_REF] Fu | Analysis of interacting BPEL web services[END_REF] or WSDL, our approach differs on several aspects. First, we want to avoid the shortcomings described in [START_REF] Gronmo | Model-driven web services development[END_REF], that is WSDL models contain too much technical details and are difficult to understand for humans. Indeed our metamodel is simpler and less expressive than WSDL or BPEL. Second, this allows SWSG to provide more support to users; the balance between flexibility and support is discussed in [START_REF] Wil | Declarative workflows: Balancing between flexibility and support[END_REF].

One of the tools featured in the OpenAPI ecosystem is Swagger Code Generator [START_REF]Swagger Code Generator[END_REF]. It aims at generating client librairies, server stubs or documentations from an OpenAPI model. It supports many languages and frameworks, but only helps developers to write new services by generated boilerplate code. This automatizes a tedious task; however they still need to add a lot of code on top of it. Moreover, when services evolve, developer need to manually propagate evolutions into the codebase because Swagger Code Generator isn't able to merge them automatically. Our approach solves this issue because it gives flexibility to developers before the code generation step, making useless editing generated code. This ensures alignment between the (OpenAPI) model and the web services.

CONCLUSION

We extended OpenAPI 3.0 with a component system to describe implementations of web services from a high-level. This establishes an intermediary artifact that fits well in the top-down development process, between the abstract OpenAPI model and the concrete implementation of web services. This allows to generate an implementation of the web services that does not need manual modifications. As a consequence, alignment between the OpenAPI model and the implementation is guaranteed. Moreover, this paves the way for enhanced support features; for example, verification or interactive visualization. We also presented SWSG, a tool that can check structural consistency of extended OpenAPI models and generate code of the web services.

This approach was tested on several case studies 7 , but needs more testing on bigger projects. Another perspective is to generalize SWSG to support more languages and technologies. This would probably lead to the formalization of interfaces to make it possible to plug external code generators. Also, it would be interesting to allow SWSG to check the compliance between implementations of atomic components and their contract. In a similar manner, calling external tools could be a solution to keep SWSG generic.

Figure 1 :

 1 Figure 1: Process of SWSG

https://www.startup-palace.com

https://github.com/OAI/OpenAPI-Specification/blob/master/IMPLEMENTATIONS. md

Minimum Viable Product: a finished but lightweight product made to test market hypotheses.

https://php.net/

https://laravel.com/

https://github.com/playframework/twirl