An iterative approach to monochromatic phaseless inverse scattering

This paper is concerned with the inverse problem to recover a compactly supported Schr{\"o}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 (1.1a) Here ∆ is the standard Laplacian in x, and the potential function v is assumed to satisfy (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 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., [9,28]: |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 [4,5,6,7,8,9,12,13,18,17,22,23,24,25,28] and references therein), much less studies were performed in the phaseless case (see [1,7,20,21,30,26,27,29]). 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., [29]. In the present work we continue studies of [29,1] 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 = wl 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 (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 [28], and first numerical tests were reported in [5].
Studies on Problem 1 in dimension d = 1 for L = 1 were started in [3], 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 [11]. Studies on Problem 1 in dimension d ≥ 2 were started in [29] and continued recently in [1]. The key result of [29] 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 [1] 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 [1] 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 (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 is surjective. Here and in the following B d r denotes the closed ball 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 γ : 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 2 √ E the setΦ −1 (p) is a piecewise smooth manifold of size Φ −1 (p) and define the averaging operator Using this mapping we can define an approximation u E to v on D by Let W n,1 (R d ) denote the Sobolev space of n-times smooth functions in the sense of L 1 (R d ): If v ∈ W n,1 (R d ), n > d, in addition to the initial assumptions (1.1b), then u E satisfies the error bound for d ≥ 2; see, for example, [28]. Essential improvements of the approximation u E in (2.6) were achieved in [23,24,28]. In particular, formula (2.8) was principally improved in [28] by constructing iteratively nonlinear approximate reconstructions u as j → +∞, (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 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 [2], is essential for studies on direct scattering (see, e.g., [8] ( §29), [28], and references therein) and will also be crucial for our analysis: (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 For inverse scattering with phase information the iterative step of [28] 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 for some s > 1 2 , where a 0 (d, s) is the constant of (2.11). (2.13) 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), k − l = p without reference potential: F = |f (k, l)| 2 scattering amplitudes with reference potentials: 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 [1].
Background potentials of type A: The first simplest case analyzed in [1] is for some w ∈ C(R d ) such that w = w, w(x) = 0 for |x| > R, and 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 [1]. One can see that Background potentials of type B: The second simplest case analyzed in [1] is 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 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 for some fixed τ ∈ (0, 1], (3.11) to an approximate reconstruction in the configuration space is as follows (see [1,Theorem 1]): 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: w1, w2 (p, E) dp, for some τ ∈ (0, 1], Z ε y is defined in (3.10), and 14) 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 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 where D, Ω 1 , Ω 2 are the domains in (1.1b) and (1.5).

23)
for some fixed δ > 0. Then: 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).
Let v * E be an approximation to v satisfying (2.12) for some A > 0, α ≥ 0, K > 0 and 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 . 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, where σ is the constant of (3.6). Let Remark 3.4. The following three conditions are equivalent:

32)
where β is defined in (3.29), n > d, α > 0, and σ > d. Iterations Let v, w 1 , w 2 satisfy the assumptions of Theorem 3.2. Let u E is similar but does not coincide with the approximate reconstruction of Theorem 1 of [1]; see formulas (3.12), (3.11) of the present work. In particular, we have that (3.33) Then, applying the iterative step described above in this subsection we construct nonlinear approximate reconstructions u The approximations u E of (2.9) for phased inverse scattering. In a similar way with (2.10), 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 (3.36) Note that for any p ∈ Z ε y there exists the unique triple (p ⊥ , z, t) such that In addition to U * * of (3.27), we also define under the assumptions that 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: Note also that if assumptions (3.39) are valid for some Under assumptions (3.6), (3.8), the iterative step is realized as follows.
Theorem 3.5 is proved in Section A.3.
where ζ * is defined by (3.20), and also Remark 3.7. The following three conditions are equivalent: 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: Iterations. Let v, w 1 , w 2 satisfy the assumptions of Theorem 3.5. Let u E is similar, but does not coincide with the approximate reconstruction of Theorem 2 of [1]; see formulas (3.16), (3.13) of the present article. In particular, we have that . (3.50) Then, applying the iterative step described above in this subsection we construct nonlinear approximate reconstructions u (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), as j → +∞, 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 .

Numerical experiments 4.1 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 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: 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 [19]. The definitions of the grids Γ s , Γ m , Γ e and of the Fourier transform matrix A are summarized in Algorithm 2. 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 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 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.

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 wl, 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. 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 [31,15].
The iterative step proceeds by computing the discrete functionf appr E :Γ m ∪ Γ e → C, which is the discrete analog of function f appr E of (3.19), as well as the function U * * :Γ m ∪ Γ e → C, which is analogous to function U * * E of (3.27). This computation starts by initializingf appr E and U * * by the zero functions and continues as follows: 1. For each point p ∈Γ m such that |p| < r j , where j ≥ 2 is the current iteration number 4 , 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. for p ∈ Γ m such that |p| < r j do Choose (l(p), l(p)) ∈ argmax (l,l)∈P Re wl ,p Im w l,p − Im wl ,p Re w l,p ;

Evaluate
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 . 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 , σ l (k m ) := n : kn=km (F † l ) n (see [16] 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 Var(F l ) m = t l (k m ) −1/2 (F † l ) m . For 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 N −1/2 p ) is a better characterization of the noise level. As the proposed method only yields an approximate solution at fixed energy even 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 [14]) 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 [14]).
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 (a) and three type-B potentials based on the Wendland functions with k = 1.  Table 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.
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  Table 2: relative L ∞ reconstruction errors in percents for different choices of background potentials. 2A means two background potentials of type A, whereas 2B and 3B refers to two or three background potentials of type B, respectively. (R) refers to indicator functions of rectangles as in Fig. 3, whereas (W) refers to Wendland functions (4.7) with k = 1 with similar scaling. All simulations are performed for E = 15 2 and N p = 10 9 .
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 [14]. 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 [14]). 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: where µ 1, §A.1 > 0. The proof of estimates (A.1) is based on the following inequalities: for l = 1, 2. Applying Lemma 2.1 to v, v * E and to v l , v * E,l , respectively, we get the estimates where 1 2 µ 1, §A.1 > 0 is the constant in the right-hand side of formula (3.14) of [28], p ∈ B d Next, using (A.1) we get the following estimate: In turn, one can rewrite (A.6) in matrix form as 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). Using the formula for the inverse matrix, we also get the following equality: E 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 Let v * E be an approximation to v satisfying (2.12). Let U * * E (p) be defined according to (3.27), 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.
Proof of Proposition A.1. Due to (2.13), (3.27) and Lemma 3.1, we have Besides, in view of (3.9), we also have that We represent v as follows: (A.14) Since v ∈ W n,1 (R d ), n > d, we have where · n,1 is defined in (2.7), and |S d−1 | is the standard Euclidean volume of S d−1 ; see [1].
where |B d 1 | denotes the standard Euclidean volume of B d 1 . It follows from (A.14), (A.15), and (A. 16 (A. 17) In addition, if r = r(E), where r(E) is the radius of (3.29), then Theorem 3.2 is proved.
In addition, we have the following estimate proved in [1] (see the proof of Theorem 2 of [1]):