A D Agaltsov

T Hohage

R G Novikov

An iterative approach to monochromatic phaseless inverse scattering

Keywords: Inverse scattering problems, phaseless inverse scattering, Schrödinger equation AMS subject classification: 35J10, 35R30, 65N21, 81U40, 78A46

This paper is concerned with the inverse problem to recover a compactly supported Schrödinger potential given the differential scattering cross section, i.e. the modulus, but not the phase of the scattering amplitude. To compensate for the missing phase information we assume additional measurements of the differential cross section in the presence of known background objects. We propose an iterative scheme for the numerical solution of this problem and prove that it converges globally of arbitrarily high order depending on the smoothness of the unknown potential as the energy tends to infinity. At fixed energy, however, the proposed iteration does not converge to the true solution even for exact data. Nevertheless, numerical experiments show that it yields remarkably accurate approximations with small computational effort even for moderate energies. At small noise levels it may be worth to improve these approximations by a few steps of a locally convergent iterative regularization method, and we demonstrate to which extent this reduces the reconstruction error.

Introduction

In quantum mechanics the interaction of an elementary particle at fixed energy E > 0 with a macroscopic object contained in a bounded domain D is described by the Schrödinger equation

-∆ψ + v(x)ψ = Eψ, x ∈ R d .
(1.1a)

Here ∆ is the standard Laplacian in x, and the potential function v is assumed to satisfy

v ∈ L ∞ (R d ), supp v ⊂ D ⊂ R d , d ≥ 2.
(1.1b)

Equation (1.1a) can be also considered as the Helmholtz equation of acoustic and electrodynamic wave propagation at fixed frequency.

For equation (1.1a) we consider the classical scattering solutions ψ + (•, k) of the form ψ + (x, k) = e ikx + ψ s (x, k) with a plane incident field e ikx such that k ∈ R d , |k| 2 = E, and a scattered field ψ s (•, k) satisfying the Sommerfeld radiation condition (1.2)

1 Max-Planck Institute for Solar Systems Research, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany (agaltsov@mps.mpg.de). 2 Institute for Numerical and Applied Mathematics, University of Göttingen, Lotzestr. [START_REF] Hohage | Inverse problems with Poisson data: statistical regularization theory, applications and algorithms[END_REF][START_REF] Isayev | Exponential instability in the inverse scattering problem on the energy interval[END_REF][START_REF] Isayev | New global stability estimates for monochromatic inverse acoustic scattering[END_REF]37083 Göttingen, Germany (hohage@math.uni-goettingen.de) and Max-Planck Institute for Solar Systems Research.

3 CMAP, Ecole Polytechnique, CNRS, Université Paris-Saclay, 91128 Palaiseau, France; and IEPT RAS, 117997 Moscow, Russia (novikov@cmap.polytechnique.fr).

with a function f E called the scattering amplitude or far field pattern at energy E. There are different conventions for the choice of the constant c(d, |k|). The one above leads to the following simple asymptotic relation between the scattering amplitude and the inverse Fourier transform of v, see, e.g., [START_REF] Faddeev | Uniqueness of the solution of the inverse scattering problem[END_REF][START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF]:

v(k -l) = f E (k, l) + O(E -1 2 ), E → +∞, (1.3) v(p) := (2π) -d R d e ipx v(x) dx, p ∈ R d (1.4)
|f E (k, l)| 2 is known as the differential scattering cross section for equation (1.1a). In quantum mechanics this quantity describes the probability density of scattering of particle with initial impulse k into direction l/|l| = k/|k|, see, for example, [10, Chapter 1, Section 6]. Typically, the differential cross section is the only measurable quantity whereas the phase of the scattering amplitude cannot be determined directly by physical experiments. The problem of finding v from |f E | 2 is known as the phaseless inverse scattering problem for equation (1.1a). Whereas the inverse scattering problem with phase information for equation (1.1a), i.e. the problem of finding v from f E , has been studied intensively for a long time (see [START_REF] Alexeenko | Solution of the three-dimensional acoustical inverse scattering problem. The modified Novikov algorithm[END_REF][START_REF] Barceló | Numerical approximation of the potential in the two-dimesional inverse scattering problem[END_REF][START_REF] Burov | Multifrequency generalization of the Novikov algorithm for the two-dimensional inverse scattering problem[END_REF][START_REF] Chadan | Inverse problems in quantum scattering theory[END_REF][START_REF] Eskin | Lectures on linear partial differential equations[END_REF][START_REF] Faddeev | Uniqueness of the solution of the inverse scattering problem[END_REF][START_REF] Grinevich | The scattering transform for the two-dimensional Schrödinger operator with a potential that decreases at infinity at fixed nonzero energy[END_REF][START_REF] Hähner | New stability estimates for the inverse acoustic inhomogeneous medium problem and applications[END_REF][START_REF] Isayev | New global stability estimates for monochromatic inverse acoustic scattering[END_REF][START_REF] Isayev | Exponential instability in the inverse scattering problem on the energy interval[END_REF][START_REF] Novikov | Multidimensional inverse spectral problem for the equation -∆ψ + (v(x) -Eu(x))ψ = 0[END_REF][START_REF] Novikov | Rapidly converging approximation in inverse quantum scattering in dimension 2[END_REF][START_REF] Novikov | The ∂-approach to approximate inverse scattering at fixed energy in three dimensions[END_REF][START_REF] Novikov | Approximate Lipschitz stability for non-overdetermined inverse scattering at fixed energy[END_REF][START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF] and references therein), much less studies were performed in the phaseless case (see [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF][START_REF] Chadan | Inverse problems in quantum scattering theory[END_REF][START_REF] Klibanov | Phaseless inverse scattering problems in three dimensions[END_REF][START_REF] Klibanov | Reconstruction procedures for two inverse scattering problems without the phase information[END_REF][START_REF] Romanov | Phaseless inverse problems that use wave interference[END_REF][START_REF] Novikov | Inverse scattering without phase information[END_REF][START_REF] Novikov | Formulas for phase recovering from phaseless scattering data at fixed frequency[END_REF][START_REF] Novikov | Explicit formulas and global uniqueness for phaseless inverse scattering in multidimensions[END_REF]).

It is well known that the phaseless scattering data |f E | 2 does not determine v uniquely even if |f E | 2 is given completely for all E > 0; see, e.g., [START_REF] Novikov | Explicit formulas and global uniqueness for phaseless inverse scattering in multidimensions[END_REF]. In the present work we continue studies of [START_REF] Novikov | Explicit formulas and global uniqueness for phaseless inverse scattering in multidimensions[END_REF][START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] assuming additional measurements of the following form: For the unknown v satisfying (1.1b) we consider additional a priori known background scatterers w 1 , . . . , w L such that w l ∈ L ∞ (R d ), supp w l ⊂ Ω l , Ω l is an open bounded domain in R d , w l = 0, w l = w l if l = l (in L ∞ (R d )), (1.5) where l, l ∈ {1, . . . , L}. In practice, we also typically have Ω l ∩ D = ∅, but this property will not be needed in our analysis. We set

S = {|f E | 2 , |f E,1 | 2 , . . . , |f E,L | 2 }, (1.6)
where f E is the scattering amplitude for v at energy E, and f E,1 , . . . , f E,L are the scattering amplitudes for v 1 , . . . , v L , where

v l = v + w l , l = 1, . . . , L. (1.7)
One can see that S consists of the phaseless scattering data |f E | 2 , |f E,1 | 2 , . . . , |f E,L | 2 measured sequentially, first, for the unknown scatterer v and then for the unknown scatterer v in the presence of known scatterer w l disjoint from v for l = 1, . . . , L. We consider the following inverse scattering problem without phase information for equation (1.1a):

Problem 1. Reconstruct coefficient v from the phaseless scattering data S for some appropriate background scatterers w 1 , . . . , w L .

In this paper we propose an iterative approach to Problem 1 with iterates u (j) E , j = 1, 2, . . . and prove error bounds of the form

u (j) E -v L ∞ = O E -αj (1.8)
with exponents α j tending to ∞ as j → ∞ for infinitely smooth potentials v.

For the inverse scattering problem with phase information such a substantial improvement of the Born approximation, which serve as first iterate u 1 (see (1.3)) has been obtained in [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF], and first numerical tests were reported in [START_REF] Barceló | Numerical approximation of the potential in the two-dimesional inverse scattering problem[END_REF].

Studies on Problem 1 in dimension d = 1 for L = 1 were started in [START_REF] Aktosun | Inverse problems on the line without phase information[END_REF], where phaseless scattering data was considered for all E > 0. Note also that a phaseless optical imaging in the presence of known background objects was considered, in particular, in [START_REF] Govyadinov | Phaseless three-dimensional optical nanoimaging[END_REF]. Studies on Problem 1 in dimension d ≥ 2 were started in [START_REF] Novikov | Explicit formulas and global uniqueness for phaseless inverse scattering in multidimensions[END_REF] and continued recently in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]. The key result of [START_REF] Novikov | Explicit formulas and global uniqueness for phaseless inverse scattering in multidimensions[END_REF] consists in a proper extension of formula (1.3) for the Fourier transform v of v to the phaseless case of Problem 1, d ≥ 2, which will be detailed in Section 3.1. The main results of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] consist in proper extensions of formula (2.8) in the configuration space to the case of Problem 1 for d ≥ 2; see also Section 3.1. However, the convergence of the approximations to v as E → +∞ in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] is slow, in particular, the exponent α in (1.8) is always ≤ 1 2 . In addition, our theoretical iterative monochromatic reconstructions for Problem 1 are illustrated numerically in Section 4.

Iterative inversion with phase information 2.1 Inverse scattering with phase information

Recall that the scattering amplitude f E is defined on the set

M E = (k, l) ∈ R d × R d : |k| 2 = |l| 2 = E . (2.1)
In view of (1.3) we assume that the scattering amplitude, and later the differential cross section is defined on some subset

M m E ⊂ M E such that the function Φ : M m E → B d 2 √ E Φ(k, l) := k -l (2.2)
is surjective. Here and in the following B d r denotes the closed ball

B d r = p ∈ R d : |p| ≤ r , r > 0. (2.3)
For d ≥ 2 we may construct a d-dimensional subset M m E ⊂ M E such that Φ is bijective as follows: Let us choose a piece-wise continuous function γ :

B d 2 √ E → R d such that |γ(p)| = 1 and γ(p)p = 0 for all p ∈ B d 2 √ E and set M m E = p 2 + E -p 2 4 1/2 γ(p), -p 2 + E -p 2 4 1/2 γ(p) : p ∈ B d 2 √ E . (2.4)
To use the Born approximation (1.3) and its refinements if Φ is not injective, we average over the set Φ-1 (p). To this end we assume that for all p ∈ B d

√

E the set Φ-1 (p) is a piecewise smooth manifold of size Φ-1 (p) and define the averaging operator

Φ : L 1 (M m E ) → L 1 B d 2 √ E , (Φf )(p) := 1 | Φ-1 (p)| Φ-1 (p) f E (k, l) d(k, l). (2.5)
Using this mapping we can define an approximation u E to v on D by

u E (x) := B d 2 √ E e -ipx (Φf )(p) dp, x ∈ D. (2.6)
Let W n,1 (R d ) denote the Sobolev space of n-times smooth functions in the sense of L 1 (R d ):

W n,1 (R d ) := u ∈ L 1 (R d ) : u n,1 < ∞ with u n,1 := max |J|≤n ∂ |J| u ∂x J L 1 (R d ) , n ∈ N ∪ {0}. (2.7) If v ∈ W n,1 (R d ), n > d, in addition to the initial assumptions (1.1b), then u E satisfies the error bound u E -v L ∞ (D) = O(E -α ) as E → +∞ with α = n -d 2n (2.8)
for d ≥ 2; see, for example, [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF]. Essential improvements of the approximation u E in (2.6) were achieved in [START_REF] Novikov | Rapidly converging approximation in inverse quantum scattering in dimension 2[END_REF][START_REF] Novikov | The ∂-approach to approximate inverse scattering at fixed energy in three dimensions[END_REF][START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF]. In particular, formula (2.8) was principally improved in [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF] by constructing iteratively nonlinear approximate reconstructions u

(j) E such that u (1)
E = u E and u (j) E -v L ∞ (D) = O(E -αj ) with α j = n-d 2d 1 -n-d n j , j ≥ 1, (2.9
)

as E → +∞ for d ≥ 2 if v ∈ W n,1 (R d ), n > d, in addition to the initial assumptions (1.1b). The point is that α j → α ∞ = n-d 2d as j → +∞, α j → j 2 as n → +∞, α ∞ → +∞ as n → +∞;
(2.10) that is the convergence in (2.9) as E → +∞ is drastically better than the convergence in (2.8), at least, for large n and j.

Iterative step for phased inverse scattering

Recall that the outgoing fundamental solution to the Helmholtz equation is given by

G + (x, k) = -(2π) -d R d e iξx dξ ξ 2 -k 2 -i0 = i 4 k 2π|x -y| ν H (1) ν (k|x -y|) with ν := d 2 -1,
where H

(1) ν denotes the Hankel function of the first kind of order ν. Let G + (k) denote the convolution operator with kernel G + (•, k). The following estimate, which goes back to [START_REF] Agmon | Spectral properties of Schrödinger operators and scattering theory[END_REF], is essential for studies on direct scattering (see, e.g., [START_REF] Eskin | Lectures on linear partial differential equations[END_REF] ( §29), [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF], and references therein) and will also be crucial for our analysis:

Λ -s G + (k)Λ -s L 2 (R d )→L 2 (R d ) ≤ a 0 (d, s)|k| -1 , k ∈ R d , |k| ≥ 1, for s > 1 2 .
(2.11)

Here Λ -s denotes the operator of multiplication by (1 + |x| 2 ) -s/2 . Let v, satisfying (1.1b), be the unknown potential and v * E be an approximation to v, and assume that there exist constants A, E * , K, α > 0 such that

v * E ∈ L ∞ (R d ), supp v * E ⊂ D (2.12a) v * E -v L ∞ (D) ≤ AE -α (2.12b) v L ∞ (R d ) ≤ K, v * E L ∞ (R d ) ≤ K (2.12c)
for all E ≥ E * .

For inverse scattering with phase information the iterative step of [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF] is based on the following lemma.

Lemma 2.1. Let v, satisfying (1.1b), be the unknown potential and v * (•, E) be an approximation to v satisfying (2.12) for some A > 0, α ≥ 0, K > 0 and E * = E * (K, D), where

E * (K, D) = 2Ka 0 (d, s) sup x∈D (1 + |x| 2 ) s/2 ,
for some s > 1 2 , where a 0 (d, s) is the constant of (2.11).

(2.13)

Let f E , f * E be the scattering amplitudes of v, v * E . Then there exists a constant µ 1, §2.2 = µ 1, §2.2 (A, K, D) > 0 such that sup (k,l)∈M E |f E (k, l) -f * E (k, l) + v * E (k -l) -v(k -l)| ≤ µ 1, §2.2 E -α-1 2 , E ≥ E * , (2.14)
where v, v * E are the Fourier transforms of v, v * E defined according to (1.4). Note that in this paper we use the notation µ k, §X , k ≥ 1, for the constants of Section X. Lemma 2.1 follows from Lemma 3.2 of [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF] for v 0 ≡ 0, where v 0 is the background potential of [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF]. The proof of Lemma 3.2 of [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF] essentially uses estimate (2.11).

In particular, due to (2.14), the function

U * * E := Φf E -Φf * E + v * E , (2.15)
where Φf E , Φf * E are defined according to (2.5), satisfies the following improved error estimate compared to (2.12b):

U * * E -v L ∞ (B d 2 √ E ) ≤ µ 1, §2.2 E -α-1 2 , E ≥ E * . (2.16) If v ∈ W n,1 (R d )
, n > d (in addition to the initial assumptions (1.1b)), and if

α < n 2d -1 2 , (2.17)
then this permits to construct an improved approximation v * * E to the unknown potential v as follows:

v * * E (x) := B d r(E) e -ipx U * * E (p) dp, x ∈ D, 0, x ∈ D, r(E) := 2τ E 1+2α 2n , τ ∈ (0, 1], (2.18)
Here U * * E is defined in (2.15) and E * = E * (K, D) is the constant of (2.13). It follows that there exists a constant B > 0 such that

v * * E -v L ∞ (D) ≤ BE -β , E ≥ E * with β := α(1 -d n ) + 1 2 -d 2n , (2.19)
Note that α < β < n 2d -1 2 and that condition (2.17) implies that r(E) ≤ 2 √ E, so that the definition (2.18) is correct.

3 Iterative inversion from phaseless data

Low order potential reconstruction formulas from phaseless data

In this subsection we extend the formulas (1.3) and (2.8) to the phaseless case. The key result of [START_REF] Novikov | Explicit formulas and global uniqueness for phaseless inverse scattering in multidimensions[END_REF] consists in the following formulas for solving Problem 1 in dimension d ≥ 2 for L = 2 at high energies E:

sup p∈B d 2 √ E | v l (p)| 2 -|Φf E,l (p)| 2 = O(E -1 2 ), E → +∞, l = 0, 1, 2, (3.1)
where

v 0 = v, v l is defined by (1.7), l = 1, 2, f E,0 = f E , f E,1 , f E,2 are the scattering amplitudes for v 0 , v 1 , v 2 , respectively; in addition, Re v Im v = 1 2
Re

w 1 Im w 1 Re w 2 Im w 2 -1 | v 1 | 2 -| v| 2 -| w 1 | 2 | v 2 | 2 -| v| 2 -| w 2 | 2 , (3.2) where v = v(p), v l = v l (p), w l = w l (p), p ∈ R d ,
(p), Φ|f E | 2 (p), Φ|f E,1 | 2 (p), Φ|f E,2 | 2 (p) , |p| ≤ 2 √ E (3.4) Algorithm 1 function U = Urec(W 1 , W 2 , F, F 1 , F 2 )
data: Fourier transforms of reference potentials at some point p: W 1 = w 1 (p), W 2 = w 2 (p); scattering amplitude at (k, l), kl = p without reference potential:

F = |f (k, l)| 2
scattering amplitudes with reference potentials:

F 1 = |f 1 (k, l)| 2 , F 2 = |f 2 (k, l)| 2
result: approximation to the Fourier transform of the unknown potential v at p: U ≈ v(p)

M := Re W 1 Im W 1 Re W 2 Im W 2 ; b := F 1 -F -|W 1 | 2 F 2 -F -|W 2 | 2 ; x y := 1 2 M -1 b; U := x + iy
On the level of analysis (e.g., error estimates), the principal complication of (3.1), (3.2) in comparison with (1.3) consists in possible zeros of the determinant ζ w1, w2 of (3.3). This complication is, in particular, essential if one tries to transform (3.1), (3.2) into an approximate reconstruction in the configuration space, applying the inverse Fourier transform to v = Re v + i Im v of (3.2). For some simplest cases, the results of transforming (3.1), (3.2) to approximate reconstructions in the configuration space, including error estimates, were given in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF].

Background potentials of type A: The first simplest case analyzed in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] is

w 1 (x) := w(x -T 1 ), w 2 (x) := iw(x -T 1 ), x ∈ R d , (3.5)
for some w ∈ C(R d ) such that w = w, w(x) = 0 for |x| > R, and

w(p) = w(p) ≥ µ 1, §3.1 (1 + |p|) -σ , p ∈ R d , (3.6)
for some fixed

T 1 ∈ R d , R > 0, µ 1, §3.1 > 0, σ > d,
where T 1 and R are chosen in such a way that w 1 satisfies (1.5) (and, as a corollary, w 1 , w 2 satisfy (1.5) with Ω 2 = Ω 1 ). In addition, a broad class of w satisfying (3.6) was constructed in Lemma 1 of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]. One can see that (3.7)

ζ w1, w2 (p) = | w(p)| 2 ≥ µ 2 1, §3.1 (1 + |p|) -2σ , p ∈ R d , if w 1 ,
Background potentials of type B: The second simplest case analyzed in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] is

w 1 (x) = w(x -T 1 ), w 2 (x) = w(x -T 2 ), x ∈ R d , for some fixed T 1 , T 2 ∈ R d , (3.8)
where w is the same as in (3.6), and T 1 , T 2 , R are chosen in such a way that w 1 , w 2 satisfy (1.5).

One can see that where

ζ w1, w2 (p) = sin(py)| w(p)| 2 , y = T 2 -T 1 = 0, p ∈ R d , |ζ w1, w2 (p)| ≥ µ 2 1, §3.1 (1 + |p|) -2σ 2ε π , p ∈ R d \ Z ε y , if w 1 ,
Z ε y = p ∈ R d : py ∈ (-ε, ε) + πZ , y ∈ R d \ 0, 0 < ε < 1. (3.10)
First consider background potentials w 1 , w 2 of type A (see (3.5), (3.6)) and assume that v satisfies (1.1b), v ∈ W n,1 (R d ) for some n > d. Then the result of transforming U w1, w2 in (3.4) by

u E (x) := B d r(E) e -ipx U w1, w2 (p, E) dp, x ∈ R d , r(E) = 2τ E α n-d
for some fixed τ ∈ (0, 1], (3.11) to an approximate reconstruction in the configuration space is as follows (see [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]Theorem 1]):

u E -v L ∞ (R d ) = O(E -α ) E → +∞ with α = 1 2 (n -d) n + σ . (3.12)
Now consider background potentials w 1 , w 2 of type B (see (3.8)) and assume again that v satisfies (1.1b), v ∈ W n,1 (R d ) for some n > d. We transform U w1, w2 in (3.4) to an approximation u E in the configuration space as follows:

u E (x) = u E,1 (x) + u E,2 (x), x ∈ R d , u E,1 (x) = B d r(E) \Z ε(E) y e -ipx U w1, w2 (p, E) dp, u E,2 (x) = B d r(E) ∩Z ε(E) y e -ipx U ε(E) w1, w2 (p, E) dp, r(E) = 2τ E α n-d , ε(E) = E -α 2 for some τ ∈ (0, 1], (3.13)
Z ε y is defined in (3.10), and

U ε w1, w2 (p, E) = 1 2 U w1, w2 (p ε -, E) + U w1, w2 (p ε + , E) , (3.14)
p ε ± = p ⊥ + πz(p) y |y| 2 ± ε y |y| 2 , p ⊥ = p -(p • y) y |y| 2 , p ∈ B d 2 √ E ∩ Z ε y , for the unique z(p) ∈ Z such that |py -πz(p)| < ε. (3.15) Then it was shown in [1, Theorem 2] that u E -v L ∞ (R d ) = O(E -α ) E → +∞ with α = 1 2 (n -d) n + σ + n-d 2 , (3.16)
The geometry of vectors p, p ⊥ , y, p ε ± is illustrated in Fig. 1 for the case when the direction of y coincides with the basis vector e 1 = (1, 0, . . . , 0).

Approximate reconstruction of phased scattering data

We consider Problem 1 for d ≥ 2, L = 2, with the unknown potential v satisfying (1.1b) and with the background potentials w 1 , w 2 satisfying (1.5). Let

D ext = D ∪ Ω 1 ∪ Ω 2 , (3.17)
where D, Ω 1 , Ω 2 are the domains in (1.1b) and (1.5). Let v * E be an approximation to v satisfying (2.12) for some A > 0, α ≥ 0, K > 0 and for

E * = E * (K, D ext )
, where E * (K, D ext ) is defined according to (2.13). In addition, we suppose that Using the scattering amplitudes

v + w l L ∞ (R d ) ≤ K, v * (•, E) + w l L ∞ (R d ) ≤ K, E ≥ E * , l = 1, 2. (3.18) (x 2 , . . . , x d ) x 1 π |y| . . . πz(p) |y| p p ε - p ε + p ⊥ y O ε |y| ε |y|
f * E , f * E,1 , and f * E,2 of the known potentials v * E , v * E,1 := v * E + w 1 , and v * E,2 := v * E + w 2
, respectively, and the phaseless scattering data S of Problem 1, we construct an approximation f appr

E (p) to (Φf E )(p) for |p| ≤ 2 √
E by the function in Algorithm 1 as follows: E if ζ w1, w2 (p) = 0, where ζ w1, w2 is the determinant of formula (3.3). This follows from the estimate

f appr E (p) := Urec Φf * E,1 (p) -Φf * E (p), Φf * E,2 (p) -Φf * E (p), Φ|f E | 2 (p), Φ|f E,1 | 2 (p), Φ|f E,2 | 2 (p) (3.19) Note that f appr E (p) is well defined if ζ * (p, E) := det Re(Φf * E,1 (p) -Φf * E (p)) Im(Φf * E,1 (p) -Φf * E (p)) Re(Φf * E,2 (p) -Φf * E (p)) Im(Φf * E,2 (p) -Φf * E (p)) = 0. ( 3
|ζ * (p, E) -ζ w1, w2 (p)| ≤ µ 1, §3.2 E -1 2 , p ∈ B d 2 √ E , E ≥ E * , (3.21)
where µ 1, §3.2 > 0 and E * = E * (K, D ext ) is defined according to (2.13), (3.17). Estimate (3.21) follows from the definition of M in Algorithm 1, the formula

M -1 = 1 ζ Im W2 -Im W1 -Re W2 Re W1
, and from the estimates

v(p) -Φf E (p) ≤ µ 2, §3.2 E -1 2 , p ∈ B d 2 √ E , E ≥ E * , v l (p) -Φf E,l (p) ≤ µ 2, §3.2 E -1 2 , p ∈ B d 2 √ E , E ≥ E * , l = 1, 2, (3.22)
where

µ 2, §3.2 = µ 2, §3.2 (K, D ext ) > 0, E * = E * (K, D ext )
; see, e.g., [START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF].

Lemma 3.1. Let v, satisfying (1.1b), be the unknown potential of Problem 1 for d ≥ 2, L = 2.
Let w 1 , w 2 be the same as in (1.5). Let v * (•, E) be an approximation to v, satisfying (2.12) for some A > 0, α ≥ 0, K > 0 and for E * = E * (K, D ext ) defined according to (2.13), (3.17

). Let p ∈ B d 2 √ E be such that |ζ w1, w2 (p)| ≥ δ, (3.23)
for some fixed δ > 0. Then:

|ζ * (p, E)| ≥ δ 2 , E ≥ E * * , (3.24
)

Φf E (p) -f appr E (p) ≤ µ 3, §3.2 δ -1 E -α-1 2 , E ≥ E * * , (3.25)
E * * = max 4µ 2 1, §3.2 δ -2 , E * , (3.26)
where

f appr E is defined by (3.19), ζ * is defined by (3.20), µ 3, §3.2 is defined in (A.10) and µ 1, §3.2 is the constant of (3.21). Lemma 3.1 is proved in Section A.1.
The point is that the right-hand side of the estimate in (3.25) decays faster than the righthand side of the estimate in (2.12b) as E → +∞. This is a crucial advantage of f appr E as an approximation to the unknown phased scattering data Φf E in comparison with Φf * E . Using Lemma 3.1 we construct the iterative step for phaseless inverse scattering, see Sections 3.3 and 3.4.

Iterations for background potentials of type A

In this subsection we consider background potentials w 1 , w 2 satisfying (3.5) and (3.6).

Iterative step. We consider Problem 1 for d ≥ 2, L = 2, with the unknown potential v satisfying (1.1b) and with the background potentials w 1 , w 2 satisfying (1.5), (3.5), (3.6).

Let v * E be an approximation to v satisfying (2.12) for some A > 0, α ≥ 0, K > 0 and

E * = E * (K, D ext ), where E * (K, D ext ) is defined according to (2.13), (3.17) (with Ω 1 = Ω 2 ).
We construct an improved approximation v * * E to the unknown potential v via the scheme of Section 2.2 with Φf E replaced by f appr

E of formula (3.19) of Section 3.2. Put U * * E (p) = f appr E (p) -Φf * E (p) + v * E (p), p ∈ B d 2 √ E , ζ * (p, E) = 0, (3.27)
where f * E is the scattering amplitude of v * E , and ζ * is defined in (3.20). Under assumptions (3.5), (3.6), the iterative step for phaseless inverse scattering is realized as follows.

Theorem 3.2. Let v satisfy (1.1b) and v ∈ W n,1 (R d ) for some n > d. Let w 1 , w 2 be the same as in (1.5), (3.5), (3.6). Let v * E be an approximation to v satisfying (2.12) for some A > 0, α ≥ 0, K > 0 and for E * = E * (K, D ext ), where E * (K, D ext ) is defined according to (2.13), (3.17) (with

Ω 1 = Ω 2 ). We suppose also that α < n-d 2(d+2σ) , (3.28)
where σ is the constant of (3.6). Let

v * * E (x) := B d r(E) e -ipx U * * E (p) dp, x ∈ D, 0,
x ∈ D, where r(E) := 2τ E β n-d and β :=

( 1 2 +α)(n-d) n+2σ
(3.29) for some τ ∈ (0, 1]. Here U * * E is defined in (3.27), and B d r is defined by (2.3). Then there exist constants

B 1 = B 1 (τ, v n,1 , A, K, D ext , d, σ, n, µ 1, §3.1 ) and E 1 = E 1 (τ, α, A, K, D, d, σ, n, µ 1, §3.1 ) defined in (A.20) such that v * * E -v L ∞ (D) ≤ B 1 E -β for all E ≥ E 1 . (3.30) Remark 3.3. Under assumptions of Theorem 3.2, v * * E (x) is well-defined for E ≥ E 1 , i.e.: ζ * (p, E) = 0 for p ∈ B d r(E) , E ≥ E 1 , r(E) ≤ 2 √ E for E ≥ E 1 . (3.31)
Remark 3.4. The following three conditions are equivalent:

α < n-d 2(d+2σ) , β < n-d 2(d+2σ) , α < β, (3.32)
where β is defined in (3.29), n > d, α > 0, and σ > d.

Theorem 3.2 is proved in Section A.2. Theorem 3.2 of the present work can be considered as an extension of Theorem 1 of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] to the case when v * = 0. However, Theorem 3.2 of the present work for v * = 0 does not coincide with Theorem 1 of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF].

Iterations Let v, w 1 , w 2 satisfy the assumptions of Theorem 3.2. Let u

(1) E = v * * E for v * = 0. Note that u (1)
E is similar but does not coincide with the approximate reconstruction of Theorem 1 of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]; see formulas (3.12), (3.11) of the present work. In particular, we have that

u (1) E -v L ∞ (D) = O(E -α1 ), E → +∞ with α 1 := 1 2 n-d n+2σ .
(3.33)

Then, applying the iterative step described above in this subsection we construct nonlinear approximate reconstructions u

(j) E , j ≥ 2, such that u (j) E -v L ∞ (D) = O(E -αj ), E → +∞ with α j := 1 2 n-d d+2σ 1 -n-d n+2σ j . (3.34)
The approximations u

(j)
E of (3.34) for phaseless inverse scattering under assumptions (3.5), (3.6) are analogs of approximations u (j) E of (2.9) for phased inverse scattering. In a similar way with (2.10),

α j → α ∞ = 1 2 n-d d+2σ as j → +∞, α j → j 2 as n → +∞, α ∞ → +∞ as n → +∞, (3.35)
so that the convergence in (3.34) as E → +∞ is much more optimal then the convergence in (3.12), at least, for large j, n.

Iterations for background potentials of type B

In this subsection we consider Problem 1 for d ≥ 2 with shifted background potentials w 1 , w 2 described by (3.6), and (3.8) and unknown potential v satisfying (1.1b).

Iterative step. We consider the set Z ε y of formula (3.10). Put

p ε (p ⊥ , z, t) = p ⊥ + πz y |y| 2 + tε y |y| 2 , p ⊥ ∈ R d , p ⊥ • y = 0, z ∈ Z, t ∈ R.
(3.36)

Note that for any p ∈ Z ε y there exists the unique triple (p ⊥ , z, t) such that

p ε (p ⊥ , z, t) = p, p ⊥ ∈ R d , p ⊥ • y = 0, z ∈ Z, t ∈ (-1, 1).
(3.37)

In addition to U * * of (3.27), we also define

U * * N,ε (p ε (p ⊥ , z, t), E) = -N ≤j≤N,j =0 U * * p ε (p ⊥ , z, j , E)L j (t), (3.38)
under the assumptions that

ε < π N +1 , |p ε (p ⊥ , z, t)| ≤ 2 √ E -π |y| , t ∈ (-1, 1), ε > 0, N ≥ 1, (3.39)
where

L j (t) = (-1) N -j j(t + j) (N -j)!(N + j)! 1≤i≤N,i =j (t 2 -i 2 ), j = ±1, . . . , ±N. (3.40)
Note that for fixed p ⊥ ∈ R d , p ⊥ • y = 0, and for fixed z ∈ Z, function U * * N,ε (p ε (p ⊥ , z, t), E) is the Lagrange interpolating polynomial of degree 2N -1 in t ∈ R for U * * (p ε (p ⊥ , z, t), E) with the nodes at t = ±1, . . . , ±N . In addition, L j (t) is the j-th elementary Lagrange interpolating polynomial of degree 2N -1:

L j (t) = 1, t = j, 0, t = ±1, . . . , ±N, t = j. (3.41)
Note also that if assumptions (3.39) are valid for some

p ⊥ ∈ R d , p ⊥ • y = 0, z ∈ Z, t ∈ (-1, 1), then p ε (p ⊥ , z, s) ∈ B d 2 √ E \ Z ε y for all s ∈ [-N, -1] ∪ [1, N ]. (3.42)
Under assumptions (3.6), (3.8), the iterative step is realized as follows.

Theorem 3.5. Let v satisfy (1.1b) and v ∈ W n,1 (R d ) for some n > d. Let w 1 , w 2 be the same as in (1.5), (3.6), (3.8). Let v * E be an approximation to v satisfying (2.12a), (2.12) for some A > 0, α ≥ 0, K > 0 and for E * = E * (K, D ext ), where E * (K, D ext ) is defined according to (2.13), (3.17) (with Ω 1 = Ω 2 ). We suppose also that

α < n-d 2(d+2σ+ n-d 2N +1 ) for some N ≥ 1, (3.43)
where σ is the constant of (3.6). Let

v * * E (x) :=      B d r(E) \Z ε(E) y e -ipx U * * E (p) dp+, B d r(E) ∩Z ε(E) y e -ipx U * * N,ε(E) (p, E) dp x ∈ D 0, x / ∈ D with r(E) = 2τ E β n-d , ε(E) = E -β 2N +1 , β = ( 1 2 +α)(n-d) n+2σ+ n-d 2N +1
, for some τ ∈ (0, 1], Then the following estimate holds:

(
v * * E -v L ∞ (D) ≤ B 2 E -β , E ≥ E 2 , (3.45)
where

B 2 = B 2 (τ, y, N, v n,1 , A, K, D ext , d, σ, n, µ 1, §3.1
) and E 2 = E 2 (τ, y, α, N, A, K, D, d, σ, n, µ 1, §3.1 ) are defined in (A.49).

Theorem 3.5 is proved in Section A.3.

Remark 3.6. Under assumptions of Theorem 3.5, v * * E,1 (x) and v * * E,2 (x) are well-defined for E ≥ E 2 , i.e.:

ζ * (p, E) = 0 for p ∈ B d r(E) \ Z ε(E) y , E ≥ E 2 , r(E) ≤ 2 √ E for E ≥ E 2 , (3.46)
where ζ * is defined by (3.20), and also

ε(E) < π N +1 , |p| ≤ 2 √ E -π |y| for p ∈ B d r(E) , E ≥ E 2 .
(3.47)

Remark 3.7. The following three conditions are equivalent:

α < n-d 2(d+2σ+ n-d 4N +1 ) , β < n-d 2(d+2σ+ n-d 4N +1 ) , α < β, (3.48)
where β is defined in (3.44), n > d, α > 0, and σ > d. In addition, each of conditions (3.48) is equivalent to the following pair of conditions:

α < n-d 2(d+2σ) , N > 1 2 α -1 - n-d 2(d+2σ) -1 -1 -1 4 . (3.49)
Iterations. Let v, w 1 , w 2 satisfy the assumptions of Theorem 3.5. Let u

E = v * * E for v * = 0. Note that u (1) (1)
E is similar, but does not coincide with the approximate reconstruction of Theorem 2 of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]; see formulas (3.16), (3.13) of the present article. In particular, we have that

u (1) E -v L ∞ (D) + O(E -α1 ), E → +∞ with α 1 = 1 2 n-d n+2σ+ n-d 2N +1
.

(3.50)

Then, applying the iterative step described above in this subsection we construct nonlinear approximate reconstructions u (j) E , j ≥ 2, such that

u j E -v L ∞ (D) = O(E -αj ) E → +∞ with α j = 1 2 n-d d+2σ+ n-d 2N +1 1 - n-d n+2σ+ n-d 2N +1 j .
(3.51) These approximations u (j) for phaseless inverse scattering under assumptions (3.6), (3.8) are analogs of approximations u (j) of (2.9) for phased inverse scattering. In a similar way with (2.10) and (3.35),

α j → α ∞ = 1 2 n-d d+2σ+ n-d 2N +1
as j → +∞,

α j → j 2 as n → +∞, N → +∞, α ∞ → +∞ as n → +∞, (3.52)
so that the convergence in (3.51) as E → +∞ is much faster then the convergence in (3.16), at least, for large j, n, N .

4 Numerical experiments

Implementation of the Fourier transform and its inverse

The iterative algorithm presented in sections 3.3, 3.4 is implemented in Matlab in the twodimensional case. In our implementation we represent potentials v and w l by discrete functions v, w l defined on the space-variable grid

Γ s := 2 N (n 1 , n 2 ) : n 1 , n 2 ∈ Z N for N ∈ 2N, N ≥ 2 √ E π
where

Z N := -N 2 , -N 2 + 1, . . . , N 2 -1 . (4.1)
In turn, the input data |f (k, l)| 2 , |f l (k, l)| 2 are measured on a grid

M m E = {(k m , l m ) : m = 1..M } ⊂ M E ,
whose precise form depends on the experimental setup. This leads to the following grid in Fourier space:

Γ m := {k m -l m : m = 1..M }
Minimal data. To approximate the inverse Fourier transform of a function, which is supported on B 2

2 √ E and sampled on the grid Γ m , by the Fast Fourier transform (FFT), Γ m has to be rectangular. If Γ s is given by (4.1), a minimal choice of Γ m is

Γ m min := B 2 2 √ E ∩ πZ 2 . ( 4.2)
A corresponding measurement grid M m E,min can be defined as in (2.4) replacing B 2

2 √ E by Γ m min . Extending a function given on Γ m to the exterior grid

Γ e := {π(n 1 , n 2 ) : n 1 , n 2 ∈ Z N , π 2 (n 2 1 + n 2 2 ) > 4E} (4.3)
by 0, we can compute an approximation to the inverse Fourier transform on Γ s by FFT.

Discrete Ewald circles. The above choice Γ m min of the set of measurement points is inconvenient both from an experimental and from a computational point of view since each point in Γ m min corresponds to a different incident wave. If a scattering experiment is performed for some incident wave or if the solution to the scattering problem is computed numerically, the resulting far field pattern can be evaluated at other points without essential additional costs. Therefore, we now consider input data for M 1 uniformly distributed incident wave vectors k where each far field pattern is evaluated at M 2 uniformly distributed scattered wave vectors l:

k(s) := √ E cos(2π s M1 ), sin(2π s M1 ) , l(s, t) := √ E cos( 2πs M1 + 2πt M2 ), sin( 2πs M1 + 2πt M2 ) , s ∈ Z M1 , t ∈ Z M2 , (4.4)
resulting in Here it is necessary to include the points in Γ e to obtain small condition numbers of A since the inverse Fourier transform is computed by numerically inverting A. Matrix-vector products with A and A * can be computed efficiently without the need to set up and store the matrix A using the Nonequispaced Fast Fourier Transform (NFFT). In our work we use the NFFT implementation of [START_REF] Keiner | Using NFFT 3-a software library for various nonequispaced fast fourier transforms[END_REF]. The definitions of the grids Γ s , Γ m , Γ e and of the Fourier transform matrix A are summarized in Algorithm 2.

M m E,M1,M2 := (k(s), l(s, t)) | s ∈ Z M1 , t ∈ Z M2 , M = M 1 M 2 , (4.5)
Γ m M1,M2 := {k -l : (k, l) ∈ M m E,M1,M2 }. ( 4
Approximation of the inverse Fourier transform. A first idea may be to approximate the continuous inverse Fourier transform by the Moore-Penrose inverse of A. If A is injective, this corresponds to the solution of a least-squares problem

A † U = argmin w Aw -U 2 2
which can be achieved by the conjugate gradient (CG) method applied to the normal equation (A * A)w = A * U . The ideal situation would be that A is isometric, i.e. Aw 2 = w 2 for all w : Γ s → C. In this case A † = A * , and the CG method would yield the exact solution in the

Algorithm 2 function [Γ m , Γ e , A, Φ] = fourier setup(Γ s , M m E ) data: spatial grid: Γ s := { 2 N (n 1 , n 2 ) : n 1 , n 2 ∈ Z, -N 2 ≤ n 1 , n 2 < N 2 } with N ∈ 2N, N ≥ 2 √ E π measurement points: M m E = {(k m , l m ) : m = 1..M } ⊂
m := {k m -l m : m = 1..M }; grid inside B 2 2 √ E Γ e := {π(n 1 , n 2 ) : n 1 , n 2 ∈ Z, -N 2 ≤ n 1 , n 2 < N 2 , π 2 (n 2 1 + n 2 2 ) > 4E}; grid outside B 2 2 √ E A := (πN ) -2 exp(ix • p) p∈Γ m ∪Γ e ,x∈Γ s ; avoid computation of A and use NFFT instead (Φf ) p := 1 #(Mp)
m∈Mp f m where M p := {m : k m -l m = p} first step. However, in our situation A is typically far from being isometric, so that A * A is far from the identity matrix, and the CG method requires a big number of iterations. The reason is that even though the continuous Fourier transform is isometric, the Euclidean norm U 2 of the sampled version

U : Γ m ∩ Γ e → C of a function U : [-N π 2 , N π 2 ] 2 → C can be a bad approximation for U L 2 .
To overcome this difficulty, we design a weight matrix

D such that D 1/2 U | 2 ≈ U L 2 .
Then we approximate the inverse Fourier transform by the Moore-Penrose inverse of A with respect to this weighted norm D 1/2 • 2 , or in matrix notation

(D 1/2 A) † D 1/2 U = argmin w D 1/2 Aw -D 1/2 U 2 2 .
Recall that by the first-order optimality conditions, which are necessary and sufficient for convex functionals, this minimization problem is equivalent to solving the normal equation

Implementation of the inversion method

Phaseless Born approximation. Recall that in our implementation the potentials w l are represented by discrete functions w l : Γ s → C, and the measured phaseless farfield data |f | 2 , |f l | 2 are represented by the discrete functions F and F l : M m E → R. Also note that in addition to the background potentials and phaseless farfield data, a cutoff radius r 1 > 0 and a threshold δ > 0 are specified as input data for the algorithm. The cutoff radius 0 < r 1 ≤ 2 √ E is analogous to the radius r(E) of formulas (3.11), (3.13), whereas the threshold δ is analogous to threshold ε(E) of formula (3.13).

The implementation of the phaseless Born approximation is shown in Algorithm 3. The algorithm is formulated for an arbitrary number L ≥ 2 of reference potentials, but reduces to the algorithm in the theoretical part of this paper if L = 2. The principal part is the computation of the reduced grid Γm ⊂ Γ m , which consists of points p ∈ Γ m meeting the threshold and the cutoff constraints, and of the function U : Γm ∪ Γ e → C, which is the discrete version of the function U w1, w2 of Subsection 3.1. This computation starts by initializing Γm by the empty grid and U by the zero function and proceeds as follows:

1. For each point p ∈ Γ m and each pair w l, w l of reference potentials, compute the determinant ζ l,l p corresponding to ζ w1, w2 (p) of formula (3.3). Since the approximate Fourier transform U p of v can be computed from any pair ( l, l) for which ζ l,l p = 0 and since the computation is the more stable the larger |ζ l,l p |, we choose the pair ( l(p), l(p)), for which |ζ l,l p | is largest.

If |ζ l(p),l(p) p

| > δ (threshold constraint) include p in Γm . In addition, if |p| < r 1 (cutoff constraint) compute U p using Algorithm 1 with appropriate parameters. Rather than implementing an explicit interpolation scheme at points where ζ l(p),l(p) p vanishes or is too small as done in (3.14) for our theoretical analysis, we use a trigonometric interpolation of the discrete Fourier transform induced by fitting to the remaining points of Γ m . This procedure is easier to implement and allows the numerical treatment of arbitrary background potentials.

The last step is to compute the function v * E : Γ s → C, which is the discrete version of the phaseless Born approximation u E of Subsection 3.1, as the inverse Fourier transform of the function represented by U . This step is explained in Subsection 4.1.

Algorithm 3 function [v * E , Γ m , Ã, D] = phaseless born inv (F 0 , F 1 , w 1 , . . . , F L , w L , Γ s , M m E , r 1 , δ) data: measured data (F 0 ) m ≈ |f (k m , l m )| 2 , (F l ) m ≈ |f l (k m , l m )| 2 for m = 1..M , l = 1..L (discrete) background potentials w 1 , . . . , w L : Γ s → C Γ s , M m E as
D := voronoi weights Γm ∪ Γ e ; v * E := (D 1/2 Ã) † D 1/2 [U p ] p∈ Γm ∪Γ e ;
may be computed by CG and NFFT Iterative algorithm. In addition to the input parameters of the phaseless Born approximation, the iterative algorithm requires cutoff radii 0 < r 1 ≤ • • • ≤ r J ≤ 2 √ E to be specified. The implementation of the iterative method is shown in Algorithm 4.

The first step is the computation of the phaseless Born approximation v * E : Γ s → C, as well as of the grids Γ m , Γ e , Γm , Fourier matrices A, Ã and Voronoi's matrix of weights D, as explained above. The main part of the algorithm is the iteration procedure producing improved approximations v * E . The iterative step starts by evaluating the scattering amplitudes f * E and f * E,l , l = 1, ..L of the potentials represented by discrete functions v * E , v * E + w l at the points of the grid

M m E yielding discrete functions f * E , f * E,l : M m E → C.
In principle, any black-box solver can be used to evaluate the scattering amplitudes. In our work we use the solver described in [START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF][START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF].

The iterative step proceeds by computing the discrete function f appr and U * * by the zero functions and continues as follows:

E : Γm ∪ Γ e → C,
1. For each point p ∈ Γm such that |p| < r j , where j ≥ 2 is the current iteration number4 , compute ( f appr

E

) p using Algorithm 1 with appropriate parameters. If L ≥ 3 reference potentials are available, we again choose the most stable pair at each point p.

Evaluate U * *

p according to formula (3.27).

The iteration ends by computing v * E , as the inverse Fourier transform of the function represented by U * * . The computation is explained in Subsection 4.1.

Algorithm 4 function v * E = phaseless iterative inv(F 0 , F 1 , w 1 ,..,F L , Γ s , M m E ,w L , r 1 ..r J , δ) data: measured data (F 0 ) m ≈ |f (k m , l m )| 2 , (F l ) m ≈ |f l (k m , l m )| 2 for m = 1..M , l = 1..L (discrete) background potentials w 1 , . . . , w L : Γ s → C Γ s , M m E as in fourier setup cutoff radii 0 < r 1 ≤ r 2 ≤ • • • ≤ r J ≤ 2 √ E threshold δ > 0 result: potential reconstruction v * E [Γ m , Γ e , A, Φ] := fourier setup(Γ s , M m E ); see Algorithm 2 [v * E , Γ m , Ã, D] := phaseless born inv (F , F 1 , w 1 , . . . , w L , F L , Γ s , M m E ,
E ) p := Urec (Φf * E, l(p) -Φf * E ) p , (Φf * E,l(p) -Φf * E ) p , (ΦF 0 ) p , (ΦF l(p) ) p , (ΦF l(p) ) p ; U * * p := (Av * E ) p + ( f appr E ) p -(Φf * E ) p ; end for v * E := (D 1/2 Ã) † D 1/2 U * * ; appr
may be computed by CG and NFFT end for

Numerical results

Reconstruction errors. We consider the reconstruction of the potential v shown in Fig. 3 (a)5 using three background potentials w 1 , w 2 , w 3 shown at Fig. 3 (b),(c),(d). The differential scattering cross-section of the potential v for the incident direction k = √ E(-1, 0) at E = 10 2 and E = 20 2 is shown at Fig. 4. One can see that for bigger energy the differential scattering cross-section is more concentrated near l ≈ k.

In the experiments 32 equidistant incident directions and 256 equidistant measurement directions (for each k) were used. Moreover we choose N ≈ c √ E with c = 5 so that the space grid discretization step is 2/N ≈ 0.4/ √ E. For simplicity, we choose uniformly increasing cutoff radii

r j = √ E 1 + j J+1 .
-1 0 1 In quantum mechanical and optical applications the dominant source of noise is often caused by the limited number N p of measured particles, leading to Poisson distributed data. More precisely, recall that the quantity |f (k, l)| 2 is proportional to the probability density of scattering of a particle with initial momentum k into direction l/|l| = k/|k|. Let K := #{k 1 , . . . , k m } denote the number of incident waves. We assume that for each background potential and each incident wave the exposure time t l (k m ) is chosen such that the same expected number of particles N p /K(L + 1) is recorded. Thus, our simulated noisy data F := [F 0 , F 1 , . . . , F L ], were generated from exact data

F † := [F † 0 , F † 1 , . . . , F † L ] with (F † 0 ) m := |f (k m , l m )| 2 and (F † l ) m := |f l (k m , l m )| 2 , l = 1, . . . , L by (F l ) m ∼ 1 t l (k m ) Pois t l (k m )(F † l ) m with t l (k m ) := N p K(L + 1)σ l (k m ) , σ l (k m ) := n : kn=km (F † l ) n
(see [START_REF] Hohage | Inverse problems with Poisson data: statistical regularization theory, applications and algorithms[END_REF] for more details). Here Pois(x) stands for a Poisson random variable with mean x. Recall that E(F l ) m = (F † l ) m and that the pointwise noiselevel is the potential in Fig. 3 with E = 15 2 the average pointwise noise level F -F † 2 / F † 2 was about 44%, 14%, 4.4% and 1.4% for total count numbers N p ∈ {10 7 , 10 8 , 10 9 , 10 10 }. However, we stress that pointwise noise levels, although frequently used, are misleading as they tend to infinity as the discretization of the data space becomes finer and finer without loss of information, and that N p (or

Var(F l ) m = t l (k m ) -1/2 (F † l ) m . For
N -1/2 p
) is a better characterization of the noise level. As the proposed method only yields an approximate solution at fixed energy for noiseless data, a natural question is how much the reconstructions of our method can be improved by iterative regularization methods. Here we choose the Newton conjugate gradient method (NewtonCG) (see [START_REF] Hanke | Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems[END_REF]) with H 1 inner product in the preimage space as a commonly used representative of this class of methods, and use the results of our method as starting point of NewtonCG. The stopping index of the NewtonCG iteration was chosen by the discrepancy principle using the estimated noise level E F -

F † 2 2 = l,m Var(F l ) m = l,m t l (k m ) -1 (F † l ) m ≈ t l (k m ) -1 (F l ) m (
even though the discrepancy principle is actually only justified for deterministic noise models, see [START_REF] Hanke | Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems[END_REF]).

The number of iterations J is chosen basing on the following observations: for small particle count such as N p = 10 7 after two or three iterations the accumulated noise in reconstruction by our method is already comparable with the reconstruction error and the method should be stopped; for bigger particle counts such as N p ≥ 10 10 the number of steps can be chosen in a relatively large range without significant impact on the reconstruction error; we use J = 8 iterations for large particle counts N p ≥ 10 8 .

Figure 5 shows cross-sections of the reconstructed potentials for the reconstructions using (a) the phaseless Born approximation, (b) our method and (c) our method combined with NewtonCG. In this experiment we used the values of parameters E = 10 2 and N p = 10 9 and the highest possible scaling factor for the potential shown in Fig. 3 (a), for which our method still works at this energy level. This corresponds to (L 1 , L 2 , L ∞ ) norms (13.5, 14.3, 37.0), respectively.

Table 1 shows the reconstruction errors using the same methods for different energy levels and different expected count numbers N p . Errors are averaged over 5 experiments. These results demonstrate that our method performs well far beyond the scope of validity of the Born approximation. It can also be seen from these tables that the best reconstruction results are achieved by combining the proposed method with an iterative regularization method such as NewtonCG.

Comparison of convergence regions with NewtonCG. We also studied the influence of scaling of background potentials on reconstructions using our method and using NewtonCG. We made tests for the unknown potential of Fig. 3 1: relative L ∞ reconstruction errors in percents for the potential v and the background potentials w 1 , w 2 , w 3 shown in Fig. 3.

First, we tried to recover the potential using NewtonCG with zero initial guess. Simulations show that the iterations do not converge to the exact potentials unless the potential is downscaled by a factor of 75 or bigger. On the other hand, our method works reasonably well, see Tab. 2, column 3B(W), and it can provide an initial guess from which NewtonCG converges.

We also noticed that the simultaneous downscaling of the unknown and background potentials, even by a factor of 1000, does not solve the convergence problem of NewtonCG: the norm ratio of the background potential to the unknown potential must be also sufficiently small to guarantee the convergence of NewtonCG from the zero initial guess.

Non-smooth potentials. Consider the potential of 6 (a). The differential scattering crosssections of this potential at energies E = 2 2 and E = 10 2 for the incident direction k = √ E(0, 1) is shown at Fig. 4 (right).

Fig. 6 (b), (c) shows reconstructions of this non-smooth potential using the Born approximation and our method without NewtonCG. In the experiment we use background potentials w 1 , w 2 of type A (i.e. w 2 = iw 1 ), the energy is E = 10 2 , and the particle count is N p = 10 15 . One can see that even though our method is theoretically justified only for sufficiently smooth potentials, it performs well also for non-smooth potentials.

Choice of background potentials. Examples of functions with compact support satisfying assumption (3.6) are Wendland's radial basis functions (see [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF]), in particular

w(x) = Φ k (|x| 2 /h)
for some h > 0, k ∈ {0, 1} and (4. More precisely, (3.6) is satisfied with σ = d + 2k + 1 (see [START_REF] Wendland | Error estimates for interpolation by compactly supported radial basis functions of minimal degree[END_REF]), and w ∈ C 2k (R d ).

We also checked to which extent our approach still works if we use background potentials which do not satisfy assumption (3.6), but may be easier to realize experimentally, such as indicator functions of squares. Our results are documented in Table 2 and Fig. 7. It can be seen that non-smooth background potentials yield better results than smooth background potential, but our method still works reasonably well for the C 2 -Wendland function in (4.7). Moreover, the best results are obtained for the nonsmooth rectangular potentials even though they do not satisfy assumption (3.6). Furthermore, our method yields good results for two background potentials of type A, but significantly worse results for two background potentials of type B. However, the results for type B background potentials can be improved to a quality comparable to type A potentials if either a subsequent NewtonCG iteration is used or if data from a third shifted potential are available.

Robustness against errors in background potentials. In practice it is usually not possible to measure the background potentials exactly. Therefore, we tested our algorithm in the case where the simulated data are generated using a potential w which is a perturbation of the potential w used in our reconstruction method. Figure 8 shows cross-sections of the potentials w and w, and a cross-section of the reconstruction of v. In this example the potential w is obtained from w by amplitude scaling (×1.3), support scaling (×0.8), translation by (0.1, 0), addition of Gaussian noise with standard deviation .22 and convolution with Gaussian kernel with standard deviation 0.5. We use the same unknown potential v of Figure 3 as before and set E = 15 2 , N p = 10 9 .

In this example the phaseless Born approximation, our method and our method combined with NewtonCG give the relative L ∞ errors (64%, 3.9%, 1.5%), respectively. This demonstrates a remarkable robustness of our method against errors in the reference potentials.

Conclusions

We have proposed a method for the solution of phaseless inverse medium scattering problems in the presence of known background potentials. Let us summarize the advantages and disadvantages of our method in comparision with iterative regularzation methods such as NewtonCG:

• Global convergence: Iterative regularization methods require a good initial approximation to the unknown potential, whereas we have shown global convergence of our method as the energy tends to infinity. Numerical experiments at fixed energy demonstrate excellent performance of our method for large potentials and/or weak background potentials where NewtonCG failed.

• Computation time. Each iteration step of our reconstruction method requires only (L + 1) solutions of a forward problem, where L is the number of background potentials, and other comparatively cheap operations. In contrast, regularized Newton methods additionally require the solution of a linearized inverse problem in each iteration step, and they typically need a larger number of iteration to achieve the accuracy of our method. In our experiments the proposed method was typically more than 20 times faster than NewtonCG, but we stress that the quotient of computation times strongly depends on the noise level, the energy, the potentials, the choice of the direct solver (setup time vs. solution time), and other parameters.

• Use of black-box solvers. As our method only requires the solution of forward scattering problems, any block-box solver for such problems can be used. In contrast, iterative regularization methods additionally use the Fréchet derivative of the forward operator and typically also the adjoint of the Fréchet derivative. The implementation of these operations may require modifications of the source code of the forward solver.

• Asymptotic exactness. At fixed energy in the absence of noise the NewtonCG is expected to converge to the exact solution under some additional assumptions, in particular the uniqueness of solution and the tangential cone condition, see [START_REF] Hanke | Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems[END_REF]. In contrast, theoretically our method will converge to the exact potential only in the limit E → +∞.

• Stopping rules. There exists a considerable literature on a-posteriori stopping rules for iterative regularization methods and their convergence properties such as the discrepancy principle for NewtonCG (see [START_REF] Hanke | Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems[END_REF]). In contrast, we have used a rather ad-hoc a-priori stopping rule in our experiments with the proposed method for the lack of a better alternative.

This discussion shows that the pros and cons of our method are rather complementary to those of iterative regularization methods. Therefore, the proposed method provides a valuable new tool for the solution of phaseless inverse medium scattering problems. Hybrid methods using the proposed method to compute an initial guess for an iterative regularization method allow to combine the advantages of both methods. But in many cases the proposed method itself may already provide a sufficiently accurate reconstruction.

On the theoretical side we have demonstrated fast convergence of our method as energy tends to infinity for exact data and two types of background potentials. It remains for future research to study the behavior of the proposed method in the presence of noise and to devise and analyze useful stopping rules.

A Proofs

A.1 Proof of Lemma 3.1

Under the assumptions of Lemma 3.1, the following estimates hold:

|Φf E (p) + Φf * E,j (p) -Φf E (p)| -|Φf E,j (p)| ≤ µ 1, §A.1 E -α-1 2 , p ∈ B d 2 √ E , E ≥ E * , j = 1, 2, (A.1)
where µ 1, §A.1 > 0. The proof of estimates (A.1) is based on the following inequalities: E , E ≥ E * . Formula (A.1) follows from (A.2), (A.3).

|Φf E + Φf * E,l -Φf * E | -|Φf E,j | ≤ |Φf E + Φf * E,l -Φf * E -Φf E,l | = (Φf E -v) -(Φf * E -v * E ) -(Φf E,l -v E,l ) -(Φf * E,l -v * E,l ) ≤ (Φf E -v) -(Φf * E -v * E ) + (Φf E,l -v l ) -(Φf * E,l -v * E,l ) (A.
Next, using (A. Using the formula for the inverse matrix, we also get the following equality: (A.40)

M -1 (p, E) F = |Φf * E,1 -Φf * E | 2 + |Φf * E,2 -Φf * E | 2

-1 2 f 2 ,

 22 s (x, k) → 0 as |x| → ∞ (1.1c)uniformly in x/|x|. This implies that ψ s has the asymptotic behaviorψ s (x, k) = c(d, |k|) e i|k||x| |x| dE k, |k| x |x| + O |x| -d+1 |x| → ∞, c(d, |k|) := -πi(-2πi) (d-1)/2 |k| (d-3)/2

  w 2 are defined by (3.5),(3.6).

  w 2 are defined by (3.6), (3.8),(3.9)

Figure 1 :

 1 Figure 1: Vectors p, p ⊥ , y and p ε ± of formula (3.15)

. 20 )

 20 Note that condition(3.20) is satisfied for sufficiently large E at fixed p ∈ B d 2 √

. 6 )M 1 = 6 , 1 = 8 , M 2 = 30 Figure 2 :

 61618302 Figure 2: Illustration of sets of transformed measurement points Γ m M1,M2 on Discrete Ewald circles for different values of the numbers M 1 of incident waves and far field points M 2 . The right panel shows the Voronoi diagram corresponding to the grid Γ e ∪ Γ m M1,M2 . The color is related to the area of a cell. Given a discrete function v : Γ s → C representing a function v : R 2 → C, the Fourier transform of v can be approximately represented by a discrete function v : Γ m ∪ Γ e → C such that v = Av, where A := (N π) -2 exp(ix • p) p∈Γ m ∪Γ e ,x∈Γ s .

  A * DAw = A * DU . To construct a matrix D such that D 1/2 U | 2 2 ≈ |U (p)| 2 dp, we use a Voronoi partition of the square [-πN 2 , πN 2 ] 2 into cells C(p) centered at points p ∈ Γ m ∪ Γ e . In Matlab this subdivision is computed by the built-in function voronoi, see Fig. 2. To approximate the integral by a Riemann sum we evaluate the area |C(p)| of each cell C(p), p ∈ Γ m ∪ Γ e and choose D as the diagonal matrix D := diag(|C(p)|) p∈Γ m ∪Γ e . The use of this matrix D drastically decreases the number of CG steps and allows the approximate evaluation of the inverse Fourier transform with just a few CG steps. The reason is that A * DA is much closer to the identity matrix as A * A.

Remark 4 . 1 .

 41 If the function U is defined only on some subgrid G ⊂ Γ m ∪ Γ e , two minor modifications are necessary: (i) The Fourier matrix A must be restricted to the grid G yielding matrix à = [A pq ] p∈Γ,q∈Γ s ; (ii) The matrix of Voronoi weights D must be computed for the grid G.

  in fourier setup cutoff radius r 1 > 0 threshold δ > 0 results: potential reconstruction: v * E reduced Fourier space grid and Fourier transform matrix: Γ m , A Voronoi diagonal weight matrix: D [Γ m , Γ e , A, Φ] = fourier setup(Γ s , M m E ); see Algorithm 2 w l := A w l for l = 1, . . . , L; use NFFT P := {( l, l) ∈ {1, .., L} 2 : l < l}; Γm := [ ]; U p := 0 for p ∈ Γ m ∪ Γ e ; for p ∈ Γ m do for ( l, l) ∈ P do ζ l,l p = Re w l,p Im w l,p -Im w l,p Re w l,p ; end for Choose ( l(p), l(p)) ∈ argmax ( l,l)∈P |ζ l,l p |; if |ζ l(p),l(p) p | > δ then Γm := Γm ∪ {p}; if |p| < r 1 then U p := Urec w l(p),p , w l(p),p , (ΦF 0 ) p , (ΦF l(p) ) p , (ΦF l(p) ) p ; see Algorithm 1 end if end if end for à := [A px ] p∈ Γm ∪Γ e ,x∈Γ s ;

  r 1 , δ); see Alg. 3 P := {( l, l) ∈ {1, .., L} 2 : l < l}; for j = 2, . . . , J do f * E := scattering amplitude(v * E ); Solve phased forward problem f * E,l := scattering amplitude(v * E + w l ) for l = 1..L; f appr E := 0; U * * p := 0 for p ∈ Γm ∪ Γ e ; grid functions on Γ m ∪ Γ e for p ∈ Γ m such that |p| < r j do Choose ( l(p), l(p)) ∈ argmax ( l,l)∈P Re w l,p Im w l,p -Im w l,p Re w l,p ; ( f

  v + w 3

Figure 3 :

 3 Figure 3: Unknown potential and background potentials

4 Figure 4 :

 44 Figure 4: Differential scattering cross-sections |f (k, l)| 2 for the potential of Fig. 3 (a) at k = √ E(-1, 0) (left) and for the potential of Fig. 6 (a) at k = √ E(0, 1) (right). The distance to the curve in direction l is equal to |f (k, l)| 2 .

  (a) and three type-B potentials based on the Wendland functions with k = 1.

Figure 5 :

 5 Figure 5: Cross-sections of the exact v and reconstructed potentials v * , and L ∞ relative reconstruction errors. Here E = 10 2 and N p = 10 9 .

7 )Φ 0

 70 (r) := max(1 -r, 0) 2 , Φ 1 (r) := max(1 -r, 0) 4 (4r + 1).(4.8)

Figure 6 :Table 2 :

 62 Figure 6: Reconstruction of a non-smooth potential for rectangular background potentials of type A, E = 10 2 , and N p = 10 15 . 2A(R) 2A(W) 3B(R) 3B(W) 2B(R) 2B(W) Born approximation 20 20 19 21 20 25 our method 2.2 3.4 1.2 5.8 6.5 16 our method+NewtonCG 0.95 1.9 0.95 2.3 2 4.4

Figure 7 :Figure 8 :

 78 Figure 7: Background potentials based on the Wendland functions Φ 0 , Φ 1 defined in (4.7) (left), cross-sections of the unknown potential v and its reconstructions v * using our method for three type-B background potentials based on Φ 0 (center) and for three type-B background potentials based on Φ 1 (right)

2 ) for l = 1 , 2 . 2 ,µ 1 ,

 21221 Applying Lemma 2.1 to v, v * E and to v l , v * E,l , respectively, we get the estimates(Φf E (p) -v(p)) -(Φf * E (p) -v * E (p)) ≤ 1 2 µ 1, §A.1 E -α-1 (Φf E,l (p) -v l (p)) -(Φf * E,l (p) -v * E,l (p)) ≤ 1 2 µ 1, §A.1 E -α-§A.1> 0 is the constant in the right-hand side of formula (3.14) of[START_REF] Novikov | An iterative approach to non-overdetermined inverse scattering at fixed energy[END_REF], p ∈ B d 2 √

1 2 |ζ 2 √E 2 √ 2 3 δ -1 µ 2 , 2 Proposition A. 1 . 1 2 √ 1 ≤ r ≤ 2 √ 1 ≤ r ≤ 2 √

 12222221121212 * (p, E)| , E ≥ E * * , (A.8) for p ∈ B d such that (3.23) holds and where • F denotes the Frobenius matrix norm: |a ij | 2 for any square matrix A = (a ij ) n i,j=1 . Using (3.24), (A.5), (A.8) we obtain the estimate M -1 (p, E) F ≤ §A.1 , E ≥ E * * . (A.9) Finally, using (A.7) and (A.9), we get the estimate (3.25), where µ 3, §3.2 = 2 3 µ 1, §A.1 µ 2 2, §A.1 . (A.10) Lemma 3.1 is proved.A.2 Proof of Theorem 3.Let v and w 1 , w 2 be the same as in (1.1b), (3.6),(3.8). Let v * E be an approximation to v satisfying (2.12). Let U * * E (p) be defined according to(3.27), Then|U * * E (p) -v(p)| ≤ µ 1, §A.2 E -α-1 2 r 2σ , µ 1, §A.2 := 4 σ µ -2 1, §3.1 µ 3, §3.2 + µ 1, §2.2 , E ≥ max{µ 2, §A.2 r 4σ , E * }, µ 2, §A.2 := 4 2σ+1 µ -4where σ is the same as in (3.6), E * = E * (K, D ext ) is defined according to (2.13), (3.17), and µ k, §X , k ≥ 1, are the constants of Section X.Proof of Proposition A.1. Due to (2.13), (3.27) and Lemma 3.1, we have|U * * E (p) -v(p)| ≤ (µ 3, §3.2 δ -1 + µ 1, §2.2 )E -α-1 2 , if p ∈ B d 2 √ E , |ζ w1, w2 (p)| ≥ δ, E ≥ E * * , E * * = max 4µ 2 1, §3.2 δ -2 , E * . (A.12)Besides, in view of (3.9), we also have that minp∈B d r |ζ w1, w2 (p)| = µ 2 1, §3.1 (1 + r) -2σ ≥ µ 2 1, §3.1 4 -σ r -2σ , 1 ≤ r. (A.13)Using (A.12) with δ = µ 2 1, §3.1 4 -σ r -2σ and (A.13), we get (A.11).Proof of Proposition A.4. As in the proof of Proposition A.1, we have formula (A.12). Besides, it follows from (3.7) that|ζ w1, w2 (p)| ≥ 4 -σ 2 π µ 2 1, §3.1 r -2σ ε, p ∈ B d r \ Z ε y , 0 < ε < 1, 1 ≤ r ≤ 2 √ E. (A.34)Using (A.34) and (A.12) with δ = 4 -σ 2 π (µ 1, §3.1 ) 2 r -2σ ε, we get (A.30). It remains to prove (A.32). Using definition (3.38), one can writeU * * N,ε (p, E) -v(p) = ϕ N,ε (p, E) + V N,ε (p), for p = p ε (p ⊥ , z, t) ∈ B d r ∩ Z ε y , r ≤ Eπ |y| with ϕ N,ε (p, E) := -N ≤j≤N,j =0 U * * p ε (p ⊥ , z, j), Ev p ε (p ⊥ , z, j) L j (t), (A.35)where V N,ε is given by (A.22).Using estimate (A.30) and formulas (3.38), (3.42), (A.35) we get|ϕ N,ε (p, E)| ≤ µ 2, §A.3 E -α-1 2 r + π |y| 2σ ε -1 -N ≤j≤N,j =0 |L j (t)| ≤ 4 σ max(1, π |y| ) 2σ µ 2, §A.3 E -α-1 2 ε -1 r 2σ -N ≤j≤N,j =0 |L j (t)|, p = p ε (p ⊥ , z, t) ∈ B d r ∩ Z ε y , t ∈ (-1, 1) and j = ±1, . . . , ±N ,|L j (t)| = |j||t + j| (Nj)!(N + j)! 1≤i≤N,i =j (i 2 -t 2 ) ,j =0 |L j (t)| ≤ 4N 1 -N N +1 N , t ∈ (-1, 1). (A.38) Using (A.36), (A.38), we get|ϕ N,ε (p, E)| ≤ µ 4, §A.3 E -α-1 2 r 2σ ε -1 , for p = p ε (p ⊥ , z, t) ∈ B d r ∩ Z ε y , Eπ |y| , E ≥ max µ 3, §A.3 r 4σ ε -2 , E * .(A.39) Using (A.23), (A.35), (A.39), we get (A.32).Proposition A.4 is proved.The final part of the proof of Theorem 3.5 is as follows. In a similar way with (A.14), we represent v as follows: v(x) = v + 1 (x, r) + v + 2 (x, r) + v -(x, r), x ∈ D, r > 0 with v + 1 (x, r) := B d r \Z ε y e -ipx v(p) dp, v + 2 (x, r) := B d r ∩Z ε y e -ipx v(p) dp, v -(x, r) := R d \B d r e -ipx v(p) dp.

  9 .

	N p \E 10 2 15 2 20 2	N p \E 10 2 15 2 20 2	N p \E 10 2 15 2	20 2
	10 7 53	20	16	10 7 15	8.3	7.9	10 7	6	3.6	5.2
	10 8 53	19	16	10 8 7.9 3.2	2.9	10 8 3.8	2	2.4
	10 9 53	19	16	10 9	6	1.2	1.3	10 9 2.5 0.95 0.97
	10 10 53	20	16	10 10 5.8 1.1 0.78	10 10	2	0.67 0.64
	(a) Born approximation	(b) Our method		(c) Our method+NewtonCG
	Table

  µ 2, §A.1 such that max |Φf E |, |Φf E,1 |, |Φf E,2 |, |Φf * E,1 -Φf * E |, |Φf * E,2 -Φf * E | ≤ 1 3 µ 2, §A.1 , . The estimate (A.4) can be rewritten as 2 Re(Φf * E,l -Φf * E ) Re Φf E (p) + 2 Im(Φf * E,l -Φf * E ) Im Φf E (p) -|Φf E,l | 2 -|Φf E | 2 -|Φf * E,l -Φf * E | 2 ≤ µ 1, §A.1 µ 2, §A.1 E -α-1 2 , l = 1, 2. Re Φf E (p) Im Φf E (p) -b(p, E) ≤ µ 1, §A.1 µ 2, §A.1 E -α-1 2 (A.7) with the matrix M and the vector b from Algorithm 1. Using (3.21), (3.23) and (3.26), we obtain the estimate (3.24).

	(A.5)
	where p ∈ B d 2 √
	(A.6)
	In turn, one can rewrite (A.6) in matrix form as
	2M (p, E) •

[START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]

we get the following estimate:

|Φf E (p) + Φf * E,l (p) -Φf * E (p)| 2 -|Φf E,l (p)| 2 ≤ µ 1, §A.1 µ 2, §A.1 E -α-1 2 , p ∈ B d 2 √ E , E ≥ E * , l = 1, 2, (A.4) for E , E ≥ E *

by convention, j = 1 corresponds to phaseless Born approximation

this potential is given by the Matlab's function peaks

Acknowledgement: TH gratefully acknowledges financial support by DFG through grant CRC 755/C02.

We represent v as follows: where • n,1 is defined in (2.7), and |S d-1 | is the standard Euclidean volume of S d-1 ; see [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]. Using (A.11) we obtain v + (x, r) -

(A.17)

In addition, if r = r(E), where r(E) is the radius of (3.29), then

Using formulas (3.29), (A.17), (A.18), (A. [START_REF] Keiner | Using NFFT 3-a software library for various nonequispaced fast fourier transforms[END_REF], we obtain

Remark A.2. We have that

A.3 Proof of Theorem 3.5

where p ε (p ⊥ , z, t), L j are defined in (3.36), (3.40). Then

Proof of Proposition A.3. Let p ⊥ ∈ R d , p ⊥ • y = 0 and z ∈ Z be fixed. The following estimate holds:

Estimate (A.24) follows from the formula

from the fact that

with the nodes at t = ±1, . . . , ±N , and from the following standard estimate for the Lagrange interpolating polynomial

In addition, the following estimate was proved in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]:

If we replace v(x) by (x • y) 2N v(x), we get

Besides, we also have that

where we have used that ε < 1 2 |y|. Using (A.24), (A.28), (A.29), we get (A.23). Proposition A.3 is proved.

Proposition A.4. Let v and w 1 , w 2 be the same as in (1.1b), (3.6), (3.8). Let v * E be an approximation to v satisfying (2.12). Let U * * E be defined according to (3.27), and let U * N,ε (p, E) be defined according to (3.36), (3.38), (3.39), (3.40). Then

where σ is the same as in (3.6), E * = E * (K, D ext ) is defined according to (2.13), (3.17), and µ k, §X , k ≥ 1, are the constants of Section X.

Using estimates (A.30), (A.32), we get:

)

, for E as in (A.33).

In addition, we have the following estimate proved in [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF] (see the proof of Theorem 2 of [START_REF] Agaltsov | Error estimates for phaseless inverse scattering in the Born approximation at high energies[END_REF]):

In addition, if r = r(E), ε = ε(E), where r(E), ε(E) are defined in (3.44), then

Next, we establish a lower bound for E for which r = r(E), ε = ε(E) satisfy conditions (A.31), (A.33) of Proposition A.4. Note that


Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Raffaella Rizzoni
email: raffaella.rizzoni@unife.it

Sylvie Ronel-Idrissi

Numerical Analysis of Two Non-linear Soft Thin Layers

In a first part, we consider a bar with extremities subject to a given displacement and made by two elastic bodies with linear stress-strain relation separated by an adhesive layer of thickness h. The material of the adhesive is characterized by a non convex (piecewise quadratic) strain energy density with elastic modulus k. After considering the equilibrium problem of the bar and determining the stable and metastable solutions, we let (h, k) tending to zero and we obtain the corresponding asymptotic contact laws, linking the stress to the jump of the displacement at the adhesive interface. The second part of the paper is devoted to the bi-dimensional problem of two elastic bodies separated by a thin soft adhesive. The behaviour of the adhesive is non associated elastic-plastic. As in the first part, we study the asymptotic contact laws.

Introduction

The modelling of adhesive bonding of elastic bodies leads to problems taking into account several parameters, typically the thickness and the rigidity of the adhesive layer. It comes natural to study the limit problem when both these parameters tend to zero. In the limit, the thin layer is replaced by a sharp interface and an asymptotic contact law is obtained, linking the stress to the jump of the displacement at the interface [START_REF] Licht | A modelling of elastic adhesive bonded joints[END_REF][START_REF] Lebon | Numerical study of soft adhesively bonded joints in finite elasticity[END_REF]. Different types of material behaviours have been considered for the adherent and the adhesive bodies. In this paper, two kind of material are considered. We present the first results of a study still in progress in which we consider an elastic bar composed by two adherent parts with linear stress-strain curve separated by an adhesive thin layer of thickness h with non convex energy density and rigidity k. We first discuss the stability and the metastability of the equilibrium configurations when a given relative displacement is imposed to the ends of the bar. Due to the non convexity of the energy, there are multiple metastable configurations made by a mixture of two phases, each phase corresponding to one ascending branch of the stress-strain curve. We then study the limit problem when the pair (h, k) tends to zero and discuss the asymptotic contact laws corresponding to the stable and metastable equilibrium configurations of the bar. The second part is devoted to some numerical results obtained considering an adhesive obeying to a non associated elastic plastic behaviour.

One of the motivations of this analysis is the complexity of the finite element modelling due to the strongly non-linear behaviour of the thin layer. This behaviour implies a very large number of computations. The cost is amplified by the thickness and the stiffness of the layer which are very small compared with the substrata characteristics. These low values increase the ill-conditioning of the tangent matrices and thus, the precision of the computations.

First Part: Non-convex Behaviour

The Equilibrium Problem

Consider a one-dimensional bar occupying a reference unstressed configuration 0 ≤ x ≤ L from which we consider displacement fields u = u(x) continuous with piecewise continuous derivative. The bar is fixed at the extremity x = 0, and subject to a prescribed displacement δ > 0 at x = L, so the following boundary conditions are assigned u(0) = 0, u(L) = δ. Let U be the class of displacements continuous with piecewise continuous derivative, and satisfying the boundary conditions. The bar is assumed to be made by two different elastic materials. In interval 0 ≤ x ≤ L it is composed by an adhesive layer characterized by the piecewise quadratic stored energy density

w a (ξ) = min ξ>0 k 2 ξ 2 ; k 2 ξ 2 -αξ + αΔ (1)
where the modulus k, and the material parameters α, Δ, are positive constants. The derivative of the energy density w a has two ascending branch, each one corresponding to a different phase. In the interval L ≤ x ≤ L the bar is composed by a material with quadratic stored energy density

w b (ξ) = K 2 ξ 2 (2)
The total energy of the bar is given by

E(u) = L 0 w a (u )dx + L L w b (u )dx (3)
with u ∈ U. We say that ū ∈ U corresponds to

| u(x) -ū(x) | + sup x∈(0,L) | u (x) -ū (x) |< η for some η > 0
For ū ∈ U to be a global or a relative minimizer, it is necessary that the first variation of E(u) vanishes at ū :

d dh E(ū + hv) h=0 = 0 (4)
for all admissible perturbations v : [0, L] → IR continuous with piecewise continuous gradient and such that v(L) = 0 = v(0). This leads to the system of equilibrium equations:

w a,ξ (ū ) = σ x ∈ (0, L) w b,ξ (ū ) = σ x ∈ ( L, L) ( 5 )
with σ the (constant) stress in the bar. In view of ( 1) and (2), this system admits the following solutions [START_REF] Lebon | Aymptotic analysis of soft thin layers with nonconvex energy[END_REF]:

i) If 0 ≤ σ < kΔ -α, then, (5) are satisfied by ū = σx/k x ∈ [0, L] σx/K + σ L(1/k -1/K) x ∈ ( L, L] (6)
with σ = k δL -1 , and k-1 = k -1 +(1-)K -1 . By using the expression of σ, it is easy to see that this solution is possibile whenever 0 ≤ δ < (kΔα)L k-1 . This solution is said to correspond to a low-strain monophase configuration because the strain in the adhesive layer (0, L) lies in the first ascending branch of the stress-strain curve w a . ii) If kΔα ≤ σ < kΔ, then, the adhesive layer is a mixture of two phases; a low-strain phase with strain on the first ascending branch of w a , and a high-strain phase with strain on the second ascending branch. Let λ ∈ (0, 1) denote the volume fraction of the low-strain solution. As λ varies in [0, 1], we obtain a a family of equilibrium solutions:

ū = ⎧ ⎨ ⎩ σx/k x ∈ [0, λ L] [(σ + α)x] /k -α λL/k x ∈ (λ L, L] σx/K + L{[σ + (1 -λ)α]/k -σ/K)} x ∈ ( L, L] (7)
with

σ = k δ L -(1 -λ) α k (8)
By using this expression for σ, we find that the two-phase solution exists whenever

(kΔ -α) k-1 + α(1 -λ)k -1 L ≤ δ ≤ kΔ k-1 + α(1 -λ)k -1 L (9)
Because the volume fraction λ ranges over the interval 0 < λ < 1, this condition can be further extended as follows

(kΔ -α) k-1 L ≤ δ ≤ (kΔ k-1 + αk -1 )L (10)
Thus, given any λ ∈ (0, 1), if the prescribed elongation δ lies in the above range, then in the adhesive there exists a two-phase solution involving a mixture of both phases. iii) If σ ≥ kΔ then, the solution is

ū = (σ + α)x/k x ∈ [0, L] σx/Kx + L[(σ + α)/k -σ/K] x ∈ ( L, L] (11)
with σ = k (δL -1α k -1 ). This solution, possibile whenever δ > kΔL k-1 + α Lk -1 , corresponds to a high-strain monophase configuration because the strain in (0, L) lies in the second ascending branch of the stress-strain curve of the adhesive.

Note that a discussion about the metastability of the equilibrium solutions is given in [START_REF] Lebon | Aymptotic analysis of soft thin layers with nonconvex energy[END_REF]. As shown in [START_REF] Lebon | Aymptotic analysis of soft thin layers with nonconvex energy[END_REF] a study of the global stability of the equilibrium solutions restricts the solutions to those corresponding to stable configurations of the bar. This has the effect to select the following configurations:

a) If 0 ≤ δ < (kΔ -α/2)L k-1 , then the low strain monophase solution (6) is energy minimizer. b) If (kΔ -α/2)L k-1 ≤ δ < (kΔ -α/2)L k-1 + αLk -1 ,
then it is minimizer the mixture of low and high strain phases (7) associated with the Maxwell stress σ M = kΔα/2. By substituting σ M into (8), we obtain that in stable configurations the volume fraction of the low strain phase is uniquely determined for the given δ:

λ sta = 1 + K α k k Δ - α 2 k - δ L (12)
Note that when δ = (kΔα/2)L k-1 we have λ sta = 1: this means that the adhesive is still in the low strain phase. When δ = (kΔα/2)L k-1 + αLk -1 , then and λ sta = 0 and the material of the adhesive has completed the transformation from the low strain phase to the high strain phase. c) If δ ≥ (kΔ-α/2)L k-1 + αLk -1 , then the high strain monophase solution [START_REF] Sagan | An introduction to the calculus of variations[END_REF] is energy minimizer.

Contact Laws

We now consider the case in which both the thickness and the elastic modulus of the adhesive are very small compared to the corresponding quantities of the adherent elastic bodies. We would then lead to study the preceding equilibrium solutions when both the parameters ( , k) tend to zero and the thin adhesive layer is replaced by a sharp interface. Accordingly, we set

k = k 0 (13)
and we rescale the material parameters so as

Δ = Δ 0 -1 α = α 0 ( 14
)
with k 0 , α 0 and Δ 0 independent on . Our aim here is to study the relations between the limits Σ = lim

→0 + σ [u] = lim →0 + ū( L) ( 15
)
where σ and ū correspond to the stable equilibrium configurations determined in the previous section. Our aim here is to identify an asymptotic contact law, linking the limit stress in the bar (Σ) to the jump of the displacement at the adhesive interface ([u]). This contact law describes the limit behaviour of the adhesive. By substituting (13) and ( 14) into the expressions of σ and ū listed in Section 2.1 in (a), (b) and (c), taking the limit → 0 + and eliminating δ between Σ and [u] we obtain the following contact law corresponding to stable configurations of the original equilibrium problem.

Σ sta = ⎧ ⎪ ⎨ ⎪ ⎩ k 0 [u] L 0 ≤ [u] < L(Δ 0 -α0 2k0 ) (k 0 Δ 0 -α0 2 ) L(k 0 Δ 0 -α0 2k0 ) ≤ [u] < L(k 0 Δ 0 + α0 2k0 ) k 0 [u] L -α 0 [u] ≥ L(k 0 Δ 0 + α0 2k0 ) (16)
In the same way, considering the σ and ū as in (i), (ii) and (iii), we can calculate the limit contact law corresponding to metastable configurations. Note that this law turns out to be undefined, due to the lack of information associated to non uniqueness of metastable equilibrium solutions. Indeed, we find

Σ meta = ⎧ ⎪ ⎨ ⎪ ⎩ k 0 [u] L 0 ≤ [u] < L(Δ 0 -α0 k0 ) Σ (k 0 Δ 0 -α0 k0 ) ≤ [u] < L(k 0 Δ 0 + α0 k0 ) k 0 [u] L -α 0 [u] ≥ L(Δ 0 + α0 k0 ) (17)
with Σ taking a value in [k 0 Δ 0α 0 , k 0 Δ 0 + α 0 ]. Therefore, the metastable configurations give rise to a multiplicity of contact laws, all contained in the interior of the parallelogram shown in Fig. 1 for

(k 0 Δ 0 -α 0 /k 0 ) ≤ [u] < L(k 0 Δ 0 + α 0 /k 0 ).

The Role of Nucleation and Kinetics

A way rule out the stress indeterminacy of metastable equilibrium solutions is to assign a criterion to select a path between metastable configurations. In [START_REF] Lebon | Asymptotic Study on a Soft Thin Layer: The Non-Convex Case[END_REF][START_REF] Lebon | Aymptotic analysis of soft thin layers with nonconvex energy[END_REF], authors suggest to follow Aberayatne and Knowles [START_REF] Aberayatne | On the driving traction acting on a surface of strain discontinuity in a continuum[END_REF][START_REF] Aberayatne | Kinetic relations and the propagation of phase boundaries in solids[END_REF][START_REF] Aberayatne | Strain-energy functions with multiple local minima: Modeling phase transformations using finite thermoelasticity[END_REF] for introducing a nucleation and a kinetic condition into the equilibrium problem of the composite bar. Here we just sketch a description of these conditions. In the context of the dynamical problem of the composite bar, Aberayatne and Knowles assume that the nucleation and the propagation of a phase boundary is governed by a relation between the normal speed on the phase boundary and the driving force f acting on it. As discussed in [START_REF] Aberayatne | Strain-energy functions with multiple local minima: Modeling phase transformations using finite thermoelasticity[END_REF], one could consider a quasi-static context in which inertial effects are neglected and the governing problem is still the equilibrium problem for the composite bar except that δ, and thus ū depend on the time as a parameter. The totality of the solutions are again described by (i), (ii) and (iii). Consider now the dissipation inequality, which states that the rate of change of the total energy (the sum of potential energy and the kinetic energy) minus the power of the work done at the boundary must be non negative for each equilibrium displacement. In the present context, the total energy coincides with the potential energy so the dissipation can be written as

D := dE(ū) dt -σ dδ dt ≥ 0 (18)
with ū given by [START_REF] Aberayatne | Strain-energy functions with multiple local minima: Modeling phase transformations using finite thermoelasticity[END_REF]. By substituting (7) into the dissipation inequality above and differentiating, we obtain

D = ṡ w a σ k -w a σ + α k + σα k ( 19
)
with ṡ = (dλ/dδ)(dδ/dt) L playing the role of the speed of the phase boundary at x = λ L. The driving force f acting on the phase boundary is defined to be the quantity which multiplicates ṡ:

f := w a σ k -w a σ + α k + σα k = α k (σ -σ M ) (20)
Then, accordingly to Aberayatne and Knowles, a nucleation condition can be stated in the following form: a phase boundary nucleates provided that the driving force on it soon after it nucleates exceeds some critical value f nuc (> 0), i.e. f ≥ f nuc . This condition determines the relative displacement δ at which the point (σ, u ( L)) leaves the first ascending branch of the non monotone stress-strain curve w a,ξ . Once that a phase boundary has nucleated, it propagates according to the kinetic (or evolution) condition, which relates the driving force with the "speed" of the phase boundary

f = φ( ṡ) (21)
with zφ(z) ≥ 0 as enforced by the dissipation inequality (18). The condition (21) makes a selection among the metastable solutions and uniquely determines a path in the stress-strain plane between them. Different paths arise from different types of kinetic function φ.

Now we consider the classical linear kinetics φ(z) = mz and the pinning kinetics φ(z) = m (ax 2b) + , which models the presence of defects that slow down the phase boundary motion [START_REF] Aberayatne | Strain-energy functions with multiple local minima: Modeling phase transformations using finite thermoelasticity[END_REF]. We obtain the evolution of the parameter s

s(t) = (s(0) + K 2 )e -K3t + K 1 t -K 2 (linear) m (a ṡ2 -b) = K 4 + K 5 δt/L + K 6 s (pinning) (22)
The constants K i are given. We suppose that m = m 0 -1 . The contact laws obtained in the metastable domain are:

Σ = k 0 [u]/L -s 0 α 0 s 0 (t) = (s 0 (0) + L 2 )e -L3t + L 1 t -L 2 (linear) m 0 (a ṡ2 0 -b) = L 4 + L 5 Δ 0 t/L + L 6 s 0 (pinning) (23)
The constants L i are given. The contact laws are presented in Fig. 1.

Part 2: Non-associated Elastic Plastic Behaviour

The objective of this part of the paper is to analyse soft thin layers in order to replace them by interface laws. The aim of our work is to study nonlinear soft materials; especially, we focus on the case of non-associated elastoplastic materials of Mohr-Coulomb and Drucker-Prager kind. In [START_REF] Lebon | Asymptotic analysis of Mohr-Coulomb and Drucker-Prager soft thin layers[END_REF], we have obtained, using matched asymptotic expansions, contact laws of Coulomb kind. In bidimensional cases we have:

Σ.n = K 0 ([u e ] -[u p ]) |τ | ≤ C -tan(ψ)σ n [u] p = -kτ Mohr-Coulomb |σ d | ≤ C -tan(ψ)σ m [u] p = -kσ d n Drucker-Prager (24)
The constants C and ψ are material coefficients. We use similar notations than in Section 2. The indices e and p correspond to the elastic and the plastic parts of the displacement. The vector τ corresponds to the direction of the shear stress, σ n is the normal stress. The tensor σ d is the deviatoric part of the stress, σ m is the pressure. Note that these laws depend on all the stress tensor (non local problem). In this part of the paper we analyse the influence of each component of the stress tensor in order to simplify the previous laws. We have shown in the theoretical part, that in the limit problem it is necessary to solve a local problem coupled with the global one. In the local problem, there are two significant quantities which intervene: the plastic yield and the plastic flow. The aim of this analyse is to quantify the level of each term in the plastic yield. In fact, we want to analyse if it is possible to replace the "real" plastic yield in which all the terms of the where σ m s and σ d s are the normal and the deviatoric parts of the tensor σn ⊗ s n. The example analysed is a dovetail assembly (Fig. 2). We compare this plastic yield with the simplified one using only the terms corresponding to the stress vector.

Figure 3 shows the differences between the two plastic yields. We observe a very low difference for the this example and that this difference does not modify the initiation of the plastic process. The gap is maximum in the elastic zone but generally remains lower than 5%. In the plastic zones, this gap decreases to 2%. As a conclusion of this study, we have shown that our simplification is valid and that it is possible to work only with the stress vector for the computation of the plastic yield. As a conclusion, the numerical results obtained in this section show that the local problem introduced in the theoretical study can be neglected, that is to say that the interface law can be written only in terms of stress vector. We obtain a compliance law (regularized Coulomb law), well known in the literature.

Fig. 1

 1 Fig.1Contact law corresponding to stable configurations (bold line), domain of contact laws corresponding to metastable configurations (shaded region) and three metastable configurations (linear, pinning, DC)

Fig. 2

 2 Fig.2The example of a dovetail assembly

Fig. 3

 3 Fig. 3 Exact and simplified plastic yields