A -C Bennis

T Chac

Ó N Rebollo

M Gomez Marmol

R Lewandowski

Stability of some turbulent vertical models for the ocean mixing boundary layer

Keywords: oceanography, turbulence models, stability, partial differential equations MSC classification : 35J60, 35K55, 76E20, 76F40

We consider four turbulent models to simulate the boundary mixing layer of the ocean. We show the existence of solutions to these models in the steady-state case then we study the mathematical stability of these solutions.

Introduction

The presence of an homogeneous layer near the surface of the ocean has been observed since a long time. The so called "mixed layer" presents almost constant profiles of temperature and salinity (or equivalently the density). The bottom of the mixed layer corresponds either to the top of the thermocline, zone of large gradients of temperature, or to the top of the zone where haline stratification is observed [START_REF] Vialard | An ogcm study for the toga decade. part i: Role of salinity in the physics of the western pacific fresh pool[END_REF]. Some attempts to describe this phenomenon can be found for example in Defant [START_REF] Defant | Schichtung und zirculation des atlantischen ozeans[END_REF] or Lewandowski [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF]. The effect of the wind-stress acting on the sea-surface was then considered to be the main forcing of this boundary layer. Observations in situ were completed by laboratory experiments [START_REF] Deardorff | Laboratory investigation of nonsteady penetrative convection[END_REF] and more recently by numerical modelizations of the mixed layer.

In this note, we consider four turbulent models to describe this homogeneous boundary layer. The first one is the Pacanowski-Philander model, and two of these models are new models. They aim to compute the velocity and the water density of a water column, are one space dimensional and the eddy viscosities depend on the Richardson number. For those model, we show the existence of a steady-state solution and we analyse the mathematical linear stability of these steady state solution, showing that only one of these model, the one we introduce in this note (model labelized as R -2 -2 -4 below), has a unique staedy state solution with a large range of stabilty. Moreover, in [START_REF] Bennis | Parametrization of the mixing layer: Comparison of four models depending on the richardson number[END_REF] we have used these models to simulate the warm pool at the equator. Numerical results confirm that R -2 -2 -4 is the most accurate parametrization.

The equations

We denote by (u, v) the horizontal water velocity and ρ its density. Since the numerical simulation performed in [START_REF] Bennis | Parametrization of the mixing layer: Comparison of four models depending on the richardson number[END_REF] concerns the equator zone, we do not take the Coriolis force into account. The closure equations are:

(2.1)

                                       ∂u ∂t - ∂ ∂z ν 1 ∂u ∂z = 0, ∂v ∂t - ∂ ∂z ν 1 ∂v ∂z = 0, ∂ρ ∂t - ∂ ∂z ν 2 ∂ρ ∂z = 0, for t 0 and -h z 0, u = u b , v = v b , ρ = ρ b at the depth z = -h, ν 1 ∂u ∂z = ρ a ρ 0 V x , ν 1 ∂v ∂z = ρ a ρ 0 V y , ν 2 ∂ρ ∂z = Q at the surface z = 0, u = u 0 , v = v 0 , ρ = ρ 0 at initial time t = 0.
In system (2.1), the coefficients ν 1 and ν 2 are the vertical eddy viscosity and diffusivity coefficients and will be expressed as functions of the Richardson number R defined as

R = - g ρ 0 • ∂ρ ∂z ∂u ∂z 2 +
∂v ∂z

2
where g is the gravitational acceleration and ρ 0 a reference density (ρ 0 ≃ 1025 kg.m -3 ).

The constant h denotes the thickness of the studied layer that must contain the mixing layer. Therefore the circulation for z < -h, under the boundary layer, is supposed to be known, either by observations or by a deep circulation numerical model. This justifies the choice of Dirichlet boundary conditions at z = -h, u b , v b and ρ b being the values of horizontal velocity and density in the layer located below the mixed layer. The air-sea interactions are represented by the fluxes at the sea-surface : V x and V y are respectively the forcing exerced by the zonal wind-stress and the meridional wind-stress and Q represents the thermodynamical fluxes, heating or cooling, precipitations or evaporation. We have

V x = C D |u a | 2 and V y = C D |v a | 2
, where U a = (u a , v a ) is the air velocity and C D a friction coefficient.

We study hereafter four different formulations for the eddy coefficients

ν i = f i (R), labeled as "R -2 -i" and/or "R -2 -i -j" . In all models, f 1 (R) = α 1 + β 1 (1 + 5R) 2 , except in model R -2 -3 below: (2.2) R -2 -1 -3 : f 2 (R) = α 2 + f 1 (R) 1 + 5R = α 2 + α 1 1 + 5R + β 1 (1 + 5R) 3 . R -2 -3 f 1 (R) = α 1 + β 1 (1 + 10R) 2 , f 2 (R) = α 2 + β 2 (1 + 10R) 3 R -2 -2 -4 : f 2 (R) = α 2 + f 1 (R) (1 + 5R) 2 = α 2 + α 1 (1 + 5R) 2 + β 1 (1 + 5R) 4 , R -2 -2 f 2 (R) = α 2 + β 2 (1 + 5R) 2 ,
Formulation R -2 -1 -3 corresponds to the modelization of the vertical mixing proposed by Pacanowski and Philander [START_REF] Pacanowski | Parametrization of vertical mixing in numericals models of the tropical oceans[END_REF]. The coefficients α 1 , β 1 and α 2 have the following values: α 1 = 1.10 -4 , β 1 = 1.10 -2 , α 2 = 1.10 -5 (units:m 2 s -1 ).This formulation has been used in the OPA code developed in Paris 6 University [START_REF] Madec | version 8.0. ocean general circulation model[END_REF] with coefficients α 1 = 1.10 -6 , β 1 = 1.10 -2 , α 2 = 1.10 -7 units: m 2 s -1 .The selection criterion for the coefficients appearing in these formulas was the best agreement of numerical results with observations carried out in different tropical areas. Formulation R -2 -3 has been proposed by Gent [START_REF] Gent | The heat budget of the toga-coare domain in an ocean model[END_REF]. Formulations R -2 -2 -4 and R -2 -2 are new as far as we know. Notice that models R -2 -1 -3 and R -2 -3 are no more physically valid respectively for R ∈ (-3.13, -0.2) and R ∈ (-2.25, -0.1) since the coefficient ν 2 becomes negative.

Steady-state solutions

Steady-state solutions to system (2.1) satisfy

(2.3) ∂ ∂z f 1 (R) ∂u ∂z = 0, ∂ ∂z f 1 (R) ∂v ∂z = 0, ∂ ∂z f 2 (R) ∂ρ ∂z = 0. Theorem 2.1 System (2.
3) has at leat one smooth solution on [0, -h] for each model in (2.2). In case of R -2 -2 -4 the solution is unique.

Proof. Integrating (2.3) with respect to z, yields

(2.4) f 1 (R) ∂u ∂z = V x ρ a ρ 0 , f 1 (R) ∂v ∂z = V y ρ a ρ 0 , f 2 (R) ∂ρ ∂z = Q. and since R = - g ρ 0 • ∂ρ ∂z ∂u ∂z 2 +
∂v ∂z

2 we deduce from (2.4) that R = - gQρ 0 ρ 2 a (V 2 x + V 2 y ) • (f 1 (R)) 2 f 2 (R) , which yields (2.5) (f 1 (R)) 2 f 2 (R) = - ρ 2 a (V 2 x + V 2 y ) gQρ 0 R
which is a fixed point equation for R. Any solution R to equation (2.5) yields a Richardson number R e corres ponding to the fluxes V x ,V y and Q and not on z as ν 1 and ν 2 are independent on the depth variable z as well asq the turbulent viscosities. The Richardson number R e being known, steady-state profiles for velocity and density are obtained by integrating (2.4) with respect to z, taking into account the boundary conditions at z = -h:

(2.6)

u e (z) = u b + V x ρ a ρ 0 f 1 (R e ) (z + h) , v e (z) = v b + V y ρ a ρ 0 f 1 (R e ) (z + h) , ρ e (z) = ρ b + Q f 2 (R e ) (z + h) .
It remains to analyse the existence of solutions of equation (2.5). These solutions can be interpreted as the intersection of the curves

k (R) = (f 1 (R)) 2 f 2 (R) and h (R) = CR with C = - ρ 2 a (V 2 x + V 2 y ) gQρ 0
. The existence and the number of solutions are controlled by the constant C and then by the parameter

V 2 Q , V 2 = V 2
x + V 2 y , depending only on the surface fluxes. The graph of function k and h for Q < 0 and Q > 0 is plotted on Figures 1 and2 below when f 1 and f 2 in case of R-2-2-4 and R -2 -2. CR meets k at only one point for Q > 0 and Q < 0. Therefore it exists one unique equilibrium Richardson number R e whatever the values of the surface fluxes V x , V y and Q. In the case of the other models, we get several solutions. The proof is finished. Notice that in [START_REF] Bennis | Parametrization of the mixing layer: Comparison of four models depending on the richardson number[END_REF] we show that the most accurate model is R -2 -2 -4 from the physical and numerical viewpoint.

Linear stability of the equilibrium solutions

In this section we analyse the time evolution of a small perturbation of one of the equilibrium states (u e , v e , ρ e ) described in the previous section.

At initial time t = 0 we set (u 0 , v 0 , ρ 0 ) = (u e , v e , ρ e ) + (u ′ 0 , v ′ 0 , ρ ′ 0 ) and we denote by (u, v, ρ) = (u e , v e , ρ e ) + u ′ , v ′ , ρ ′ the solution of equations (2.1) at time t where (u e , v e , ρ e ) are solution to the steady-state system (2.3), and ν e 1 = f 1 (R e ) and ν e 2 = f 2 (R e ) are two positive constants. Introducing the new variables ψ = ∂ρ ∂z , θ = ∂u ∂z and β = ∂v ∂z , the Richardson number can be expressed as

R = - g ρ 0 ψ (θ 2 + β 2 ) = R (θ, β, ψ)
Applying the Taylor formula, we get The equations satisfied by the perturbation (u ′ , v ′ , ρ ′ ) are deduced from equations (2.1):

F = (θ -θ e ) ∂ν 1 ∂θ (θ e ,
(2.7)

     ∂u ′ ∂t - ∂ ∂z (ν 1 (θ, β, ψ) (θ e + θ ′ )) = 0, ∂v ′ ∂t - ∂ ∂z (ν 1 (θ, β, ψ) (β e + β ′ )) = 0, ∂ρ ′ ∂t - ∂ ∂z (ν 2 (θ, β, ψ) (ψ e + ψ ′ )) = 0.
We now replace ν 1 and ν 2 by expresions deuced from the Taylor's development and retain only the first order terms. The approximated equations for (u ′ , v ′ , ρ ′ ) then are (2.8)

                   ∂u ′ ∂t - ∂ ∂z ν e 1 + θ
            , V =   u ′ v ′ ρ ′   ,
Equations (2.8) can be written

(2.9) ∂V ∂t - ∂ ∂z A ∂V ∂z = ∂V ∂t -A ∂ 2 V ∂z 2 = 0.
Let (λ 1 , λ 2 , λ 3 ) be the eigenvalues of matrix A. Assuming the eigenvalues dinstincts,matrix A is equal to P -1 DP , where D is diagonal, and such that d 11 = λ 1 , d 22 = λ 2 and d 33 = λ 3 .

Set now W = P V. The vector W verifies the system ∂W ∂t -D ∂ 2 W ∂z 2 = 0, i.e.

(2.10)

∂w 1 ∂t -λ 1 ∂ 2 w 1 ∂z 2 = 0, ∂w 2 ∂t -λ 2 ∂ 2 w 2 ∂z 2 = 0, ∂w 3 ∂t -λ 3 ∂ 2 w 3 ∂z 2 = 0.
Stability of the equilibrium solution (u e , v e , ρ e ) means that any perturbation (u ′ 0 , v ′ 0 , ρ ′ 0 ) imposed at initial time t = 0 is damped as t → ∞. This is verified if the eigenvalues λ 1 , λ 2 , λ 3 are such that Re (λ 1 ) > 0, Re (λ 2 ) > 0 and Re (λ 3 ) > 0. These three conditions are equivalence to det (A) > 0, tr (A) > 0 and tr (Adj(A)) > 0. From these conditions, we build the graph below (see figure 1), obtained thanks an analytical computation (we skip the technical details here): The results are summarized in Figure 1.

The circle zone represents a zone where the solution is physically not valid. It is the case for the R-2-3 and R-2-1-3 formulation. The rectangular zone is a unstability zone. All formulations have a unstability zone. Nevertheles, one observes that for each model, mathematical stability holds for non negative R.

Conclusion

All the models have a steady-state solution, unique in the case of R-2-2-4. Each one is linearly stable for non negative R, which corresponds to physical stability. All these models present a linear unstable zone, located in a region where R is non positive. They all presents a linear stability zone for some non positive values of R, situation that can arise in real situation, as reported in [START_REF] Bennis | Parametrization of the mixing layer: Comparison of four models depending on the richardson number[END_REF] (physical unstability). All these models have been tested in [START_REF] Bennis | Parametrization of the mixing layer: Comparison of four models depending on the richardson number[END_REF]. The simulation confirms the existence of stable linear steady-state solutions and the ability of these models to describe a boundary mixing layer. However, the numerical study in [START_REF] Bennis | Parametrization of the mixing layer: Comparison of four models depending on the richardson number[END_REF] confirms that R -2 -2 -4 yields better numerical results.

  Figure A

  β e , ψ e ) + (ββ e ) ∂ν 1 ∂β (θ e , β e , ψ e ) + (ψψ e ) ∂ν 1 ∂ψ (θ e , β e , ψ e ) + • • • G = (θθ e ) ∂ν 2 ∂θ (θ e , β e , ψ e ) + (ββ e ) ∂ν 2 ∂β (θ e , β e , ψ e ) + (ψψ e ) ∂ν 2 ∂ψ (θ e , β e , ψ e ) + • • • We set for k = 1, 2 : F = ν 1 (θ, β, ψ)ν 1 (θ e , β e , ψ e ), G = ν 2 (θ, β, ψ)ν 2 (θ e , β e , ψ e ) and ν e k = ν k (θ e , β e , ψ e ), θ ′ = θθ e , β ′ = ββ e , ψ ′ = ψψ e , ∂ν k ∂θ e = ∂ν k ∂θ (θ e , β e , ψ e ) , ∂ν k ∂β e = ∂ν k ∂β (θ e , β e , ψ e ) , ∂ν k ∂ψ e = ∂ν k ∂ψ (θ e , β e , ψ e ) .

Figure 1 :

 1 Figure 1: Numerical stability

  e ∂ν 1 ∂θ

								e	θ ′ -	∂ ∂z	θ e ∂ν 1 ∂β	e	β ′ -	∂ ∂z	θ e ∂ν 1 ∂ψ	e	ψ ′ = 0,
	∂v ′ ∂t	-	∂ ∂z	β e ∂ν 1 ∂θ	e	θ ′ -		∂ ∂z		ν e 1 + β e ∂ν 1 ∂β	e	β ′ -	∂ ∂z	β e ∂ν 1 ∂ψ	e	ψ ′ = 0,
	∂ρ ′ ∂t	-	∂ ∂z	ψ e ∂ν 2 ∂θ	e	θ ′ -	∂ ∂z	ψ e ∂ν 2 ∂β	e	β ′ -	∂ ∂z	ν e 2 + ψ e ∂ν 2 ∂ψ	e	ψ ′ = 0.
	We set
			 	ν e 1 + θ e ∂ν 1 ∂θ	e	θ e ∂ν 1 ∂β	e			∂ψ θ e ∂ν 1	e
			
			
	A =	    	β e ∂ν 1 ∂θ	e			ν e 1 + β e ∂ν 1 ∂β	e	∂ψ β e ∂ν 1	e
			
			 	ψ e ∂ν 2 ∂θ	e			ψ e ∂ν 2 ∂β	e		ν e 2 + ψ e ∂ν 2 ∂ψ	e


Houssem Haddar
email: haddar@cmap.polytechnique.fr

Rainer Kress
email: kress@math.uni-goettingen.de

A conformal mapping algorithm for the Bernoulli free boundary value problem

 for the solution of a related inverse boundary value problem. For this we interpret the free boundary Γ0 as the unknown boundary in the inverse problem to construct Γ0 from the Dirichlet condition on Γ0 and Cauchy data on the known boundary Γ1. Our method for the Bernoulli problem iterates on the missing normal derivative on Γ1 by alternating between the application of the conformal mapping method for the inverse problem and solving a mixed Dirichlet-Neumann boundary value problem in D. We present the mathematical foundations of our algorithm and prove a convergence result. Some numerical examples will serve as proof of concept of our approach.

Introduction

Both free and inverse boundary value problems for partial differential equations are important contemporary topics in applied analysis. Typically in both of them the boundary of the solution domain consists of a known and an unknown part and the latter has to be determined from overdetermined Cauchy data. These Cauchy data are given on the known boundary part in the case of inverse problems and, to the contrary, on the unknown boundary part in the case of free boundary problems.

A classical example for a free boundary value problem is the Bernoulli problem from fluid dynamics. The velocity field V of an incompressible irrotational fluid flow satisfies div V = 0 and curl V = 0 and therefore, at least locally, there exists a harmonic function v, called velocity potential, such that V = grad v. In two dimensions, in addition to the potential also a stream function w can be introduced as a conjugate harmonic of v. Since the velocity V = grad v is orthogonal to the equipotential lines v = const and since the lines v = const and w = const of the two conjugate harmonic functions are orthogonal, the latter represent the streamlines of the potential flow.

In boundary value problems for such potential flows the boundaries will be either rigid, in which case they are known, or they will be free streamlines, i.e., free boundaries, in which case their shape is not known. On rigid boundaries the fluid flow is assumed to be tangential to the boundary, i.e., the rigid boundary is also a streamline. On a free boundary, to compensate for the fact that its shape is not known, a second boundary condition in addition to w = const has to be imposed. In cases where the free boundary is the interface between the fluid and a surrounding gas the free streamline is in equilibrium with the gas and this means that the pressure p of the fluid at the free streamline is equal to the constant pressure in the gas. Finally conservation of energy, i.e., Bernoulli's law 1 2 |V | 2 + p = const allows us to reformulate this condition into the form | grad w| = const in terms of the stream function w. This finally implies that the normal derivative of the stream function w must be constant on the free streamline.

For the concise formulation of the Bernoulli problem we assume that D 0 and D 1 are two simply connected bounded domains in R 2 with C 2 smooth boundaries Γ 0 := ∂D 0 and Γ 1 := ∂D 1 such that D 1 ⊂ D 0 and denote by D the doubly connected domain D := D 0 \ D 1 . The unit normal vector ν to both Γ 0 and Γ 1 is assumed to be directed into the complement of D. The Bernoulli free boundary value problem consists in determining the exterior boundary curve Γ 0 such that the unique solution w ∈ H

and w = 1 on Γ 1 (1.3) satisfies - ∂w ∂ν = λ on Γ 0 (1.4)
where Γ 1 is known and λ is a given positive constant. We note that by Hopf's lemma for harmonic functions (see e.g. [START_REF] Evans | Partial Differential Equations[END_REF]), the positivity of the constant λ is necessary for the existence of a solution to (1.1)-(1.4). Since the free boundary Γ 0 is the exterior component of ∂D, the above problem is referred to as the exterior Bernoulli free boundary value problem.

Together with the corresponding interior Bernoulli problem where the roles of the free and the known boundary component are interchanged, it has been extensively investigated (see e.g. [START_REF] Acker | On the geometric form of Bernoulli configurations[END_REF][START_REF] Alt | Existence and regularity for a minimum problem with free boundary[END_REF] and the literature referenced therein). Existence of a solution has been established by Beurling [5] more than fifty years ago. Contemporary existence proofs are based on variational methods in the sense of shape optimization. Analyticity of the free boundary has been shown by Lewy [START_REF] Lewy | A note on harmonic functions and a hydrodynamic problem[END_REF]. If Γ 1 is convex, uniqueness of the solution has been shown by Tepper [START_REF] Tepper | Free boundary problem[END_REF]. Counter examples show that convexity for Γ 1 is necessary for uniqueness (see [START_REF] Flucher | Bernoullis free-boundary problem, Qualitative theory and numerical approximation[END_REF]). Tepper [START_REF] Tepper | Free boundary problem[END_REF][START_REF] Tepper | On a free boundary problem, the starlike case[END_REF] has also shown that if Γ 1 is convex or starlike, then so is Γ 0 .

Bernoulli's free boundary problem also occurs in electro-and magnetostatics. For example, the optimal design problem to construct an insulation layer that minimizes the current leakage from coaxial cables leads to a Bernoulli problem. A further example is galvanization where in order to produce a layer of constant thickness it is required that the applied electric field has constant strength on the surface of the workpiece. For these examples we refer to [START_REF] Flucher | Bernoullis free-boundary problem, Qualitative theory and numerical approximation[END_REF] and the references therein.

The prominent numerical methods for solving free boundary value problems are so-called trial methods from shape optimization. The basic structure of such methods is to first make an initial guess for the free boundary. Then, firstly, the boundary value problem with one of the conditions on the free boundary omitted is solved and, secondly, the remaining boundary condition is used to update the free boundary. These two steps are iterated until both free boundary conditions are satisfied up to some specified accuracy (see [START_REF] Flucher | Bernoullis free-boundary problem, Qualitative theory and numerical approximation[END_REF][START_REF] Harbrecht | Stabilization of the trial method for the Bernoulli problem in case of prescribed Dirichlet data[END_REF][START_REF] Harbrecht | Improved trial methods for a class of generalized Bernoulli problems[END_REF] and further references therein).

Here, with the conformal mapping method in mind that was proposed and analyzed by Akduman, Haddar and Kress over the last decade (see [START_REF] Akduman | Electrostatic imaging via conformal mapping[END_REF][START_REF] Haddar | Conformal mappings and inverse boundary value problems[END_REF] among others) we proceed differently and involve the inverse boundary value problem to determine the unknown boundary Γ 0 from the homogeneous Dirichlet condition on Γ 0 and the Cauchy data on Γ 1 . For the solution of this inverse problem we introduce the annulus B bounded by two concentric circles C 0 with radius ρ > 1 and C 1 with radius one centered at the origin. In the sequel, we will identify the Euclidean space R 2 and the complex plane C in the usual way. By the Riemann conformal mapping theorem for doubly connected domains there exists a uniquely determined radius ρ > 1 and a holomorphic function Ψ that maps B bijectively onto D such that the boundaries C 0 and C 1 are mapped onto Γ 0 and Γ 1 , respectively, with all boundary curves in counter clockwise orientation. The function Ψ is unique up to a rotation of the annulus B. We parameterize Analogous to D, the unit normal vector ν to both C 0 and C 1 is directed into the complement of B.

The main ingredient of the conformal mapping method is a nonlinear differential equation for the boundary correspondence function ϕ in terms of the given Cauchy data on Γ 1 . In the first step of the algorithm this differential equation, in combination with an equation for the radius ρ, is solved by successive approximations. Then in the second step, knowing ρ and ϕ a Cauchy problem for Ψ is solved by a Laurent expansion and the unknown boundary is obtained as Γ 0 = Ψ(C 0 ). Thus we obtain a a parameterization of Γ 0 by

χ(t) := Ψ(ρe it ), t ∈ [0, 2π]. (1.6)
For the purpose of iteratively solving the Bernoulli problem, we make an initial guess for the unknown normal derivative g := ∂ ν w on Γ 1 and apply our conformal mapping method to construct an approximation for the unknown Γ 0 . Then we update g by the normal derivative on Γ 1 of the solution to the mixed Dirichlet-Neumann problem (1.2) and (1.4) and iterate these two steps. We will present the mathematical foundation of this algorithm and prove a convergence result. Numerical examples will confirm the feasibility of the method.

The plan of the paper is as follows. In Section 2 we will outline the basic properties of the conformal mapping algorithm for the inverse boundary value problem as required for our purpose. Then we will present the details of our proposed algorithm for the free boundary value problem in Section 3 followed by a convergence analysis in Section 4. The numerical examples in Section 5 serve as proof of concept.

To conclude this introduction we briefly sketch a further possibility to solve the Bernoulli problem with the aid of the conformal mapping Ψ. For this we note that the solution v := w • Ψ to the Dirichlet problem in B with boundary values v = 1 on C 1 and v = 0 on C 0 is given by

v(x) = 1 - ln |x| ln ρ , x ∈ B, (1.7)
and has normal derivatives

∂v ∂ν = - 1 ρ ln ρ on C 0 (1.8) and ∂v ∂ν = 1 ln ρ on C 1 . (1.9)
We denote by w a conjugate harmonic of w and set

v = w • Ψ, that is, v is a conjugate harmonic of v.
Then in view of (1.4) and (1.8) from the Cauchy-Riemann equations, both for w and w and for v and v we obtain that

-λ |χ (t)| = |χ (t)| ∂w ∂ν (χ(t)) = |χ (t)| ∂ w ∂s (χ(t)) = d dt w(χ(t)) = d dt v(ρe it ) = ρ ∂v ∂ν (ρe it ) = - 1 ln ρ for 0 ≤ t ≤ 2π, that is, |χ (t)| = 1 λ ln ρ , 0 ≤ t ≤ 2π.
(1.10) By Cauchy's integral formula, we can represent the holomorphic function Ψ by its boundary values χ and ϕ as

Ψ(z) = ρ 2π 2π 0 χ(τ ) e iτ ρ e iτ -z dτ - 1 2π 2π 0 γ(ϕ(τ )) e iτ e iτ -z dτ, z ∈ B.
From this, by the Sokhotski-Plemelj jump relations (see [START_REF] Kress | Integral Equations[END_REF]) for the Cauchy integral by letting z tend to C 0 and C 1 we obtain the two integral equations

χ(t) = 1 2π 2π 0 χ(τ ) e iτ e iτ -e it dτ - 1 2π 2π 0 γ(ϕ(τ )) e iτ e iτ -ρ e it dτ, 0 ≤ t ≤ 2π, and
γ(ϕ(t)) = ρ 2π 2π 0 χ(τ ) e iτ ρ e iτ -e it dτ - 1 2π 2π 0 γ(ϕ(τ )) e iτ e iτ -e it dτ, 0 ≤ t ≤ 2π,
for χ and ϕ. In our future research on the Bernoulli problem we intend to explore the possibility to solve these two nonlinear integral equations using the additional condition (1.10).

The conformal mapping method adjusted

We note that as opposed to our previous papers on the conformal mapping method the unknown boundary Γ 0 is the exterior part of the boundary ∂D.

With this reversion of the roles of the exterior and interior boundaries the conformal method solves the inverse boundary value problem to find the unknown boundary Γ 0 from the Cauchy data f = u| Γ1 and g = ∂ ν u| Γ1 on Γ 1 of a harmonic function u ∈ H 1 (D) satisfying u = 0 on Γ 0 . We proceed by a brief adjustment of the conformal mapping method for this special case. Proceeding analogous to the derivation of (1.10), from the Neumann conditions ∂ ν u = g on Γ 1 and (1.9) we obtain that

g(γ(ϕ(t))) d dt γ(ϕ(t)) = 1 ln ρ (2.1)
for 0 ≤ t ≤ 2π. Since the boundary correspondence function ϕ is strictly monotonically increasing, this implies

ϕ = 1 ln ρ |γ • ϕ| g • γ • ϕ (2.2)
together with the boundary conditions

ϕ(0) = 0, ϕ(2π) = 2π. (2.3)
Integrating (2.1) over the interval [0, 2π] we obtain

1 ln ρ = 1 2π Γ1 g ds (2.4)
as formula for the radius ρ. We note that in the general case for arbitrary Dirichlet data f the formula for the radius corresponding to (2.4) contains the unknown ϕ and the differential equation corresponding to (2.2) is nonlocal (see [START_REF] Akduman | Electrostatic imaging via conformal mapping[END_REF][START_REF] Haddar | Conformal mappings and inverse boundary value problems[END_REF]).

To deal with the two boundary conditions (2.3) for the differential equation of order one, integrating (2.2) we obtain

1 2π 2π 0 1 ln ρ |γ • ϕ| g • γ • ϕ dt -1 = 0.
Subtracting this from (2.2) and inserting (2.4) for ρ leads to the modified differential equation

ϕ = Γ1 g ds 2π |γ • ϕ| g • γ • ϕ -Γ1 g ds 4π 2 2π 0 1 |γ • ϕ| g • γ • ϕ dt + 1 (2.5)
that the boundary correspondence function ϕ has to satisfy together with the initial condition ϕ(0) = 0.

(

Obviously, each solution ϕ to the initial value problem (2.5)-(2.6) automatically also satisfies the second boundary condition ϕ(2π) = 2π. Provided g is positive on all of Γ 1 , the classical Picard-Lindelöf theorem guarantees uniqueness and existence for the initial value problem (2.5)-(2.6) and this solution can be obtained by successive approximations. By Hopf's lemma combined with the maximum-minimum principle, the required positivity is satisfyed for the Neumann trace

g = ∂ ν u| Γ1 of a harmonic function in H 1 (D) satisfying u = 1 on Γ 1 and -∂ ν u = λ on Γ 0 .
From the equivalent form of the initial value problem (2.5)-(2.6) as a nonlinear Volterra integral equation, by the implicit function theorem it can be seen that the operator T : H 1 (Γ 1 ) → H 1 [0, 2π] given by T : g → ϕ which maps the Neumann data g onto the solution ϕ of (2.5)-(2.6) is Fréchet differentiable and that the Fréchet derivative dT (g; h) of T at g in the direction h is given by the unique solution to the linear differential equation

[dT (g; h)] = α(g)dT (g; h) + β(g; h) - 1 2π 2π 0 {α(g)dT (g; h) + β(g; h)} dt (2.7)
satisfying the initial condition dT (g; h)(0) = 0.

(2.8)

Here we abbreviated

α(g) := - (γ • ϕ) • (γ • ϕ) Γ1 g ds 2π |γ • ϕ| 3 g • γ • ϕ - (g • γ) • ϕ Γ1 g ds 2π |γ • ϕ| [g • γ • ϕ] 2 (2.9)
and

β(g; h) := - h • γ • ϕ Γ1 g ds 2π |γ • ϕ| [g • γ • ϕ] 2 + Γ1 h ds 2π |γ • ϕ| g • γ • ϕ (2.10)
with ϕ = T g.

For the Cauchy problem part of the conformal mapping method we expand γ•ϕ in a Fourier series with Fourier coefficients a n and obtain from the definition (1.5) of ϕ that

Ψ(e it ) = ∞ n=-∞ a n e int , t ∈ [0, 2π].
(2.11)

By the Laurent series for Ψ, in view of (1.6), this implies

χ(t) := ∞ n=-∞ a n ρ n e int , t ∈ [0, 2π], (2.12)
for the parametrization of Γ 0 . In order to numerically cope with the instability of this parametrization with respect to numerical errors in the Fourier coefficients a n for n > 0 we incorporate a regularization by truncating the series in the form

χ N (t) := N +1 n=-N +1 a n ρ n e int , t ∈ [0, 2π], (2.13)
where N ∈ N serves as a regularization parameter. This may be interpreted as a singular value cut-off. We note that also a Tikhonov type regularization can be used (see [START_REF] Akduman | Electrostatic imaging via conformal mapping[END_REF][START_REF] Haddar | Conformal mappings and inverse boundary value problems[END_REF]). We denote by R, R N :

H 1 [0, 2π] → C 2 [0, 2π
] the linear operators that map the series (2.11) onto the series (2.12) and (2.13) for χ and χ N , respectively. Then we can summarize the conformal mapping method by the operators

M := RγT and M N = R N γT from H 1 (Γ 1 ) into C 2 [0,
2π] that map the given normal derivative g onto the parameterizations χ and χ N , respectively. By the chain rule, the Fréchet derivative of M is given by

dM (g; h) := R[γ • (T g) dT (g; h)] -2πρ Γ1 h ds Γ1 g ds 2 R (γ • T g), (2.14)
where R denotes the derivative of the series (2.12) with respect to ρ. For dM N (g; h) we have an analogous expression with R replaced by R N .

The free boundary problem algorithm

The main idea of our iterative method for the free boundary value problem is to apply the conformal mapping algorithm of the previous section with the Cauchy pair f = 1 and g = w| Γ1 for the solution of (1.1)- (1.3). For this we make an initial guess for the unknown normal derivative g and then iterate by alternating between the conformal mapping method for obtaining an approximation for Γ 0 and updating g by solving the mixed Dirichlet-Neumann problem (1.1), (1.2) and (1.4).

The Algorithm 1 Choose an initial guess g.

2 Apply the conformal mapping method to obtain an approximation Γ 0 by (2.12), or for the sake of stability by (2.13).

3 After introducing the operator F : C 2 [0, 2π] → H 1 (Γ 1 ) that maps the parameterization χ of Γ 0 onto the normal derivative g of the solution to (3.1)-(3.2) on Γ 1 we note that the iteration scheme without regularization in the second step can be expressed in terms of the operator

A := F M (3.3) as g n+1 = Ag n , n = 0, 1, 2, . . . . (3.4)
Later on we will find in necessary for convergence to modify (3.4) by applying a relaxation

g n+1 = ωg n + (1 -ω)Ag n , n = 0, 1, 2, . . . , (3.5)
with a relaxation parameter ω ∈ [0, 1) adjusted to the value of λ for the Neumann boundary condition on Γ 1 .

In Step 3 alternatively one could solve the corresponding Dirichlet-Neumann problem in the annulus B via Fourier series. However, this leads to an update for g • ϕ whereas Step 2 requires g on Γ 1 in order to obtain an update for ϕ via the iteration step for (2.5)-(2.6). Using the boundary integral equation in D itself avoids the issue of a numerical inversion of ϕ to obtain g| Γ1 from g • ϕ.

In principle, for Step 2 the conformal mapping method could be replaced by any other solution method for the inverse boundary value problem. We have chosen the conformal mapping method for two reasons. Firstly, complex analysis methods have played an important role in the history of free boundary problems for the Laplace equation (see e.g. [START_REF] Garabedian | Partial Differential Equations[END_REF]). Secondly, and more importantly, for this particular choice for Step 2 we were able to analyze the convergence based on our techniques from the conformal mapping method for the inverse problem.

For completeness we briefly outline the use of boundary integral equations for the solution of (3.1)-(3.2). In terms of the fundamental solution

Φ(x, y) := 1 2π ln 1 |x -y| , x = y,
to the Laplace equation, Green's integral formula applied to a solution w of (3.1)-(3.2) reads

w(x) = Γ1 Φ(x, •)g ds - Γ0 λΦ(x, •) + ∂Φ(x, •) ∂ν ω ds (3.6)
for x ∈ D where we have set g := ∂νw| Γ1 and ω := w| Γ0 and used the fact that the double-layer potential with constant density on Γ 1 vanishes in D. Letting x approach Γ 0 and Γ 1 we obtain the system of two integral equations

ω 2 + Γ0 ∂Φ(x, •) ∂ν ω ds - Γ1 Φ(x, •)g ds = -λ Γ0 Φ(x, •) ds, x ∈ Γ 0 , (3.7) and Γ0 ∂Φ(x, •) ∂ν ω ds - Γ1 Φ(x, •)g ds = -λ Γ0 Φ(x, •) ds -1, x ∈ Γ 1 , (3.8)
for the two unknowns ω and g. Without loss of generality, we assume that there exists a point z in D 1 such that |x -z| = 1 for all x ∈ D 1 . Then Theorem 3.16 in [START_REF] Kirsch | An Introduction to the Mathematical Theory of Inverse Problems[END_REF] (see also [START_REF] Kress | Integral Equations[END_REF]Theorem 7.38]) guarantees that the single-layer operator from H -1/2 (Γ 1 ) into H 1/2 (Γ 1 ) is bijective. Note that by a change of variables the setting always can be scaled such that this condition is satisfied. For simplicity, we have chosen to just work with the single-layer operator and point out that alternatively in order to obtain injectivity we could modify it by adding an additional term as in [START_REF] Kress | Integral Equations[END_REF]Theorem 7.41].

We note that plugging in ω = 0 for the zero Dirichlet condition on the free boundary makes (3.7)-(3.8) a system of integral equations for the unknown free boundary Γ 0 and the unknown normal derivative ∂ ν w| Γ1 as a slip variable. With this interpretation, the (3.7)-(3.8) serves as basis for an integral equation method for the solution of the Bernoulli problem due to Trefftz [START_REF] Trefftz | Über die Kontraktion kreisförmiger Flüssigkeitsstrahlen[END_REF][START_REF] Wegmann | Trefftz' integral equation method for free boundary problems of potential theory[END_REF].

Using the jump relations and taking care of the logarithmic behavior of the single-layer potentials as in the proof of Theorem 7.38 in [START_REF] Kress | Integral Equations[END_REF] it can be shown that the homogeneous form of (3.7)-(3.8) only has the trivial solution and that for any solution ω and g of the inhomogeneous equation the function w defined in D by the right hand side of (3.6) provides a solution of (3.1)-(3.2). Finally, after the regularization of equation (3.8) with the inverse of the single-layer operator on Γ 1 the Riesz theory for compact operators can be employed to show existence of a unique solution ω ∈ H 1/2 (Γ 0 ) and g ∈ H -1/2 (Γ 1 ). Using the regularity of the right hand side and the mapping properties of the involved boundary operators higher regularity of the solution can be achieved.

For the approximate solution of (3.7)-(3.8) we recommend the efficient Nyström and collocation method for boundary integral equations for harmonic functions based on trigonometric approximations as described in [START_REF] Kress | Integral Equations[END_REF].

A convergence result

We begin our considerations on the convergence of the above iteration scheme by some preliminary observations on the simple case when Γ 1 is the unit circle. In this case the solution to the Bernoulli free boundary problem is given by a concentric circle Γ 0 and from (1.7) we observe that its radius is given by the nonlinear equation

1 ρ ln ρ = λ (4.1)
and that g = 1 ln ρ .

We start the iteration with a constant and assume the n-th iteration to be given by a constant g n . Then from (2.4) we have

ρ n = exp 1 g n
and the solution in Step 3 becomes w(x) = 1 -ρ n λ ln |x|. Hence, the update equation is given by

g n+1 = λ exp 1 g n . (4.3)
The function q : t → λ exp(1/t) has derivative q (1/ ln ρ) = -ln ρ. Therefore, provided g > 1, that is, ρ < e, the iterations (4.3) converge locally, i.e., we have convergence if the initial guess is not too far from the correct g. From (4.1) we observe that ρ < e is satisfied if λ > e -1 , since clearly ρ > e implies that λ < e -1 . For a more general convergence result, in principle, we need to estimate the norm of the derivative of the iteration operator A as given by (3.3) via dA(g; h) = dF (M g; dM (g; h)).

(4.4)

Therefore, in addition to the Fréchet derivative of M given by (2.7) and (2.14) we also need the Fréchet derivative of F . According to Allaire [START_REF] Allaire | Shape Optimization by the Homogenization Method[END_REF]Corollary 6.36] the Fréchet derivative of F at χ in the direction η is given by dF (χ; η) = ∂ ν v| Γ1 where v is the unique solution to the mixed Dirichlet-Neumann problem in D with boundary condition v = 0 on Γ 1 (

and

∂v ∂ν = -κλ ν • η + ∂ ∂s ν • η ∂w ∂s on Γ 0 (4.6)
in terms of the solution w of (3.1)-(3.2). Here, s is the arc length on Γ 0 and κ denotes the curvature of Γ 0 . Since deriving sharp estimates on the norm of dA based on the representations of dM and dF appears out of reach in the general case, we confine ourselves to the simple situation when Γ 1 is the unit circle with the canonical parameterization γ(t) = e it , that is, when D = B. In this case g is a constant and the exact values for ρ and g are given by (4.1) and (4.2). We proceed by evaluating the Fréchet derivative at the exact solution g and consequently also the exact ϕ given by ϕ(t) = t. Since γ • γ = 0 and g = 0 from (2.9) we obtain that α(g) = 0. From (2.10), using (4.2), we find that

β(g; h) = ln ρ -h • γ • ϕ + 1 2π Γ1 h ds .
With this, solving (2.7) and (2.8) for h = 1 leads to dT (g; 1) = 0 and from (2.14) we obtain

dM (g; 1)(t) = -ρ(ln ρ) 2 e it . (4.7)
Solving (2.7) and (2.8) for the trigonometric monomials c n (t) := cos nt and s n (t) := sin nt for n ∈ N we obtain

dT (g; c n • γ -1 ) = - ln ρ n s n and dT (g; s n • γ -1 ) = ln ρ n (c n -1) whence dM (g; c n • γ -1 )(t) = ln ρ 2n ρ 1-n e (1-n)it -ρ 1+n e (1+n)it (4.8)
and

dM (g; s n • γ -1 )(t) = i ln ρ 2n ρ 1-n e (1-n)it + ρ 1+n e (1+n)it -2ρe it (4.9)
follow by (2.12) and (2.14). Now we need to consider the mixed Dirichlet-Neumann problem (4.5)-(4.6) in D for η given by the right hand sides of (4.8) and (4.9) and interpreted as a vector in R 2 for the inner product with the normal vector ν on Γ 0 . Since the solution w of (3.1)-(3.2) is constant and κλ = λ/ρ = 1/ρ 2 ln ρ, the Neumann condition (4.6) reduces to

∂v ∂ν = - 1 ρ 2 ln ρ ν • η on Γ 0 .
This becomes ∂v ∂ν = ln ρ ρ on Γ 0 for η = dM (g; 1),

∂v ∂ν (ρe it ) = 1 2nρ (ρ -n -ρ n ) cos nt, t ∈ [0, 2π],
for η = dM (g; c n • γ -1 ), and

∂v ∂ν (ρe it ) = 1 2nρ (ρ -n -ρ n ) sin nt, t ∈ [0, 2π],
for η = dM (g; s n • γ -1 ). The correponding solutions to (4.5)-(4.6) are

v(x) = ln ρ ln |x|, 1 ≤ |x| ≤ ρ, for h = 1, v(re it ) = - 1 2n 2 ρ 2n -1 ρ 2n + 1 (r n -r -n ) cos nt, 1 ≤ r ≤ ρ, t ∈ [0, 2π], for h = c n • γ -1 and v(re it ) = - 1 2n 2 ρ 2n -1 ρ 2n + 1 (r n -r -n ) sin nt, 1 ≤ r ≤ ρ, t ∈ [0, 2π], for h = s n • γ -1 . Its normal derivatives on Γ 1 become ∂v ∂ν = -ln ρ on Γ 1 for h = 1, ∂v ∂ν (e it ) = 1 n ρ 2n -1 ρ 2n + 1 cos nt, t ∈ [0, 2π],
for h = c n • γ -1 , and

∂v ∂ν (e it ) = 1 n ρ 2n -1 ρ 2n + 1 sin nt, t ∈ [0, 2π], for h = s n • γ -1 .
In view of (4.4) and M g = ργ, this finally implies that dA(g; 1) = µ 0 with the eigenvalue

µ 0 = -ln ρ,
a result that we also can obtain more easily from (4.3), and further

dA(g; c n • γ -1 )) = µ n c n • γ -1 and dA(g; s n • γ -1 )) = µ n s n • γ -1
with the eigenvalues

µ n = 1 n ρ 2n -1 ρ 2n + 1
for n ∈ N.

For n ≥ 2 consider the function p : [1, ∞) → R defined by

p(t) := n(t n + 1) - t + 1 t -1 (t n -1).
With the aid of a geometric sum it can be seen that

p(t) = n(t n + 1) -1 -2 n-1 k=1 t k -t n
and for t ≥ 1 this implies

p (t) = n(n -1)t n-1 -2 n-1 k=1 kt k-1 ≥ n(n -1)t n-1 -2t n-1 n-1 k=1 k ≥ 0.
Therefore p is nonnegative on [0, ∞) since p(1) = 0. Setting t = ρ 2 , from this we now can deduce that sup

n∈N µ n = µ 1 .
From the above eigenfunctions and eigenvalues we now can conclude that for Γ 1 a circle the Fréchet derivative dA(g; •) : H 1 (Γ 1 ) → H 1 (Γ 1 ) at the exact g is a self adjoint compact operator with spectral radius

σ = max ln ρ, ρ 2 -1 ρ 2 + 1
which is less than one if ρ < e. As already observed above, λ > e -1 implies that ρ < e as consequence of (4.1).

If we consider the regularized iteration operatorA N where R is replaced by the truncated Cauchy operator R N , then µ n , n = 0, 1, 2, . . . , N, are also eigenvalues of dA N (g; •) and the remaining eigenvalues become zero. Hence, the norms of dA N (g; •) and dA(g; •) coincide.

From the presentation of the Fréchet derivative of A we can observe that all its components except the ill-posed Cauchy operator R depend continuously on g and Γ 1 . Hence, the Fréchet derivative of the regularized operator A N depends continuously on g and Γ 1 . Therefore, for fixed N , if Γ 1 is sufficiently close to a circle, there exists a ball centered at the correct g such that within this ball the Fréchet derivative of A N has norm less than one. Therefore, making use of the mean value theorem for nonlinear operators, the following theorem can be proven.

Theorem 4.1 Assume that Γ 1 is sufficiently close to a circle and λ > e -1 . Then the iterations g n+1 = A N g n , n = 0, 1, 2, . . . , converge provided the initial guess g 0 is sufficiently close to the exact normal derivative g for the circle.

The eigenvalues ω + (1 -ω)µ n , n ∈ N , of the Fréchet derivative of the relaxation operator B(ω) := ωI + (1 -ω)A (for the circle at the exact solution) has absolute value smaller than one for ω ∈ [0, 1) and, depending on λ, the parameter ω can be chosen such that ω -(1 -ω) ln ρ also has absolute value smaller than one. Therefore the restriction λ > e -1 of our theorem is not essential, since it can be circumvented by the relaxation iterations g n+1 := B(ω)g n .

As we will illustrate by the numerical examples in the following final section, convergence also occurs for boundaries Γ 1 with shapes quite different from circles.

Numerical examples

As proof of concept rather than a documentation of a fully developed code, in this final section we present some numerical constructions for free boundaries. For the known rigid boundary Γ 1 we choose the parameterization γ(t) = (1.5 cos t + ε cos 2t, 2 sin t), 0 ≤ t ≤ 2π, (5.1) which satisfies our geometric assumption that ensures invertibility of the singlelayer operator. For the shape parameter ε in (5.1) and the parameter λ for the Neumann condition we have chosen the values ε = 0.2, 0.4, 0.6 and λ = 0.5, 2.

In the figures the rigid boundary Γ 1 is given as dotted (green) curve and the free boundary Γ 0 as full (red) curve.

For the solution of the initial value problem (2.5)-(2.6) we approximated the iterations ϕ n by a trigonometric polynomial of degree N in the form For the outer iterations we used the relaxation factor ω = 0.5 and they were stopped when the normal derivative on Γ 1 changed less than δ = 0.00001 in the maximum norm. The increase of the deviation form the zero Dirichlet condition on the free boundary for geometries with larger deviations from an circular annulus, obviously, is due to the fact that our algorithm does not impose this Dirichlet condition explicitly. Nevertheless, the examples indicate satisfactory performance of our algorithm although we violated the principle not to destroy stability of a problem by incorporating an ill-posed component in its numerical solution.

This research was initiated while R.K. was visiting INRIA Saclay Ile de France and Ecole Polytechnique, Palaiseau. The hospitality and the support are gratefully acknowledged.

  1 (D) to the Dirichlet problem for the Laplace equation ∆w = 0 in D (1.1) with boundary values w = 0 on Γ 0 (1.2)

Γ 1 =

 1 {γ(t) : t ∈ [0, 2π)} with a continuously differentiable 2π periodic function γ : R → C such that γ| [0,2π) is injective and fix Ψ uniquely B by prescribing Ψ(1) = γ(0). Then we define a boundary correspondence function ϕ : [0, 2π] → [0, 2π] by setting ϕ(t) := γ -1 (Ψ(e it )), t ∈ [0, 2π].(1.5)

1 (

 1 Solve the mixed Dirichlet-Neumann problem for ∆w = 0 in D with boundary conditions w = 1 on Γ integral equation method, i.e., an integral equation for the unknowns ∂ ν w| Γ1 and w| Γ0 . If g -∂ ν w| Γ1 ∞ < δ for a given tolerance δ terminate the iteration. Otherwise update g := ∂ ν w| Γ1 and go back to Step 2.

  2π |γ • ϕ n | g • γ • ϕ n n | g • γ • ϕ n dt + 1at the M equidistant points t m := 2πm/M , m = 1, 2, . . . , M, to obtain a linear system of M +1 equations, including the initial condition, for the 2N +1 Fourier coefficients of ϕ n+1 which we solved by least squares. We used trigonometric polynomials of degree N = 8 and M = 62 collocation points and kept the number N also for the truncation of the Laurent series. For the number Z i of inner iterations to solve (2.5)-(2.6) we made the a priori choice Z i = 5. The iterations where started with g 0 = 2λ resulting in ϕ 0 (t) = t.For the solution of the Dirichlet-Neumann problem we numerically solved the integral equations (3.7)-(3.8) by Nyström and collocation methods based on trigonometric polynomial interpolation as described in[START_REF] Kress | Integral Equations[END_REF]. We again used M = 64 equidistant collocation points for each of the two equations. The results did not change significantly for M = 128 collocation points.

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: Free boundary for ε = 0.2 and λ = 2 (left) and λ = 0.5 (right)

Figure 5 . 3 : 1 :

 531 Figure 5.3: Free boundary for ε = 0.6 and λ = 2 (left) and λ = 0.5 (right) ε λ Z o w ∞ 0.2 0.5 12 0.00003 0.2 2 17 0.00092 0.4 0.5 13 0.00016 0.4 2 19 0.00350 0.6 0.5 16 0.00881 0.6 2 21 0.00671


Houssem Haddar
email: houssem.haddar@inria.fr

Jacob D Rezac
email: jrezac@udel.edu

A quasi-backscattering problem for inverse acoustic scattering in the Born regime

We consider the problem of detecting three-dimensional inclusions from quasi-backscattering far field data generated by an incident field of timeharmonic fixed frequency plane waves modeled with the Born approximation. We assume only partial far field data is known and use a sampling-type method to reconstruct small obstacles and extended spherical obstacles. In particular, at the location of a device transmitting an incident wave, we assume far field data is collected only along a line extending a short distance from the transmitting device. Several numerical examples are provided to demonstrate the effectiveness of the approach.

Introduction

In this paper we propose a data collection geometry in which to frame the inverse scattering problem of locating unknown obstacles from far field measurements of time-harmonic scattering data. The measurement geometry, which we call a quasibackscattering set-up, requires less data than traditional multi-static configurations. We demonstrate that the data collected can be used to locate inhomogeneities whose physical properties are such that the Born approximation applies. In particular, we are able to image a two-dimensional projection of the location of a small obstacle by checking if a test function which corresponds to a point in R 2 belongs to the range of a measurable operator. Combining several projections then allows us to identify the location of the small inclusions in R 3 . We also show how this algorithm can be extended to the case of extended spherical inclusions.

Our reconstruction algorithm uses ideas from sampling-type methods [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF][START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF][START_REF] Collino | Numerical and analytical studies of the linear sampling method in electromagnetic inverse scattering problems[END_REF][START_REF] Colton | The linear sampling method in inverse electromagnetic scattering theory[END_REF]] and MUSIC-type (MUltiple-SIgnal-Classification) algorithms for small obstacles as in [START_REF] Cheney | The linear sampling method and the MUSIC algorithm[END_REF][START_REF] Devaney | Super-resolution processing of multi-static data using time reversal and MUSIC[END_REF][START_REF] Kirsch | The MUSIC algorithm and the factorization method in inverse scattering theory for inhomogeneous media[END_REF][START_REF] Ammari | Two numerical methods for recovering small inclusion from the scattering amplitude at a fixed frequency[END_REF][START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF][START_REF] Haddar | Mdimagh Identification of small inclusions from multistatic data using the reciprocity gap concept[END_REF]. However, the data setting and analysis are different. Typical scattering experiments collect either backscattering data, in which one device acts as both transmitter and receiver, or multi-static data, in which many transmitting/receiving devices surround the inhomogeneity. The algorithms in the above-mentioned references work only for multi-static data. The inversion scheme we describe in this paper makes use of a type of quasi-backscattering data; one device acts as a transmitter and a line of receivers extends in one-dimension a small distance from the transmitter. Figure 1 demonstrates the difference between a usual multi-static Figure 1: A Comparison of multi-static data (left) and quasi-backscattering data (right). Red circles correspond to device locations. The thick blue line in the right figure denotes where the quasi-backscattering set-up is moved and where each transmitting device is located. set-up and the quasi-backscattering geometry. As the figure demonstrates, the quasibackscattering geometry requires significantly less data than the multi-static geometry, which can be beneficial in practical applications. While the proposed scheme can be applied to other types of problem, we focus here on the inverse acoustic scattering problem.

The outline of our paper is as follows: in Section 2, the direct scattering problem is formulated and the quasi-backscattering data setting is explicated. In Section 3, we introduce and analyze the inversion procedure capable of identifying two-dimensional projections of small objects' locations. We then extend the algorithm to the case of extended spherical inclusions. In Section 4, extensive numerical experimentations are presented in order to show the performance of this new algorithm. We end this section by explaining how one can obtain three-dimensional locations from two-dimensional projections.

Direct scattering problem

We begin by discussing the mathematical formulation for the problem of acoustic incident plane waves scattering against inhomogeneous media in three-dimensions. These details can be found in many places, e.g. [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], and we include them here for completeness.

Assume a plane wave incident field with a fixed wave number k is generated in the area of an inhomogeneity. Such an incident field is described by u i (x, d) = e ik d•x for x ∈ R 3 and d ∈ S 2 where S 2 is the boundary of the unit ball in R 3 and where • indicates a normalized vector. The total field u(x)

satisfies ∆u + k 2 n(x)u = 0 in R 3 , (1a) u(x) = u i (x, d) + u s (x), ( 1b
) lim r→∞ r ∂u s ∂r -iku s = 0, (1c)
where u s (x) is the scattered field, r = |x| is the Euclidean magnitude of x, n(x) is the refractive index of the inhomogeneous medium which is a bounded function with non-negative imaginary part, and (1c) is the Sommerfeld radiation condition which holds uniformly with respect to x = x/|x|. For wave numbers such that Imk ≥ 0 and refractive indices which are bounded with compact support, it is known that (1a)-(1c) has a unique solution in H 1 loc (R 3 ). This problem has been studied extensively and more information about the related direct and inverse problems can be found in, e.g., [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF][START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF].

In what follows, define the contrast function as m(x) = 1 -n(x) and assume m(x) is non-zero only on a compact set D ⊂ R 3 which contains the inhomogeneity. The Lippmann-Schwinger equation,

u(x) = e ikx• d -k 2 D m(y)Φ(x, y)u(y) dy, (2)
gives an exact expression for the unique solution to (1a)-(1c) where

Φ(x, y) = 1 4π e ik|x-y| |x -y| , x = y
is the fundamental solution to the Helmholtz equation in R gives the Born approximation

u B (x) = e ikx• d -k 2 D m(y)Φ(x, y)e ik d•y dy. (3)
The inverse problem in which we are interested is to find information about D given data about the asymptotic behavior of u s B (x), the Born approximation to the scattered field. We are able to explicitly characterize the asymptotic behavior of the scattered field through the Sommerfeld radiation condition, (1c). Specifically,

u s B (x) = e ik|x| |x| u B ∞ (x, d) + O 1 |x| 2 , |x| → ∞
where u ∞ B (x, d) is the Born approximation to what is known as the far field pattern. Using (3),

u s B (x) = -k 2 D m(y)Φ(x, y)e ik d•y dy and we conclude that u B ∞ (x, d) = - k 2 4π D e ik( d-x)•y m(y) dy, x ∈ S 2 . (4)

Forward problem for quasi-backscattering data

The above derivations have not fixed the measurement geometry. We now restrict x and d to the quasi-backscattering experimental set-up. In what follows, let x = -d + ηê where η ∈ [-δ, δ] for a small constant δ and ê ∈ S 2 is a fixed unit vector which is orthogonal to d. The traditional backscattering set-up corresponds with δ = 0. Using the orthogonality of d with ê and the fact that both are unit vectors, a Taylor expansion about η = 0 yields x = -d + ηê

-d + ηê = -d + ηê 1 + η 2 = -d + ηê + O(η 2 ).
As such, we chose x in this way as an approximation to x = -d+ηê -d+ηê up to O(η 2 ). For this reason we continue to use the notation x, although it is no longer normalized. Substituting this choice of x into (4) gives

u B ∞ (-d + ηê, d) = - k 2 4π D e 2ik d•y e -ikηê•y m(y) dy, d ∈ S 1 (ê) (5)
where S 1 (ê) := { d ∈ S 2 ; d • ê = 0}. Following the typical approach of sampling methods in inverse scattering problems, we introduce the quasi-backscattering far field operator, F :

L 2 ([-δ, δ]) → L 2 (S 1 (ê)
) which we will use extensively in solving the inverse problem. In particular, F is defined as

(F g)( d) = δ -δ u B ∞ (-d + ηê, d)g(η) dη, d ∈ S 1 (ê). ( 6
)

Quasi-backscattering inverse problem

We now turn our attention to the inverse problem of reconstructing the location of inhomogeneities from the quasi-backscattering far field data. We first consider the case of obstacles which are small compared to the wavelength of the incident wave and which are sufficiently far from one another. In Section 3.4 , we use the analysis for this case as the basis for finding the centers of extended spherical obstacles. The key result of this section is Theorem 1 which will allow us to locate obstacles by testing if a specific function is in the range of the quasi-backscattering far field operator.

Inverse problem for small obstacles

Assume there are M obstacles with supports described by D j ⊂ R 3 , j = 1, . . . , M , embedded in a homogeneous background. Let the contrast be defined by the weighted sum of characteristic functions m(x) = M j=1 m j 1 Dj where m j are constants. If D j = z j + R j Ω j are small obstacles centered at a point z j ∈ R 3 with size and shape described by R j and Ω j respectively, then using [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] we obtain that up to O(max(R j ) 4 ) error terms,

u B ∞ (-d + ηê, d) M j=1 τ j e 2ik d•zj e -ikηe•zj , d ∈ S 1 (ê), η ∈ [-δ, δ]. (7)
Here, τ j = -k 2 4π m j R 3 j |Ω j |, where |Ω j | indicates the volume of Ω j , are constants related to the strength of each scatterer. Combining ( 6) and ( 7) reduces the quasibackscattering operator to

(F g)( d) = δ -δ u B ∞ (-d + ηê, d)g(η) dη = M j=1 τ j e 2ik d•zj δ -δ e -ikηe•zj g(η) dη.
To further simplify the far field operator, we write each obstacle's location in terms of its components parallel to ê and perpendicular to ê. For a fixed ê, we write z j = Π ⊥ (z j ) + L(z j )ê, j = 1, . . . , M where Π ⊥ maps onto the plane orthogonal to ê and where L isolates the component of a vector which is parallel to ê. For example, if ê = (0, 0, 1) and z 1 = (1, 2, 3), we would have Π ⊥ (z 1 ) = (1, 2, 0) and L(z 1 ) = 3. Note that for the sake of conciseness we will sometimes treat Π ⊥ (z) as a vector in R 2 .

Decomposing the locations of obstacles in this way, we can write the far field operator as

(F g)( d) = M j=1 τ j e 2ik d•Π ⊥ (zj ) δ -δ e -ikL(z j )η g(η) dη. ( 8
)
Since F can be computed from the measurable far field pattern data, we use it to solve the inverse scattering problem. Indeed, Theorem 1 gives conditions under which we can relate the range of F to the location of a small obstacle. Such a characterization is typical for sampling-type methods such as the linear sampling or factorization schemes, as well as the MUSIC algorithm. Before stating this theorem, we prove two short lemmas which are required for its proof.

Lemma 1. Assume ê ∈ S 2 is fixed and let z j ∈ R 3 , j = 1, . . . , M be distinct points whose components in the direction of ê differ (i.e., L(z

i ) = L(z j ), i = j). Then A = {η → e -ikL(zj )η , j = 1, . . . , M } is a linearly independent sequences of functions for η ∈ [-δ, δ].
Proof. We would like to show that the Wronskian matrix of A, denoted by W , is non-singular. A short calculation shows that det

(W ) = c(η)det(V ) where c(η) = exp -ikη M j=1 L(z j ) is a function which never vanishes and V (i,j) = ω j-1 i for ω = -ikL(z i ).
Since V is a Vandermonde matrix, it has a non-zero determinant so long as ω i = ω j for each i = j, which is true by the assumption on L(z i ), i = 1, . . . , M . Lemma 2. Assume ê ∈ S 2 is fixed, let z j ∈ R 3 , j = 1, . . . , M be distinct points, and let z * ∈ R 3 be any point perpendicular to ê and distinct from each

Π ⊥ (z j ). Then B = { d → e 2ik d•z , z = z * , Π ⊥ (z 1 ), . . . , Π ⊥ (z M )} is a linearly independent sequences of functions for d ∈ S 1 (ê).
The proof of this lemma follows the idea of Theorem 4.1 in [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF]. This theorem of Kirsch and Grinberg implies that the above can also be proven for a finite number of dj , η j ∈ S 1 with a similar but more technical argument.

Proof. To show that B is linearly independent, assume

β 0 e 2ik d•z * + M j=1 β j e 2ik d•Π ⊥ (zj ) = 0 for d ∈ S 1 (ê).
The left-hand-side of the above equation is, up to a constant multiple, the far field pattern of the function

x → β 0 Φ(x, z * ) + M j=1 β j Φ(x, Π ⊥ (z j )) where Φ(x, z) = i 4 H (1) 0 (2k|x -z|), x = z (9)
is the (radiating) fundamental solution of the Helmholtz equation in R 2 with wave number 2k and H

(1) 0 is a Hankel function. As such, since the far field pattern vanishes, Rellich's lemma and unique continuation show that

β 0 Φ(x, z * ) + M j=1 β j Φ(x, Π ⊥ (z j )) = 0 for x / ∈ {z * , Π ⊥ (z 1 ), . . . , Π ⊥ (z M )}.
Taking the limit as x approaches each of z * and Π ⊥ (z j ), j = 1, . . . , M shows immediately that B is a linearly independent sequence of functions for each d ∈ S 1 (ê).

With these lemmas in hand, we are ready to prove the key theorem for small obstacles.

Theorem 1. Assume ê ∈ S 2 is fixed and d ∈ S 1 (ê). Let z j ∈ R 3 for j = 1, . . . , M and let z ∈ R 3 be orthogonal to ê. If the components of each z j parallel to ê are not equal (i.e., L(z

i ) = L(z j ), i = j), then φ z ( d) = e 2ik d•z ∈ R(F ) if and only if z ∈ {Π ⊥ (z j ), j = 1, . . . , M }.
Proof. Let z / ∈ {Π ⊥ (z j ), j = 1, . . . , M } be orthogonal to ê. Assume by contradiction that there exists some g(η)

∈ L 2 ([-δ, δ]) such that (F g)( d) = e 2ik d•z . From the definition of F , this would imply e 2ik d•z = M j=1 c j e 2ik d•Π ⊥ (zj ) ,
where c j = τ j δ -δ e -ikL(zj )η g(η) dη are constants. However, this is a contradiction with the linear independence of

{ d → e 2ik d•Π ⊥ (ζ) , ζ = z, Π ⊥ (z 1 ), . . . , Π ⊥ (z M )}, which shows that if φ z ∈ R(F ) then z ∈ {Π ⊥ (z j ), j = 1, . . . , M }.
To prove the second half of the theorem, assume L(z i ) = L(z j ), i = j and z ∈ {Π ⊥ (z j ), j = 1, . . . , M }. We will show that φ z ∈ N (F * ) ⊥ = R(F ) which gives the result since F is a finite rank operator with closed range. A short calculation gives that

(F * h)(η) = M j=1 τ j e ikL(zj )η S 1 (ê) h( d)e -2ik d•Π ⊥ (zj ) ds( d). If h ∈ N (F * ), then M j=1 τ j e ikL(zj )η S 1 (ê) h( d)e -2ik d•Π ⊥ (zj ) ds( d) = 0.
The linear independence of {e ikL(zj )η , j = 1, . . . , M } proven in Lemma 1 gives that for each j = 1, . . . , M , 0 =

S 1 (ê) h( d)e -2ik d•Π ⊥ (zj ) ds( d) = h( d), φ Π ⊥ (zj ) ( d) L 2 (S 1 (ê))
where (•, •) L 2 (S 1 (ê)) indicates the inner-product on L 2 (S 1 (ê)). As such, φ z ∈ N (F * ) ⊥ for each z ∈ {Π ⊥ (z j ), j = 1, . . . , M }, which gives the result.

The proof of Theorem 1 in fact implies a slightly stronger result.

Corollary 1. With no restrictions on L(z j ) and the same hypotheses on ê as in Theorem 1, if φ z ∈ R(F ) then Π ⊥ (z) ∈ {Π ⊥ (z j ), j = 1, . . . , M }.

Another corollary to Theorem 1 is that, for the appropriate restrictions on Π ⊥ (z j ), P φ z = 0 if and only if z ∈ {Π ⊥ (z j ), j = 1, . . . , M } where P : L 2 (S 1 (ê)) → R(F ) ⊥ is the orthogonal projection onto the orthogonal complement of the range of F . This suggests that the function I(z) = P φ z -1 for each z perpendicular to ê within a region of interest will be large when z is near Π ⊥ (z j ), j = 1, . . . , M and small otherwise. This is exactly the MUSIC-type algorithm which we will use to locate the centers of small objects.

To construct the imaging function

I(z), let (u k , σ k , v k ), k = 1, 2, .
. . be the singular system for F where the left singular functions are u k ∈ L 2 (S 1 (ê)) and the right singular functions are v k ∈ L 2 ([-δ, δ]). Since R(F ) is spanned by the left singular functions u k which correspond to singular values σ k = 0, we can write

I(z) = ∞ k=r+1 (φ z , u k ) L 2 (S 1 (ê)) 2 -1 , ( 10
)
where r is the number of non-zero singular values. Numerical results showing that I(z) is large near obstacles are given in Section 4.

A second range test for three-dimensional reconstructions

Assume that the range test described above has been performed so that {Π ⊥ (z j ), j = 1, . . . , M } are known. From Theorem 1,

φ Π ⊥ (z k ) ( d) ∈ R(F ) for a given k = 1, . . . , M . As such, there is a g Π ⊥ (z k ) ∈ L 2 ([-δ, δ]) such that (F g Π ⊥ (z k ) )( d) = M j=1 c j e 2ik d•Π ⊥ (zj ) = e 2ik d•Π ⊥ (z k )
where, as before, c j = τ j δ -δ g Π ⊥ (z k ) (η)e -ikL(z k )η dη. By linear independence, c j = τ j δ jk , where δ jk is the Kronecker delta function. This suggests a second indicator function which can be used to find L(z j ) when Π ⊥ (z j ) are already known. Formally,

J Π ⊥ (z k ) (z) = δ -δ g Π ⊥ (z k ) (η)e -ikL(z)η dη -1 (11)
is arbitrarily large when z = L(z j ), j = k. This argument is formal and we have no guarantee that J Π ⊥ (z k )(z) is small away from z = z j , j = k. Nevertheless, in the numerical examples in Section 4.2 below, J Π ⊥ (z k ) indicates the location of L(z j ), as expected when Π ⊥ (z j ) are known accurately and g z (η) is calculated using Tikhonov regularization and the Morozov discrepancy principle. As the numerical simulations will demonstrate, however, calculating L(z j ) in this manner is not robust to noise.

Inverse problem for coplanar small obstacles

Due to the hypotheses on Theorem 1, the algorithm outlined above does not necessarily locate an object in the case that L(z i ) = L(z j ) for some i = j ∈ 1, . . . , M . This problem can be easily alleviated: in all proofs we have assumed a fixed ê ∈ S 2 . Since L(z j ) is a function of ê, we can perform multiple quasi-backscattering experiments with different ê directions to solve the problem. Indeed, we recommend this for purely geometric reasons as well. Since the quasi-backscattering technique gives only twodimensional projections of the locations of scatterers, two obstacles which lie on top of each other with respect to ê (i.e., Π ⊥ (z i ) = Π ⊥ (z j ) but L(z i ) = L(z j )) will appear as the same obstacle in the reconstruction. Multiple experiments corresponding to different ê directions helps to fix this problem as well. In Section 4.3, we outline a technique for using data from multiple experiments with different ê directions to reconstruct obstacles in three-dimensions.

Before continuing, we note that the algorithm outlined above does not necessarily identify obstacles if L(z i ) = L(z j ) for all i = j. In particular, under these conditions, we show below that there is no obvious reason which suggests that I(z) will be arbitrarily large at z = Π ⊥ (z j ), j = 1, . . . , M . Indeed, the numerical simulations in Section 4 indicate that the reconstruction of co-planar obstacles is sensitive to noise.

Assume that L(z i ) = L(z j ) for each i, j = 1, . . . , M . Since we can shift the origin with a change of variables, we set each L(z i ) = 0 without loss of generality. In this case, the far field operator becomes

(F g)( d) = M j=1 τ j e 2ik d•Π ⊥ (zj ) δ -δ g(η) dη, d ∈ S 1 (ê),
and Fubini's Theorem gives that for η ∈ [-δ, δ],

(F * h)(η) = M j=1 τ j S 1 (ê) e -2ik d•Π ⊥ (zj ) h( d) ds( d) =   h( d), M j=1 τ j φ Π ⊥ (zj )   L 2 (S 1 (ê)) . Let h( d) = u k ( d) for a fixed k where u k ( d) is a left singular function of F corresponding to a singular value σ k = 0. Since u k ∈ N (F * ),   u k ( d), M j=1 τ j φ Π ⊥ (zj )   L 2 (S 1 (ê)) = (F * u k )(η) = 0.
However, we cannot conclude from the above equation that

u k ( d), φ Π ⊥ (zj ) L 2 (S 1 (ê))
= 0, j = 1, . . . , M.

Inverse problem for extended obstacles

We now adapt the arguments given in the previous section to the problem of finding extended obstacles. We show that, for small δ, the arguments given in Theorem 1 apply directly to locating the center of extended spherical obstacles. While the spherical nature of the extended obstacles does not seem to be required, it is not clear that we can uncover more information than the location of the center of these obstacles.

In this section, assume there are M obstacles D j again of the form D j = z j +R j Ω j , where z j are the obstacles' center, Ω j their shape, and R j their size. Now, however, assume each Ω j = B(0; 1) is a ball centered at zero of radius one and that R j is of similar size as the wavelength or larger. Assume, as before, that the contrast is defined by m = M j=1 m j 1 Dj where m j are constants. We will begin our discussion by calculating u B ∞ (x, d) under these assumptions. From (4) we have

u B ∞ (x, d) = - k 2 4π M j=1 m j Dj e ik( d-x)•y dy = - k 2 4π M j=1
m j e ik( d-x)•zj B(0;Rj ) e ik( d-x)•y dy.

To simplify this expression into a more useful one, we state the following lemma.

It is known that S 2 e irx•ŷ ds(ŷ) = 4πj 0 (rx) where j 0 is the spherical Bessel function of order zero [START_REF] Martin | Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles[END_REF]. Since j 0 (x) = sin(x)

x , an integration-by-parts gives the result.

We are interested in the above result for x = k( dx) and R = R j . With this choice of parameters,

u B ∞ (x, d) = -k 2 M j=1 m j e ik( d-x)•zj sin(kR j | d -x|) -kR j | d -x| cos(kR j | d -x|) (k| d -x|) 3 .
Returning to the quasi-backscattering approach and letting x = -d+êη, a Taylor expansion about η = 0 gives that sin(kR

j | d -x|) -kR j | d -x| cos(kR j | d -x|) (k| d -x|) 3 = sin(2kR j ) -2kR j cos(2kR j ) (2k) 3 + O(η 2 ), j = 1, . . . , M.
If we again define the quasi-backscattering far field operator as (F g)

( d) = δ -δ u B ∞ ( d, η)g(η) dη, we find (F g)( d) = M j=1 τ L j e 2ik d•zj δ -δ e -ikηê•zj g(η) dη + o(δ 2 ) = M j=1 τ L j e 2ik d•Π ⊥ (zj ) δ -δ e -ikL(z j )η g(η) dη + o(δ 2 )
where τ L j = -m j (sin(2kR j ) -2kR j cos(2kR j ))/8k. Here, the asymptotic analysis follows from an application of Cauchy-Schwarz. Up to constants and o(δ 2 ), the quasi-backscattering operator for extended spheres is identical to [START_REF] Cheney | The linear sampling method and the MUSIC algorithm[END_REF], the quasibackscattering operator for small obstacles. As such, using the same technique described in Section 3.1, we can find two-dimensional projections of the centers of extended spherical obstacles. Note, incidentally, that τ L j = -1 3 k 2 m j R 3 j + O(R 5 j ), which matches the expression for τ j used in the case of small spheres.

Numerical experiments

In this section, we give numerical results demonstrating the effectiveness of the above technique. In all experiments, we approximate obstacles by spheres with a small radius. Specifically, the radius for each obstacle is 1/500 units. We will use simulated forward data which is corrupted by random noise. Using the formula given in (4) we simulate u ∞ (x, d) using numerical integration. Numerically integrating [START_REF] Colton | The linear sampling method in inverse electromagnetic scattering theory[END_REF] gives a discrete representation of the far field matrix, F ij , which is corrupted by F ij (1 + γξ) where ξ is a uniform random variable in [-1, 1] and γ is a constant related to the level of noise. To calculate the indicator function, we compute the singular value decomposition of F = U SV * and use U to calculate a discrete version of [START_REF] Devaney | Super-resolution processing of multi-static data using time reversal and MUSIC[END_REF]. The approximate imaging function is regularized by computing with all but the first ten singular vectors (i.e., r = 9 in [START_REF] Devaney | Super-resolution processing of multi-static data using time reversal and MUSIC[END_REF]). In all examples, we take k = 15 to be the wave number. Other parameters are given for each experiment.

The experimental parameters discussed above merit a few comments. The first is related to our use of ten singular vectors in reconstructions. Typically when using a MUSIC-type algorithm, the number of singular vectors is related to the number of unknown obstacles, which is estimated by the numerical rank of the far field operator. In our numerical experiments we have found such a technique to be sensitive to added noise. As such, we took the number of singular vectors as an upper bound of the number of obstacles. This method for picking the number of singular values could result in false positives in reconstructions; though seemingly unlikely in practice, there is no guarantee that this modification of the indicator function in [START_REF] Devaney | Super-resolution processing of multi-static data using time reversal and MUSIC[END_REF] will be large only at obstacle locations. A more robust technique for picking the number of singular vectors would be useful in practice. Nevertheless, in the experiments below, results do not change noticeably when using the same number of singular vectors as there are obstacles.

The second comment is related to the relatively-high wave number used in these experiments. In the case of extended obstacles, low wave numbers are used to ensure that the Born approximation of the far field is valid. However, because we assume our objects are very small (a radius of 1/500 units), we are justified in using a higher wave number.

We present three types of numerical inversions. In Section 4.1, we show twodimensional projections of small obstacles. In Section 4.2, we generate the third unknown coordinate assuming the first two are known. Finally, we use multiple êdirections to generate full three-dimensional reconstructions for small obstacles in Section 4.3.

Two-dimensional projections of small obstacles

We give several numerical examples in this section which help demonstrate both the strengths and weaknesses of the quasi-backscattering technique. In all reconstructions, darker colors correspond to higher values of the imaging function which correspond with the predicted locations of the obstacles. A small red circle in each picture corresponds to the true location of each obstacle. Note that the size of the dark areas near obstacle locations do not correspond to an estimate of obstacle size, but are merely an artifact of the way in which reconstructions are displayed. In all reconstructions, we use 80 2 sampling points uniformly chosen in the unit-square. In all experiments we use 95 incident directions and for each incident direction we use 95 locations for x between -δ and δ. We call these points between -δ and δ the observation points. This is a large number of both incident directions and observation points and, indeed, acceptable results are achievable with far fewer. However, we prefer to focus these experiments on the affect geometric and physical parameters have on reconstructions.

The first example, given in Figures 2 and3, shows the algorithm differentiating between multiple small obstacles, added one at a time. The obstacles are located at z 1 = (-0.25, -0.25, -0.5), z 2 = (0.25, 0.25, -0.25), z 3 = (0.25, -0.25, 0.25), and z 4 = (-0.25, 0.25, 0.5). For the next example, we show the affect of δ on reconstructions. While the motivation for the quasi-backscattering set-up comes from a Taylor expansion about η = 0 (and hence small δ), the experiments in Figure 4 show that in the presence of noise, the reconstruction technique is not stable for too small of δ, in particular when many obstacles are present. All three figures have z 1 = (0, -0.5, 0.25), z 2 = (0, 0.5, -0.75), z 3 = (0.5, 0, -0.25), and z 4 = (-0.5, 0, 0.75). In the experiment in Figure 5 we show the resolution achievable by the quasibackscattering technique. Often, inversion schemes based on the Born approximation or a Fourier transform are limited to a half-wavelength resolution. Though we do not show this rigorously, the numerical example in Figure 5 suggests such a limitation for the quasi-backscattering technique. Indeed, we see that the method is unable to differentiate between obstacles once they are within half a wavelength of each other. In this case, the technique gives a large range of possible locations, containing the true centers of the obstacles. In this experiment, there is a constant 0.2 unit distance between the z-coordinate of the obstacles. The final two experiments of this type show the need to take multiple experiments with different ê directions when the underlying geometry of the obstacles is complicated. Figure 10 shows reconstructions from three different ê directions of three small obstacles which would form an approximate "L"-shape if they were connected with straight lines. In particular, z 1 = (-0.25, -0.25, -0.25), z 2 = (0.25, -0.24, -0.25), and z 3 = (-0.25, -0.26, 0.25). Due to the geometry of the obstacles, taking ê = (0, 0, 1) or ê = (1, 0, 0) only gives reconstructions of two of the three obstacles. By taking ê = (0, 1, 0), however, we are able to find all three obstacles.

Finally, we apply the quasi-backscattering algorithm to the reconstruction of coplanar obstacles -that is, obstacles which violate the assumptions in Theorem 1. As Figure 11 shows, in the absence of noise, reconstructions are acceptable. However, under the addition of noise, the reconstructions become less clean. Changing ê so that L(z i ) = L(z j ) results in more acceptable reconstructions. The figures are located at z 1 = (0.75, 0.75, 0.25) and z 2 = (-0.25, -0.25, 0.25).

Reconstruction of third coordinate

We now show two reconstructions of the third coordinate of a small obstacle, assuming the other two coordinates are known. We use the indicator function given by [START_REF] Haddar | Mdimagh Identification of small inclusions from multistatic data using the reciprocity gap concept[END_REF] where g z (η) is calculated using Tikhonov regularization plus the Morozov discrepancy principle. In both reconstructions, we take δ = π/50 and 377 observation points and incident directions. Though this is an unrealistically-large number of observation points and incident directions, we will show that the indicator function is still sensitive to noise. In both reconstructions, ê = (0, 0, 1) so that we are generating the zcoordinate in a typical Cartesian plane. For this reason, we explore another technique for three-dimensional reconstructions in Section 4.3 below.

In Figure 12, let z 1 = (-0.24, -0.24, -0.75), z 2 = (0.26, -0.24, 0), and z 3 = (0.26, 0.26, 0.75). Adding no noise and assuming the two-dimensional projections of z j , j = 1, 2, 3 are known exactly, the figure demonstrates we are able to construct L(z j ) under ideal circumstances. We consider a more realistic scenario in Figure 13.

Let z 1 = (-0.25, -0.25, -0.75), z 2 = (0.25, -0.25, 0), and z 3 = (0.25, 0.25, 0.75). However, we have added 0.1% noise and assume we guess Π ⊥ (z 1 ) = (-0.24, -0.24), Π ⊥ (z 2 ) = (0.26, -0.24), and Π ⊥ (z 3 ) = (0.26, 0.26). We see that even under small perturbations, the accuracy of the reconstructions is dramatically decreased.

Three-dimensional reconstructions

The inversion schemes described above do a good job locating obstacles within the twodimensional plane perpendicular to the selected ê direction. Given data from multiple experiments with multiple ê directions, we are better able to find the full threedimensional coordinates of an obstacle or set of obstacles. As discussed above, there are many scenarios in which reconstructing obstacles with multiple ê is encouraged. In this section, we show that multiple ê directions can be used to calculate threedimensional reconstructions of obstacle locations.

The creation of three-dimensional images from a selection of two-dimensional projections has been thoroughly studied in the image processing literature and we do not attempt to use state-of-the-art techniques here. Instead, we perform multiple quasi-backscattering experiments on the same obstacle set-up, interpolate the results from each experiment onto a fixed sampling grid, and average the results. We regularize our results for each ê before computing the averaged result. In particular, we apply a total variation minimization algorithm (see [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]) which emphasizes changes in gradient and hence sharpens edges. We next locally normalize each twodimensional projection over a 5 × 5 grid of sampling points to further sharpen edges. After these regularization steps are performed, we average on a sampling grid as described.

We compute forward and inverse data for this section as we did in Section 3.1. Here, however, we vary ê. Specifically, we take 30 values of ê from a circle in the XY -plane. The results are given as three-dimensional contour plots of the imaging function. The contour which is plotted is α max I(z) where α is a value between 0 and 1.

In Figure 14, we demonstrate the techniques described above to compute threedimensional object reconstructions. In particular, we consider three small objects located at the points z 1 = (-0.5, -0.5, -0, 5), z 2 = (0.5, -0.5, -0.5), and z 3 = (-0.5, -0.5, 0.5). Notice that these are in a geometry which forms an "L"-shape, as in Figure 11. As demonstrated above, when we use two-dimensional projection techniques, we require at least three ê directions to locate all objects for such a geometry. By taking more ê directions, however, we are able to give a full threedimensional image of the geometry.

Lemma 3 .r 2 S 2 e

 322 For a constant R > 0 and any two vectors x, y ∈ R 3 , B(0;R) e ix•y dy = 4π |x| 3 (sin(R|x|) -R|x| cos(R|x|)) . Proof. Under the change of coordinates y → rŷ where r = |y|, B(0;R) e ix•y dy = R 0 irx•ŷ ds(ŷ) dr.

Figure 2 :

 2 Figure 2: Symmetric addition of new objects, δ = π/50 (no noise).

Figure 3 :

 3 Figure 3: Symmetric addition of new objects, δ = π/50 (approximately 1% noise).

Figure 4 :

 4 Figure 4: Decreasing the quasi-parameter δ. Figures have δ = π (top-left), δ = π/50 (top-right), and δ = π/100 (bottom). Approximately 5% noise.

Figure 5 :

 5 Figure 5: Two objects moving closer to each other, δ = π/100 (no noise). Thick bar at bottom corresponds to half of wavelength.

Figure 10 :

 10 Figure 10: An L-Shaped geometry which requires 3 views to see all obstacles, δ = π/30 (approximately 1% noise).

Figure 11 :

 11 Figure 11: Results for co-planar obstacles, δ = π/60. Figures on left are noise-free and figures on right have approximately 1% noise.

JFigure 12 :

 12 Figure12: When J Π ⊥ (zj ) is computed, peaks appear for i = j. Red circles show exact location of L(z j ).

J

  ζ (z), ζ = Π ⊥ (z3)

Figure 13 :

 13 Figure 13: Reconstructions lose accuracy in the presence of 0.1% noise. Red circles show exact location of L(z j ).

Figure 14 :

 14 Figure 14: Three-dimensional noise-free reconstructions of point obstacles based on multiple experiments. We take 95 incident directions, 95 observation points, and use δ = π/60. The top figure is noise-free while the bottom figure has approximately 1% noise. In both figures, we display isovalues of 0.6 times the maximum value of the imaging function.

Acknowledgement

This research was initiated while JDR was visiting École Polytechnique where he was supported by INRIA-DeFI. The hospitality of that institution is gratefully acknowledged. Fioralba Cakoni is also thanked for her many helpful suggestions.