#  $Z_{1}$-oscillation in the retarding force of metals for slow ions:
Comparative study of theoretical modelings

I. Nagy  Department of Theoretical Physics, Institute of Physics,
Budapest University of Technology and Economics,
H-1521 Budapest, Hungary  Donostia International Physics Center, P. Manuel de
Lardizabal 4,
E-20018 San Sebastián, Spain  I. Aldazabal  Centro de Física de Materiales
(CSIC-UPV/EHU)-MPC, P. Manuel de Lardizabal 5,
E-20018 San Sebastián, Spain  Donostia International Physics Center, P. Manuel
de Lardizabal 4,
E-20018 San Sebastián, Spain

###### Abstract

Theoretical calculations, applying elements of the orbital-based density
functional method for self-consistent screening in an electron gas to the
stopping power of metals for slow projectiles are reviewed, by focusing on
differences in modelings. New, two-channel-based results on the average
retarding force are presented. All theoretical results are compared with
pioneering experimental data, obtained for carbon target, on
$Z_{1}$-oscillation in stopping at random collisional condition.

###### pacs:

34.50.Bw

##  I Introduction and Motivation

According to the profound classification of Hawking Hawking10 , a modeling of
reality in theoretical physics is good if it contains few adjustable elements,
agrees with and explain existing observations, and makes detailed predictions
about future observations. For instance, within our personal selection, in the
famous treatment of Landau on extended, interacting Fermi systems the
determination of the low-energy excitation spectrum requires the introduction
of just one adjustable element Kadanoff62 ; Migdal77 , the quasiparticle
effective mass ($m^{*}$).

Hawking adopted a view that we can call model-dependent realism, which
provides a framework, generally of mathematical nature, with which to
interpret physical phenomena. But different theories can successfully describe
the same phenomenon Hawking10 . Moreover, many theories that proven successful
were later replaced by other ones based on new concepts.

Motivated by such a profound view on model-dependent realism Hawking10 , here
we turn to a particular phenomenon well-known in different subfield of nature.
It is the energy loss of charged projectiles Landau58  in matter, and the
characterization of the slowing down by an average retarding force experienced
by heavy intruders. Evidently, the phenomenon prescribes a quantum (wave)
mechanical consideration of electrons in the target material. Furthermore, a
retarding-force interpretation of the energy loss per unit path length
($dE/dx$) for ions requires the calculation of expectation values of the force
operator (gradient of an external potential) considering complete sets for
electron states in order to quantify an observable.

Pioneering experiments Ormrod63 ; Ormrod64  with slow ($v=0.41$) ions with
$Z_{1}\in[1,20]$ on solid carbon target ($Z_{2}=6$, with four valence
electrons) at random collisional condition revealed a pronounced
$Z_{1}$-oscillation in the observable stopping power (we use Hartree atomic
units throughout this work). For an account on $Z_{1}$-oscillation in gases we
refer to Teplova10 .

Our present, quantum mechanical, understanding of $Z_{1}$-oscillations in
metallic targets with slow charged projectiles dates back to Puska83 ;
Echenique86 ; Nagy89 , as we will outline in the next Section. New elements to
these single-channel theoretical results (obtained for a Kohn-Sham-type
degenerate electron gas model) were already considered and quantified in
Arista94 ; Pitarke05  and, quite recently, in Nagy20  via two (phase-shift-
based) nonlinear channels to the average retarding force. For a discussion of
statistical, quasiclassical theoretical attempts Firsov59 ; Lindhard61  on the
electronic stopping power we refer to a summary Sigmund08 . Our work here is
dedicated to a comparative study on different modelings Puska83 ; Echenique86
; Nagy89 ; Arista94 ; Pitarke05 ; Nagy20  of the observed $Z_{1}$-oscillations
in carbon.

##  II Comparison and discussion

We follow here a time-ordered path in order to outline the above-selected
Puska83 ; Echenique86 ; Nagy89 ; Arista94 ; Pitarke05 ; Nagy20  modelings on
low-velocity ($v$) stopping power. In order to concentrate on concepts, and
associated new adjustable elements, we will avoid extensive (repeated)
mathematical details (where it is possible without distorting the reasoning)
by referring for them to the original papers.

In a scattering interpretation the energy loss per unit path length is a
simple product

$\frac{dE}{dx}\,=n_{0}\,[(vv_{F})\,\sigma_{tr}(v_{F})]$ (1)

where $\sigma_{tr}$ is the transport (momentum transfer) integrated cross
section. The electron number density (a part of constituents of metals) is
$n_{0}$ and the velocity in the last occupied one-electron momentum state
(plane wave) is $v_{F}$. In this semiclassical interpretation $\sigma_{tr}$
can be based on a probability (instead of a probability amplitude) and is
defined by

$\sigma_{tr}\,=\,2\pi\,\int_{0}^{\pi}\,d\theta\,\sin\theta\,(1-\cos\theta)\,d\sigma(\theta)$
(2)

in terms of the differential cross section $d\sigma(\theta)$. Quantum
statistics for fermions is encoded in $v_{F}\equiv{(3\pi^{3}n_{0})^{1/3}}$.
The product $vv_{F}$ is a simple arithmetic average of energy change
($\omega$) in two-particle (binary) collision, i.e.,
$<\omega>=v[(2v_{F}+0)/2]=vv_{F}$.

The wave character of system electrons appears in the quantum mechanical
description of $\sigma_{tr}$ by using a complex scattering amplitude
$[f(\theta,v_{F})]$ to $d\sigma\propto{|f(\theta,v_{F})|^{2}}$. Thus, in terms
of Bessel phase shifts in partial waves for regularized (screened) interaction
one arrives at

$\sigma_{tr}(v_{F})=\frac{4\pi}{v_{F}^{2}}\,\sum_{l=0}^{\infty}\,(l+1)\,\sin^{2}[\delta_{l}(v_{F})-\delta_{l+1}(v_{F})].$
(3)

Pioneering papers Puska83 ; Echenique86  are based on the above-outlined
kinetic framework. Their remarkable contribution (an adjustable element) was
to apply phase shift outputs of the self-consistent procedure of auxiliary-
orbital-based (Kohn-Sham) DFT for screening of an embedded bare charge
$Z_{1}$, independently of incoming-ion charge states. In such a buildup-
procedure of complete screening, where we respect the impact of a metallic
medium in all one-electron states a priori, one has the Friedel sum rule of
neutrality with charge $Q$ on Mahan81  impurity

$Q\equiv{Z_{1}}\,=\,\frac{2}{\pi}\,\sum_{l=0}^{\infty}\,(2l+1)\,\delta_{l}(v_{F}),$
(4)

and the Levinson theorem of scattering theory with local Burke77  interaction
is satisfied.

Notice, that the possible role of relative (r) kinematics within an
independent-particle-scattering picture beyond the impurity limit
($v\rightarrow{0}$) was investigated in thesis-papers Remi06 ; Remi07  for
$v\in{(0,v_{F})}$, by considering prefixed (at $v\rightarrow{0}$) DFT
potentials Puska83 ; Echenique86 ; Nagy89  as plausible ones. Notable
reduction in ratios, obtained at $v\rightarrow{0}$, of maximum/minimum in
$Z_{1}$-oscillation was observed for $v\rightarrow{v_{F}}$. As expected on
physical grounds, and due to mathematics in
$\sigma_{tr}(v_{r})=(4\pi/v_{r}^{2})\sum_{l=0}^{\infty}(l+1)\sin^{2}[\delta_{l}(v_{r})-\delta_{l+1}(v_{r})]$,
boundary values in oscillations are the most sensitive to tuning the relative
(r) wave number. The comparison with surface-experimental data Winter03 ,
which are almost free from complications due to lattice-ions, gave an improved
character for $v=0.5$ at Al surface ($v_{F}\simeq{0.87}$) in
$Z_{1}$-oscillation.

Of course, as we know from an insightful book Peierls79 , the above-discussed
(mean-field) consideration of screening is one possibility. An alternative way
uses the consideration that the matrix elements of the external potential
which connect occupied states described by Slater determinant has no effect on
the many-body wave function since their effect cancels out in a determinant
wave function. This fact suggests starting from the Schrödinger equation in
which the potential is replaced by a truncated potential which is defined like
the Bethe-Goldstone potential, except that it is concerned with one particle
only, and therefore simpler. However, such a pseudopotential is nonlocal since
in its definition for a given value of the space coordinate the value of one-
particle wave function at other space-points enter.

We are not aware of such an implementation for screening to stopping
calculations via phase shifts, but we should keep in mind this alternative way
for cases of partially-ionized projectiles since there are similar values in
$\sigma_{tr}$ for small phase shift differences or differences close to $\pi$.
With nonlocal Burke77  interaction the modified Levinson theorem becomes

$[\delta_{l}(0)-\delta_{l}(\infty)]\,=\,\pi(n_{b}+n_{p})$ (5)

where $n_{b}$ is the number of true bound states and $n_{p}$ is the number of
”bound states” excluded by the Pauli exclusion principle for fermions. For
instance, in electron-helium atom scattering in vacuum $n_{b}=0$ and
$n_{p}=1$, thus $\delta_{0}\rightarrow{\pi}$ since there is no stable negative
ion. In orbital-based construction of DFT with $Z_{1}=2$ one has $n_{b}=1$ and
$n_{p}=0$ thus $\delta_{0}\rightarrow{\pi}$ for vanishing scattering wave
number, i.e., at very low density Nagy94  of the electron gas. However, a
similar tendency to $\pi$ does not imply, to our best knowledge, the same
finite values for the scattering length and thus for cross sections related to
observables.

Results based on Puska83 ; Echenique86  for $(1/v)(dE/dx)$, obtained within an
orbital-based (Kohn-Sham) treatment for screening and scattering of an
embedded $Z_{1}$ are plotted in Fig. 1 by dashed curve for $r_{s}\simeq{1.6}$
of the Wigner-Seitz parameter, $r_{s}=[3/(4\pi
n_{0})]^{1/3}=(9\pi/4)^{1/3}/v_{F}$. This parameter represents the valence
electron density of carbon target . It was pointed out in the pioneering
Echenique86  that the data (solid circles) and theory (dashed curve) show the
same overall trend with good quantitative agreement for $Z_{1}<7$. However, as
was emphasized in the other pioneering work Puska83 , experimental values tend
to increase more rapidly as $Z_{1}$ increases. According to the qualitative
discussion of physics in Puska83 , in a kinetic picture, the increase can be
due to several effects. For instance, the ionic radius increases with $Z_{1}$
and therefore heavier ions see a larger effective electron density than the
lighter ones.

The closely related problem of charge inhomogeneity under random collisional
condition in carbon, was addressed and quantified in Arista94 . There, as a
necessary extension of Lindhard61  with hidden details Sigmund08 , an
averaging procedure was defined based on a statistical local-density
approximation, with $[dE/dx](r)$ local inputs constrained by Eq.(4), via
$r_{s}(r)$ as follows

$\langle{\frac{dE}{dx}}\rangle\,\equiv{N_{a}\,\int_{0}^{r_{a}}\,\frac{dE}{dx}(r)\,4\pi
r^{2}dr}.$ (6)

Here $N_{a}$ is the number of atoms per unit volume and $r_{a}$ is the atomic
cell radius. The dash-dotted curve in Fig. 1 shows the such-obtained averaged
results. There is a remarkable improvement especially at around the
experimental minimum and beyond it for higher $Z_{1}$.

However, in a more recent analysis Arista10  on averaging procedures it was
pointed out that one of the shortcomings of the above approach is the fact
that it includes in the integration all the atomic electrons ($r_{min}=0$)
without consideration of the important binding effects that will tend to
cancel the contribution of inner shell(s) to the stopping in the case of slow
ions. In the case of carbon with four valence electrons this shell is the
doubly occupied $1s$ shell. Clearly by excluding this (compact in its
extension) shell from averaging one goes back to the dashed curve. In our
opinion, this alarming observation on the realistic role of target-atom inner
shells remains valid for metallic targets with heavier Arista19  atomic
constituents as well. Intuitively, swift projectiles Nagy02  without bound
electrons around them, say $Z_{1}=\pm{1}$, could be better candidates to a
statistical averaging procedure as the one in Eq.(6).

Figure 1: Stopping power of carbon ($Z_{2}=6$) target ($r_{s}\simeq{1.6}$) for
ions with atomic number $Z_{1}$. Transmission data (denoted by solid circles),
obtained at random collisional condition, are taken from Ormrod63 ; Ormrod64
for $v=0.41$. The dashed curve corresponds to the conventional (single-
channel) result obtained in DFT at embedding condition Puska83 ; Echenique86 .
The dash-dotted, and dash-dot-dot-dashed curves (with new elements beyond this
conventional description) refer to Arista94  and Pitarke05 , respectively. The
solid curve is based on Eqs.(7-9), and refers to our two-channel approach
developed recently Nagy20 .

.

To DFT calculations Puska83 ; Echenique86  the electrons (without impurity)
are considered as noninteracting in their plane-wave Kohn-Sham one-electron
states. The momentum distribution function of them is an ideal Fermi-Dirac
distribution with (invariant) Fermi surface defined from the density $n_{0}$.
However, as many-body methods Kadanoff62 ; Mahan81 ; Pines61  indicate, the
excitation spectrum above such a surface may get modulations (due to dynamical
exchange-correlation Nagy85  between electrons) beyond the one prescribed by
the imaginary part of the noninteracting ($m^{*}=m$) density-density response
function [which is $\propto{\omega}$, independently of the characteristics
($v_{F}$) of the filled ground-state], i.e., of the Lindhard function of the
time-dependent linearized Hartree-type approximation for screening. In Nagy89
we made an attempt to combine static screening nonlinearities with many-body
modulation in the electron-hole excitation spectrum (which should depends on
$m^{*}$, according to Landau’s physics-based Kadanoff62 ; Pines61  intuitive
treatment) responsible for dissipation. Only small changes in stopping,
depending slightly on $Z_{1}\in{[5,40]}$, were observed in Echenique89  for an
electron gas with $r_{s}=2$.

Therefore, it was a surprise when a sophisticated approach Pitarke05  based on
time-dependent density-functional theory (TDDFT), for the prototype model
system of an embedded $Z_{1}$ and an interacting electron gas, within its
local-density approximation gave a much higher modulating effect [a density-
gradient-dependent additive term to Eq.(1)] in $Z_{1}$-oscillation of results
plotted by dashed curve. This prediction is exhibited in Fig. 1 by a dash-dot-
dot-dashed curve. The theoretical enhancement was interpreted tentatively as
excitonic effect in electron-hole pair excitation and good agreement with
experimental data on carbon target (see, Fig. 1) at random collisional
condition was concluded. Later the above-outlined adjustable element of
modeling was reviewed Pitarke07  within the physically more consistent time-
dependent current density functional theory (TDCDFT). Briefly, this new
approximation, due to its tensorial character, can satisfy important
conservation laws and associated sum rules. Owing to such a physical
consistency, the re-adjusted numerical results show only very small changes to
pioneering DFT results, i.e., to the dashed curve in Fig. 1.

Viewing metals as systems of degenerate electron gases (with corresponding
$r_{s}$ parameters) and fixed lattice ions, a recent theoretical attempt
Nagy20 , with the goal to extend Puska83 ; Echenique86 ; Nagy89 , considered
the role of sudden charge-change cycles, as an adjustable new element to the
average retarding force experienced by projectiles in their slowing down under
random collisional conditions. There, the force interpretation was implemented
by using quantum mechanical calculations for matrix elements, instead of
direct probability interpretation behind a differential cross section needed
in the classical definition in Eq.(2).

We believe that, within the wave mechanical description of one-electron
states, such a matrix-element-based Nagy20  treatment with an additional
sudden time-dependent process is a physically consistent one. Clearly, and in
agreement with an enlightening analysis Lindhard71  of an esteemed expert in
the field of stopping, the major (from stopping point of view) change from a
free Fermi gas occurs because the strong Coulomb field around lattice nuclei
allow an increase in the rate of local (sudden) processes of electron jumping
into and away from bound states on the projectile. Band-related further,
presumably $l$-dependent, effects are not addressed here in our projectile-
centered construction of one-particle scattering states.

Technically, we used earlier established works Gaspari72 ; Bonig89 ; Tang98
for screened interaction energy in order to model a new (regularized
Coulombic) adjustable dissipative channel which could model the sudden
ionizing role of lattice ions in the many-body system of electrons at random
collisional condition in transmission. Concretely, in Nagy20  we derived

$\frac{1}{v}\,\frac{dE}{dx}\,=\,\left[Q^{(1)}(v_{F})\,+Q^{(2)}(v_{F})\right],$
(7)

where the two coefficients ($q\neq{0}$) are given by the following expressions

$Q^{(1)}(v_{F})\,=\,\frac{4}{3\pi}\,v_{F}^{2}\,\sum_{l=0}^{\infty}(l+1)\sin^{2}[\delta_{l}(v_{F})-\delta_{l+1}(v_{F})],$
(8)
$Q^{(2)}(v_{F})\,=\,\frac{4}{3\pi}\,q^{2}\,\sum_{l=0}^{\infty}\frac{1}{l+1}\,\sin^{4}\left[\frac{\delta_{l}(v_{F})-\delta_{l+1}(v_{F})}{2}\right].$
(9)

We used phase shifts values from earlier works Puska83 ; Echenique86 ; Nagy89
in this two-channel modeling and the results ($q=1$) obtained from Eq.(7) is
plotted by solid curve in Fig. 1. The agreement with data, especially beyond
the first maximum in oscillation, is very reasonable. Thus, we believe, that
the underlying complete set of scattering states generated in DFT via an
iterative procedure by embedding a bare charge $Z_{1}$ is, a posteriori, the
proper relaxed one even for singly ionized external projectiles. This was
expected on physical grounds by considering the immediate screening action of
the charged many-body system of metals.

In order to provide a template-like figure to future new developments, we do
not plot data (for instance for Al target as in Nagy20  for $Z_{1}\leq{20}$,
at random collisional condition) in our second, illustrative Fig. 2. This
shows, in its clinical form, our prediction with a new adjustable element to
characterize an observable. We suggest further attempts. Especially,
experiments with $Z_{1}\in{[20,40]}$ on Al at random collisional condition are
recommended.

Figure 2: Illustrative theoretical stopping power at $r_{s}=2$ for slow ions
with atomic number $Z_{1}$. The dashed and dotted curves refer to Eq.(8) and
Eq.(9), respectively. Their sum is plotted by solid curve. See the text, and
Nagy20 , for further details on our force-interpretation behind channels.

.

We motivated our two-channel approach for an average retarding force by
referring to sudden local charge changing processes behind $q=1$ at random
collisional situation in metals. Can one characterize such an average within
the framework of an effective single-channel modeling via
$\bar{\sigma}_{tr}(v_{F})$ in Eq.(1)? The rest of this comparative study on
capabilities of adjustable elements is motivated by such a question. Our
discussion below on this (challenging) question is based on quantum physics
but has an intuitive character.

We return to an important message put forward in Lindhard71 . For slow moving
ions there is stationary flow (at far distances) of electrons. It was pointed
out that a total balance in that flow could replace the concept of local
capture-loss processes. But no quantification of this idea was presented.
Later, in a combined paper, we quantified Nagy00  such an idea based on
conservation laws, considering the case of slow He ions where charge-change
may occur. We applied, as constraint to a one-parametric model potential, the
form derived in Zwerger97  for the dipolar backflow ($\mathcal{D}_{b}$)
amplitude in terms of phase shifts solely at the Fermi energy

$D_{b}=\frac{1}{\pi}\sum_{l=0}^{\infty}(2l+1)\sin{(2\bar{\delta_{l}})}+\frac{4}{\pi}\sum_{l=0}^{\infty}(l+1)^{2}\sin{\bar{\delta}_{l}}\sin{\bar{\delta}_{l+1}}\sin{(\bar{\delta}_{l}-\bar{\delta}_{l+1})}$
(10)

For the very clear, pedagogically enlightening discussion of conservation laws
behind the above form we refer to the original detailed research paper. This
expression reduces to the familiar Friedel sum rule of charge-screening
($D_{b}\Rightarrow{Z_{1}}$) for small phase shifts, i.e., in the first-order
Born limit where $(Z_{1}/v_{F})<<{1}$ for the Sommerfeld parameter. In
perturbative cases with small charges or high densities, where we have no
bound states, the second term becomes third-order in the perturbing charge.

Table 1:  Ratios ($\mathcal{R}_{i}$) of He and H stopping powers for the metallic range of $r_{s}$ in a time-ordered representation. The correspondences are $i=ENAR$ Echenique86 , $i=NES$ Nagy00 , and $i=NA$ Nagy20 . The 3rd column refers to the amplitude of the dipolar backflow ($\mathcal{D}_{b}$) with DFT phase shifts based on complete screening (Friedel sum rule) of an embedded charge $Z_{1}=2$. The 5th column shows the dipolar backflow amplitude of a forced calculation Nagy00  at the $\mathcal{D}_{b}\equiv{Z_{1}}=2$ constraint in the longitudinal current of the slowly moving charge. See the text for further details. $r_{s}$ | $\mathcal{R}_{ENAR}$ | $\mathcal{D}_{b}^{ENAR}$ | $\mathcal{R}_{NES}$ | $\mathcal{D}_{b}^{NES}$ | $\mathcal{R}_{NA}$
---|---|---|---|---|---
$1.5$ | $2.44$ | $1.12$ | $2.45$ | $2$ | $2.56$
$2.0$ | $1.59$ | $0.80$ | $2.14$ | $2$ | $2.23$
$3.0$ | $0.81$ | $0.27$ | $1.77$ | $2$ | $2.37$

Table I contains the essence of our message, by considering the problem of
$He^{+}$ and proton [described solely by its $Q^{(1)}$] stopping ratio
$\mathcal{R}$ in metals at low ion velocities. At $r_{s}=1.5$, which
correspond to a comparatively high density, the ratios are similar and all
ratios are below the first-order Born value, i.e., four. At $r_{s}=2$ our
$\mathcal{R}_{NA}$ ratio is in accord with TDDFT calculation in Zeb13 . At low
density, $r_{s}=3$, our ratio is in nice agreement with experimental
prediction Markin09  for $v\leq{0.2}$ and is below the theoretical (TDDFT)
value of Zeb12 , which is about $4.7$, but is essentially larger than
$\mathcal{R}_{ENAR}\simeq{0.8}$ obtained in DFT. We believe that the high
($4.7$) value is due to underestimation of proton stopping in TDDFT
simulation. This fact was stated in Zeb12  devoted to He/H anomaly in Au. For
completeness, we note that at $r_{s}=3$ our two-channel-based He stopping is
roughly equivalent to proton stopping calculated (in DFT) at $r_{s}=1.5$. We
stress that this (separate) statement is also in harmony with experimental
data Markin09  obtained on gold target for low ion velocity, $v\leq{0.2}$, at
which we consider only an $s$-like state to define $r_{s}$. In more
quantitative terms, we have
$[Q^{(1)}+Q^{(2)}](r_{s}=3)|_{Z_{1}=2}\simeq{0.38}$ and
$Q^{(1)}(r_{s}=1.5)|_{Z_{1}=1}\simeq{0.31}$ based on Table I of Nagy20 .

Quite independently of $r_{s}$ we have a reasonable similarity between
$\mathcal{R}_{NES}$ and $\mathcal{R}_{NA}$ values, which may indicate the
applicability of an adjustable element discussed in Lindhard71 , and
quantified in Nagy00 . Here we found by comparison that, in a statistical
sense for random collisional situations, we can model (reinterpret) the local
capture-loss processes by using a constraint in the electronic flow around a
slow heavy ion which represents a longitudinal current as well, not only a
static charge. Clearly, with moving charges in a charged many-body system one
has to satisfy a fundamental constraint prescribed by continuity equation
between induced charges and the longitudinal current. Remarkably, a recent
large-scale orbital-based TDDFT simulation emphasized Correa15  that
application of current-based (not only charge-based) implicit version could be
a promising way for future more sophisticated extensions within the framework
of first-principles methods. In the orbital-free, i.e., explicit, version of
TDDFT an attempt to include current via a dynamical kinetic-energy functional
has already been presented recently Ding18 ; White18 . These modelings of
reality within TDDFT would fit, of course, to Hawking profound classification
outlined in the Introduction.

##  III summary

We have reviewed earlier approaches and presented new results on the average
retarding force of metallic targets for slow projectiles, using pioneering
experimental data on $Z_{1}$-oscillation in carbon at random collision
condition. All theoretical approaches, selected to our comparative work, are
based on different adjustable elements. These are analyzed in details in order
to arrive at established conclusions. We pointed out challenging problems
which need future developments in the important field of ion stopping in
metals.

###### Acknowledgements.

This work was supported partly by Project PID2019-105488GB-I00 of the Spanish
Ministry of Science, Innovation, and Universities, MICINN.

## References

  * (1) S. Hawking and L. Mlodinov, The Grand Design (Bantam Books, London, 2010), Chap. 3.
  * (2) L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics (Perseus Books, Cambridge, 1962), Chap. 11.
  * (3) A. B. Migdal, Qualitative Methods in Quantum Theory (Benjamin, London, 1977), Chap. 5.
  * (4) L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Pergamon Press, London, 1958).
  * (5) J. H. Ormrod and H. E. Duckworth, Can. J. Phys. 41, 1424 (1963).
  * (6) J. H. Ormrod, J. R. Macdonald, and H. E. Duckworth, Can. J. Phys. 43, 275 (1964).
  * (7) Yu. A. Belkova and Ya. A. Teplova, Bull. Russian Acad. Sci. 74, 153 (2010).
  * (8) M. J. Puska and R. M. Nieminen, Phys. Rev. B 27, 6121 (1983).
  * (9) P. M. Echenique, R. M. Nieminen, J. C. Ashley, and R. H. Ritchie, Phys. Rev. A 33, 897 (1986), and references therein.
  * (10) I. Nagy, A. Arnau, and P. M. Echenique, Phys. Rev. A 40, 987 (1989).
  * (11) J. Calera-Rubio, A. Gras-Marti, and N. R. Arista, Nucl. Instrum. Methods B 93, 137 (1994).
  * (12) V. U. Nazarov, J. M. Pitarke, C. S. Kim, and Y. Takada, Phys. Rev. B 71, 121106(R) (2005).
  * (13) I. Nagy and I. Aldazabal, Phys. Rev. A 102, 032814 (2020), and references therein.
  * (14) O. B. Firsov, Soviet Phys. JETP, 36, 1076 (1959).
  * (15) J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961).
  * (16) P. Sigmund, Bull. Russian Acad. Sci. 72, 153 (2008).
  * (17) G. Mahan, Many-Particle Physics (Plenum Press, New York, 1981).
  * (18) P. G. Burke, Potential Scattering in Atomic Physics (Plenum Press, New York, 1977).
  * (19) R. Vincent and I. Nagy, Phys. Rev. B 74, 073302 (2006).
  * (20) R. Vincent and I. Nagy, Nucl. Instrum. Methods B 256, 182 (2007).
  * (21) H. Winter, J. I. Juaristi, I. Nagy, A. Arnau, and P. M. Echenique, Phys. Rev. B 67, 245401 (2003), and references therein.
  * (22) R. Peierls, Surprises in Theoretical Physics (Princeton University Press, Princeton, 1979).
  * (23) I. Nagy and A. Arnau, Phys. Rev. B 49, 9955 (1994).
  * (24) E. A. Figueroa and N. R. Arista, J. Phys.: Condens. Matter 22, 015602 (2010).
  * (25) F. Matias, P. L. Grande, M. Vos, P. Koval, N. E. Koval, and N. R. Arista, Phys. Rev. A 100, 030701(R) (2019).
  * (26) I. Nagy, Phys. Rev. A 65, 014901 (2002), and references therein.
  * (27) D. Pines, The Many-Body Problem (Benjamin, New York, 1961), pp. 260-278.
  * (28) I. Nagy, J. László, and J. Giber, Z. Phys. A 321, 221 (1985).
  * (29) P. M. Echenique, I. Nagy, and A. Arnau, Int. J. Quantum Chem. 23, 521 (1989).
  * (30) V. U. Nazarov, J. M. Pitarke, Y. Takada, G. Vignale, and Y.-C. Chang, Phys. Rev. B 76, 205103 (2007).
  * (31) J. Lindhard and A. Winter, Nucl. Phys. A166, 413 (1971).
  * (32) G. D. Gaspari and B. Györffy, Phys. Rev. Lett. 28, 801 (1972).
  * (33) L. Bönig and K. Schönhammer, Phys. Rev. B 39, 7413 (1989).
  * (34) J.-M. Tang and D. J. Thouless, Phys. Rev. B 58, 14179 (1998).
  * (35) I. Nagy, Ch. Eppacher, and D. Semrad, Phys. Rev. B 62, 5270 (2000), and references therein.
  * (36) W. Zwerger, Phys. Rev. Lett. 79, 5270 (1997).
  * (37) M. Ahsan Zeb, J. Kohanoff, D. Sánchez-Portal, and E. Artacho, Nucl. Instrum. Methods B 303, 59 (2013).
  * (38) S. N. Markin, D. Primetzhofer, M. Spitz, and P. Bauer, Phys. Rev. B 80, 205105 (2009).
  * (39) M. Ashan Zeb, J. Kohanoff, D. Sánchez-Portal, A. Arnau, J. I. Juaristi, and E. Artacho, Phys. Rev. Lett. 108, 225504 (2012).
  * (40) A. Schleife, Y. Kanai, and A. A. Correa, Phys. Rev. B 91, 014306 (2015).
  * (41) Y. H. Ding, A. J. White, S. X. Hu, O. Certik, and L. A. Collins, Phys. Rev. Lett. 121, 145001 (2018).
  * (42) A. J. White, O. Certik, Y. H. Ding, S. X. Hu, and L. A. Collins, Phys. Rev. B 98, 144302 (2018).


# Multiscale Phase Retrieval

David J. Brady
Wyant College of Optical Sciences
University of Arizona
Tucson, Az 85721
<EMAIL_ADDRESS>
&Timothy J. Schulz
Department of Electrical and Computer Engineering
Michigan Technological University
Houghton, MI 49931
<EMAIL_ADDRESS>
&Chengyu Wang
Department of Electrical and Computer Engineering
Duke University
Durham, NC 27705
<EMAIL_ADDRESS>

###### Abstract

While characterization of coherent wavefields is essential to laser, x-ray and
electron imaging, sensors measure the squared magnitude of the field, rather
than the field itself. Holography or phase retrieval must be used to
characterize the field. The need for a reference severely restricts the
utility of holography. Phase retrieval, in contrast, is theoretically
consistent with sensors that directly measure coherent or partially coherent
fields with no prior assumptions. Unfortunately, phase retrieval has not yet
been successfully implemented for large-scale fields. Here we show that both
holography and phase retrieval are capable of quantum-limited coherent signal
estimation and we describe phase retrieval strategies that approach the
quantum limit for >1 megapixel fields. These strategies rely on group testing
using networks of interferometers, such as might be constructed using emerging
integrated photonic, plasmonic and/or metamaterial devices. Phase-sensitive
sensor planes using such devices could eliminate the need both for lenses and
reference signals, creating a path to large aperture diffraction limited laser
imaging.

## Introduction

Characterization of a coherent wave field under the constraint that only the
squared modulus of the wave is detectable has been addressed by holography [1,
2] and by phase retrieval [3, 4, 5]. Holography utilizes a reference
signal—something Gabor called a "coherent background." When the reference is
sufficiently intense, the measured intensities can be linearly processed to
recover the field. In contrast, without the benefit of a reference signal,
phase-retrieval requires estimation of the field from a quadratic form [6].
Here we utilize interconnected arrays of low rank transformations to reduce
computational complexity and enable practical reference-free measurement of
megapixel fields. We specifically show that $5n\log n$ appropriately coded
measurements are sufficient to estimate the phase and amplitude of megapixel
scale signals to within 10% of the quantum limited mean square error (MSE).

Phase retrieval, under various names and methods, has been explored in many
different contexts. Classically defined "phase retrieval" as pioneered in the
above cited studies by Gerchberg, Saxton, Fienup and Gonslaves began with the
challenge of reconstructing a signal from the modulus of its Fourier
transform, but has recently expanded to include analysis of more general
quadratic forms [6]. From another perspective, wavefront estimation from
shearing interferometry [7, 8] has long been used in adaptive optics and
astronomy. Phase estimation also lies at the center of two mode
interferometry, as explored extensively in gravity wave detection and quantum
state analysis [9]. This paper applies algorithms from conventional phase
retrieval to physical measurements related to shearing and two mode
interferometry. The advantage of this approach is that low rank measurements
reduce the physical and computational complexity of phase retrieval while
enabling near quantum limited field estimation. The main contributions of this
paper are:

  1. 1.

to show that the quantum limited lower bound on the mean square error for
estimating a coherent field consisting of $N$ modes is $N$;

  2. 2.

to show that phase retrieval can achieve this bound; and

  3. 3.

to describe and simulate practical phase retrieval measurements that approach
this bound.

While we expect that improved results may be achieved by applying signal
priors and regularization, such methods are not necessary for the results
presented here and are thus reserved for future analysis. Similarly, we
anticipate that improved results might be obtained using prepared quantum
states, but here we limit our attention to semiclassical analysis using a
Poisson noise model.

Phase retrieval consists of estimating the complex field $\mathbf{x}$ from the
magnitude or intensity of the linearly transformed field
$\mathbf{y}=\mathbf{A}\mathbf{x}$ [10]. While early studies focused on
$\mathbf{A}$ as a Fourier kernel, recent studies have modified the forward
model by adding coding elements in the Fourier plane [11, 12, 13]. Fourier
coding is closely related to multiband filtering, which is called "Fourier
ptychography" [14]. Expanding this approach, numerous recent mathematical
studies have considered phase retrieval under the assumption that $\mathbf{A}$
is a random matrix. Estimating $\mathbf{x}$ given the measurements
$|y_{m}|^{2}$ and the matrix $\mathbf{A}$ is known to be NP-hard [15] and is
typically solved by iterative algorithms. While the diversity of algorithms is
large, novel algorithms have not dramatically exceeded the signal fidelity of
classic iterative approaches [16]. However, recent studies have led to better
understanding of the role of the coding matrix $\mathbf{A}$. Design of
measurement matrices has a long history in linear measurement systems, as
evidenced by Hadamard spectroscopy [17] and uniformly redundant arrays [18],
but to date has been little considered for quadratic forms.

Shearing interferometry [19] and wavefront sensors [20] can also be used for
phase retrieval. The matrix $\mathbf{A}$ is related to the gradient operator
in this case, and numerous algorithms have been developed to estimate the
field from measurement of the gradient [21, 22, 23]. Unfortunately,
integration of the gradient is not well conditioned. One can improve
conditioning by making the local sampling more sophisticated [24]. One can
view this approach as subgroup testing using nearest neighbors, here we show
that the use nonlocal groups further improves estimation.

Phase estimation strategies for two mode interferometry have been explored by
many fundamental and applied studies in the context of gravity wave detection
and quantum state analysis [9]. Most such studies do not consider estimation
of the relative amplitudes of the modes. With respect to phase estimation, the
"shot noise limit" yields an estimator deviation inversely proportional to
$\sqrt{\bar{K}}$, where $\bar{K}$ is the expected number of photons in the
phase measurement. Using entangled states of light, it is possible to improve
on shot noise to reach the "Heisenberg" limit, where the deviation is
inversely proportional to $\bar{K}$ [9, 25]. The success of phase measurement
in two mode systems may be viewed as an exitence proof that quantum limited
phase retrieval is possible. Here we use such systems to scale to large scale
fields.

We seek to combine insights from Fourier phase-retrieval, wavefront sensing
and two-mode phase estimation in the design of measurement matrix
$\mathbf{A}$. We are inspired in this regard by more than a half century of
work considering holographic wavefront reconstruction as a coding and
communication problem [26]. Our approach focuses on "multiscale" analysis,
under which networks of low rank interferometers are interconnected to
reconstruct high dimensional fields. We are able to prove that as the
dimension rank of the interferometer and the field increase this approach
asymptotically reaches quantum limited performance. We then show in simulation
that this approach yields error within 10% of the quantum limit for more
practical measurements. We then further extend the multiscale concept by
simulating block-based phase retrieval, which linearizes the computational
cost per pixel at the cost of increasing MSE by a constant factor.

Focal imaging systems rely on the concept of a "focal plane array," e.g. a
generic sensor that can digitize an image field. An analogous device for
coherent field characterization is the vision of this paper. In imagining such
a device, we must consider what would be analogous to a the focal plane’s
pixel, e.g. what is the atomic component of the measurement. We propose that
joint measurement of small subsets of the unknown field provides this atomic
kernel. We find below that arrays of 2 to 5 mode interferometers can be used
to scale phase retrieval to arbitrarily large fields. Appropriate arrays
processing a modest number of modes have been implemented via integrated
optics[27, 28]. Use of such systems for interferometric imaging is described
by Su, et al. [29]. We hope that the mathematical results presented here will
spur further development of such systems. Among many potential applications,
such sensor arrays might enable lensless imaging over large apertures without
the need for a reference beam. This vision is contained in the idea of
"holographic aperture ladar," [30, 31] with the exception that here we suggest
that the "coherent background" of holographic—or whole field — measurement is
not necessary.

## Measurements and Fisher Information

This section addresses the general problem of estimating a coherent field from
the intensities of an affine transformation of the field. In doing so, we
arrive at four significant results: 1) a lower bound for the squared
estimation-error that any system can attain; 2) a necessary and sufficient
condition to determine if a measurement system induces this optimal bound; 3)
a simple proof that phase-quadrature holography is optimal; and 4) an
existence proof that phase-retrieval is asymptotically optimal.

A continuous coherent field $g(u)$ in an aperture plane can be described by a
discrete vector of modal amplitudes $\mathbf{x}$ by projection onto the basis
of orthonormal mode distributions $\\{\psi_{n}(u)\\}$ according to

$g(u)=\sum_{n=1}^{N}x_{n}\psi_{n}(u).$ (1)

The intensity of the field is $\|\mathbf{x}\|^{2}$, and—when using the
semiclassical theory for detection—this intensity specifies the expected
number of photo-detections we associate with the field:

$\displaystyle\bar{K}$ $\displaystyle=\|\mathbf{x}\|^{2}$
$\displaystyle=\sum_{n=1}^{N}\left|x_{n}\right|^{2}.$ (2)

Characterization of the vector of mode amplitudes, $\mathbf{x}$, requires
intensity measurements that encode both its amplitude and phase. Holography
and phase retrieval accomplish this by measuring the intensity of an affine
transformation of the modes:

$\mathbf{y}={\boldsymbol{\rho}}+\mathbf{A}\mathbf{x},$ (3)

where $\mathbf{A}$ is an $M\times N$ orthonormal matrix with $M>N$, and
${\boldsymbol{\rho}}$ is an $M$-element non-negative reference signal.
Orthonormality of $\mathbf{A}$ ensures that the measurements conserve the
field’s intensity.

The intensities $|y_{m}|^{2}$ may be detected consistent with the
semiclassical theory for detection with

$\Pr\left[\mathbf{d}=\mathbf{k}|\mathbf{x}\right]=\prod_{m=1}^{M}e^{-\left|y_{m}\right|^{2}}\left(\left|y_{m}\right|^{2}\right)^{k_{m}}\frac{1}{k_{m}!}.$
(4)

We assume that the $m$th row of the measurement matrix $\mathbf{A}$ is
$\mathbf{a}_{m}^{\dagger}$, where we use the superscript † to denote the
conjugate transpose of a matrix or vector, and the superscript T to denote the
standard transpose. If we parameterize the complex-valued mode vector by its
real and imaginary parts:

$\mathbf{x}=\mathbf{r}+j\mathbf{i},$ (5)

then the squared error for an estimate $\widehat{\mathbf{x}}$ is

$\displaystyle{\cal E}^{2}$
$\displaystyle=\sum_{n=1}^{N}\left|\widehat{x}_{n}-x_{n}\right|^{2}$
$\displaystyle=\sum_{n=1}^{N}\left(\widehat{r}_{n}-r_{n}\right)^{2}+\sum_{n=1}^{N}\left(\widehat{i}_{n}-i_{n}\right)^{2},$
(6)

and the conditional mean-square error is:

$\displaystyle\text{MSE}(\mathbf{x})$ $\displaystyle=E\left[{\cal
E}^{2}|\mathbf{x}\right]$
$\displaystyle=\sum_{n=1}^{N}E\left[\left|\widehat{x}_{n}-x_{n}\right|^{2}|\mathbf{x}\right]$
$\displaystyle=\sum_{n=1}^{N}E\left[\left(\widehat{r}_{n}-r_{n}\right)^{2}|\mathbf{x}\right]+\sum_{n=1}^{N}E\left[\left(\widehat{i}_{n}-i_{n}\right)^{2}|\mathbf{x}\right].$
(7)

Whereas this mean-square error is, in general, different for each estimator,
we can use the Fisher information and the associated Cramér-Rao lower
bound[32] to place a lower bound on the mean-square error for any unbiased
estimator.

### Cramér-Rao Lower Bound

The Cramér-Rao lower bound on the mean-square error for the estimation of the
mode vector $\mathbf{x}=\mathbf{r}+j\mathbf{i}$ is:

$\text{MSE}(\mathbf{x})\geq\text{Trace}\left\\{\mathbf{J}^{-1}\right\\}$ (8)

where the matrix $\mathbf{J}$ is the Fisher information matrix for the
estimation of the parameters $\mathbf{r}$ and $\mathbf{i}$ based on the
statistical model for the measurements in Eq. (4). Because quantum-limited
measurements of the intensities are Poisson distributed, the Fisher
information matrix is [33]:

$\mathbf{J}=\left[\begin{array}[]{cc}\mathbf{J}_{r,r}&\mathbf{J}_{r,i}\\\
\mathbf{J}_{i,r}&\mathbf{J}_{i,i}\end{array}\right],$ (9)

where

$\mathbf{J}_{\beta,\zeta}=\sum_{m=1}^{M}\frac{\partial\left|y_{m}\right|^{2}}{\partial{\boldsymbol{\beta}}}\frac{\partial\left|y_{m}\right|^{2}}{\partial{\boldsymbol{\zeta}}^{T}}\frac{1}{\left|y_{m}\right|^{2}},$
(10)

with $\boldsymbol{\beta}$ and $\boldsymbol{\zeta}$ taking all possible
combinations of $\mathbf{r}$ and $\mathbf{i}$. Because

$\frac{\partial\left|y_{m}\right|^{2}}{\partial\mathbf{r}}=\mathbf{a}_{m}^{*}y_{m}^{*}+\mathbf{a}_{m}y_{m},$
(11)

and

$\frac{\partial\left|y_{m}\right|^{2}}{\partial\mathbf{i}}=j\mathbf{a}_{m}^{*}y_{m}^{*}-j\mathbf{a}_{m}y_{m},$
(12)

the blocks of the Fisher information are:

$\mathbf{J}_{r,r}=2\mathbf{I}+\mathbf{C}_{R},$ (13)
$\mathbf{J}_{i,i}=2\mathbf{I}-\mathbf{C}_{R},$ (14)

and

$\mathbf{J}_{r,i}=\mathbf{J}_{i,r}=\mathbf{C}_{I},$ (15)

where

$\displaystyle\mathbf{C}_{R}$
$\displaystyle=2\mbox{Re}\left\\{\sum_{m=1}^{M}\mathbf{a}_{m}\mathbf{a}_{m}^{T}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}\right\\},$
(16)
$\mathbf{C}_{I}=2\mbox{Im}\left\\{\sum_{m=1}^{M}\mathbf{a}_{m}\mathbf{a}_{m}^{T}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}\right\\},$
(17)

and the orthonormality of $\mathbf{A}$ results in:

$\displaystyle\sum_{m=1}^{M}\mathbf{a}_{m}\mathbf{a}_{m}^{\dagger}$
$\displaystyle=\mathbf{A}^{\dagger}\mathbf{A}$ $\displaystyle=\mathbf{I}.$
(18)

Accordingly, the Fisher information matrix is

$\mathbf{J}=2\mathbf{I}+\left[\begin{array}[]{cc}\mathbf{C}_{R}&\mathbf{C}_{I}\\\
\mathbf{C}_{I}&-\mathbf{C}_{R}\end{array}\right],$ (19)

and its inverse satisfies:

$\mbox{Trace}\left\\{\mathbf{J}^{-1}\right\\}\geq N,$ (20)

with equality when $\mathbf{C}_{R}=\mathbf{C}_{I}=\mathbf{0}$, or,
equivalently:

$\displaystyle\mathbf{C}$ $\displaystyle=\mathbf{C}_{R}+j\mathbf{C}_{I}$
$\displaystyle=\sum_{m=1}^{M}\mathbf{a}_{m}\mathbf{a}_{m}^{T}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}$
$\displaystyle=\mathbf{0}.$ (21)

To see this, we note that all the eigenvalues of the $2N\times 2N$ symmetric
matrix

$\mathbf{C}_{RI}=\left[\begin{array}[]{cc}\mathbf{C}_{R}&\mathbf{C}_{I}\\\
\mathbf{C}_{I}&-\mathbf{C}_{R}\end{array}\right],$ (22)

must be greater than $-2$ because the eigenvalues of the Fisher information
matrix $\mathbf{J}$ must be greater than $0$. Furthermore, if $\gamma$ is an
eigenvalue for $\mathbf{C}_{RI}$, then $-\gamma$ is also an eigenvalue because

$\left[\begin{array}[]{cc}\mathbf{C}_{R}&\mathbf{C}_{I}\\\
\mathbf{C}_{I}&-\mathbf{C}_{R}\end{array}\right]\left[\begin{array}[]{c}\mathbf{u}_{1}\\\
\mathbf{u}_{2}\end{array}\right]=\gamma\left[\begin{array}[]{c}\mathbf{u}_{1}\\\
\mathbf{u}_{2}\end{array}\right],$ (23)

implies that

$\left[\begin{array}[]{cc}\mathbf{C}_{R}&\mathbf{C}_{I}\\\
\mathbf{C}_{I}&-\mathbf{C}_{R}\end{array}\right]\left[\begin{array}[]{c}-\mathbf{u}_{2}\\\
\mathbf{u}_{1}\end{array}\right]=-\gamma\left[\begin{array}[]{c}-\mathbf{u}_{2}\\\
\mathbf{u}_{1}\end{array}\right].$ (24)

This means the $2N$ eigenvalues for the Fisher information are
$\lambda_{n}=2\pm\gamma_{n}$, and, consequently, the trace for the inverse of
this matrix is

$\displaystyle\mbox{Tr}\left\\{\mbox{CRLB}\right\\}$
$\displaystyle=\sum_{n=1}^{N}\left(\frac{1}{2+\gamma_{n}}+\frac{1}{2-\gamma_{n}}\right)$
$\displaystyle=\sum_{n=1}^{N}\frac{4}{4-\gamma_{n}^{2}}$ $\displaystyle\geq
N,$ (25)

with equality when all the eigenvalues for $\mathbf{C}_{RI}$ are equal to
zero, or when $\mathbf{C}_{R}=\mathbf{C}_{I}=\mathbf{0}$.

### Holographic Measurement

One way to accomplish this condition is to make the reference signal large
($\rho_{m}\gg\max\\{\mathbf{x}\\}$) and select the matrix $\mathbf{A}$ so that
$\mathbf{A}^{T}\mathbf{A}=\mathbf{0}$. Suppose, for example, the matrix
$\mathbf{A}$ directs each mode to four measurements, and each measurement
contains only one mode. That is, $\mathbf{A}$ is a $4N\times N$ matrix:

$\mathbf{A}=\left[\begin{array}[]{cccc}\mathbf{a}&&&\\\ &\mathbf{a}&&\\\
&&\ddots&\\\ &&&\mathbf{a}\end{array}\right],$ (26)

with

$\mathbf{a}=\frac{1}{\sqrt{4}}\left[\begin{array}[]{c}1\\\ j\\\ -1\\\
-j\end{array}\right].$ (27)

This corresponds to phase quadrature holography [34, 35], and, when the
reference amplitudes $\\{\rho_{m}\\}$ are sufficiently large, we have

$\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}\rightarrow 1,$ (28)

and

$\displaystyle\sum_{m=1}^{M}\mathbf{a}_{m}\mathbf{a}_{m}^{T}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}$
$\displaystyle\rightarrow\mathbf{A}^{T}\mathbf{A}$ $\displaystyle=\mathbf{0},$
(29)

so the system attains the quantum limit for field estimation.

### Phase Retrieval

When we use a reference signal for the field estimation, all the phases for
the modes $\mathbf{x}$ are defined relative to the phase for the reference.
Without a reference signal, we typically use the phase for one of the modes as
our reference phase. To accomplish this, we can arbitrarily define the
imaginary part for some mode to be zero, and reduce the estimation from $2N$
parameters to $2N-1$ parameters.

Even without the benefit of a reference signal, we can still get arbitrarily
close to the quantum-limit for the estimation error when we have a large
number of modes $N$. This result is achieved by using a group-testing
strategy. In contrast with holography—where we take 4 measurements of each
mode—for phase retrieval we take $Q$ measurements of each group of $L\geq 2$
modes. To demonstrate that phase retrieval measurements can induce the optimal
bound, we assume we make measurements on every such group. Because the total
number of groups is $T=\binom{N}{L}$, the overall measurement matrix takes the
form

$\mathbf{A}=\frac{1}{\sqrt{K}}\left[\begin{array}[]{c}\mathbf{W}\mathbf{S}_{1}\\\
\mathbf{W}\mathbf{S}_{2}\\\ \vdots\\\
\mathbf{W}\mathbf{S}_{T}\end{array}\right]$ (30)

where $\mathbf{W}$ is a $Q\times L$ matrix with $Q\geq L$. We ensure that
$\mathbf{W}$, and by extension $\mathbf{A}$, are orthonormal by setting

$W_{q,l}=\frac{1}{\sqrt{Q}}e^{j\frac{2\pi}{Q}ql},$ (31)

and $\mathbf{S}_{t}$ as an $L\times N$ binary matrix that selects $L$-element
combinations from the mode vector $\mathbf{x}$. For example, if $N=4$ and
$L=2$, the matrix $\mathbf{S}$ that selects modes $1$ and $3$ is

$\mathbf{S}=\left[\begin{array}[]{cccc}1&0&0&0\\\ 0&0&1&0\\\
\end{array}\right]$ (32)

The factor $1/\sqrt{K}$ accounts for the fact that each mode splits its
intensity across $K=\binom{N-1}{L-1}$ of the $T$ combinations. Whereas each
mode interferes with itself in $QK$ of these measurements, it interferes with
each of the other modes in $QP$ of the measurements, where
$P=\binom{N-2}{L-2}$ so the ratio of self- to cross-interference is

$\frac{P}{K}=\frac{L-1}{N-1}.$ (33)

Now, because the elements of $\mathbf{A}$ are such that

$|a_{m,n}|=\frac{1}{\sqrt{QK}},$ (34)

the magnitude for

$C_{n,r}=\sum_{m=1}^{M}a_{m,n}a_{m,r}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}},$
(35)

satisfies:

$\displaystyle\left|C_{n,r}\right|$
$\displaystyle=\left|\sum_{m=1}^{M}a_{m,n}a_{m,r}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}\right|$
$\displaystyle\leq\sum_{m=1}^{M}\left|a_{m,n}\right|\left|a_{m,r}\right|$
$\displaystyle=\frac{1}{QK}\sum_{m=1}^{M}\alpha_{m,n}\alpha_{m,r}$
$\displaystyle=\left\\{\begin{array}[]{cc}1&n=r\\\
\displaystyle\frac{L-1}{N-1}&n\neq r\end{array}\right.,$ (38)

where $\alpha_{m,n}$ is equal to $1$ if $x_{n}$ contributes to the $m$th
measurement, and equal to $0$ otherwise. For fixed $L$, then, the matrix
$\mathbf{C}$ becomes diagonal as $N\rightarrow\infty$, and the trace of the
$(2N-1)\times(2N-1)$ Cramér-Rao bound matrix is

$\mbox{Trace}\left\\{\mathbf{J}^{-1}\right\\}=\frac{1}{2+u_{1}}+\sum_{n=2}^{N}\frac{4}{4-u_{n}^{2}-v_{n}^{2}},$
(39)

where

$u_{n}=2\mbox{Re}\left\\{C_{n,n}\right\\},$ (40)
$v_{n}=2\mbox{Im}\left\\{C_{n,n}\right\\},$ (41)

and, to account for the overall phase, we’ve assumed that the imaginary part
for $x_{1}$ is zero.

To analyze the performance bound, then, we can consider

$\displaystyle|C_{n,n}|^{2}$
$\displaystyle=\sum_{m=1}^{M}\sum_{p=1}^{M}a_{m,n}^{2}(a_{p,n}^{*})^{2}\frac{y_{m}^{2}}{\left|y_{m}\right|^{2}}\frac{(y_{p}^{*})^{2}}{\left|y_{p}\right|^{2}}$
$\displaystyle=\sum_{m=1}^{M}\sum_{p=1}^{M}a_{m,n}^{2}(a_{p,n}^{*})^{2}e^{j2(\theta_{m}-\theta_{p})},$
(42)

where $\theta_{m}$ is the phase for $y_{m}$. $|C_{n,n}|^{2}$ will, in general,
depend on both the measurement matrix $\mathbf{A}$ and the values for the
modes $\mathbf{x}$. Suppose, though, we model the uncertainty about
$\mathbf{x}$ as a zero-mean random vector with an average intensity per mode
equal to $N_{p}$:

$E\left[\frac{1}{N}\sum_{n=1}^{N}|x_{n}|^{2}\right]=N_{p}.$ (43)

The maximum-entropy distribution for this constraint is a zero-mean, complex
Gaussian random vector with diagonal covariance

$E\left[\mathbf{x}\mathbf{x}^{\dagger}\right]=N_{p}\mathbf{I},$ (44)

and, in this situation, $\mathbf{y}$ is a zero-mean, complex Gaussian random
vector with the covariance:

$\displaystyle E\left[\mathbf{y}\mathbf{y}^{\dagger}\right]$
$\displaystyle=\mathbf{R}$
$\displaystyle=N_{p}\mathbf{A}\mathbf{A}^{\dagger}.$ (45)

Then, if we denote

$e^{j2\theta_{m}}=\frac{y_{m}^{2}}{|y_{m}|^{2}},$ (46)

this phase signal will be zero-mean with the covariance [36]:

$\displaystyle E[e^{j2(\theta_{m}-\theta_{p})}]$ $\displaystyle=\Gamma_{m,p}$
$\displaystyle=e^{j2\phi_{m,p}}\left[1+\frac{(1-\rho_{m,p}^{2})\ln(1-\rho_{m,p}^{2})}{\rho_{m,p}^{2}}\right],$
(47)

where $\phi_{m,p}$ is the phase and $\rho_{m,p}^{2}$ is the amplitude of

$\frac{R^{2}_{m,p}}{R_{m,m}R_{p,p}}.$

Because

$\left[1+\frac{(1-\rho_{m,p}^{2})\ln(1-\rho_{m,p}^{2})}{\rho_{m,p}^{2}}\right]\leq\rho_{m,p}^{2},$
(48)

we have that

$\displaystyle E\left[|C_{n,n}|^{2}\right]$
$\displaystyle=\sum_{m=1}^{M}\sum_{p=1}^{M}a_{m,n}^{2}(a_{p,n}^{*})^{2}\Gamma_{m,p}$
$\displaystyle\leq\frac{1}{(QK)^{2}}\sum_{m=1}^{M}\sum_{p=1}^{M}\alpha_{m,n}\alpha_{p,n}\rho_{m,p}^{2}.$
(49)

We’ve utilized the inequality here to form an upper bound on the lower bound.
If, instead, we wanted to analyze the bound for a particular value of $L$, we
could use the exact form for this expectation.

Now, the form for our measurement matrix $\mathbf{A}$ is such that

$\displaystyle\mathbf{R}$ $\displaystyle=\mathbf{A}\mathbf{A}^{\dagger}$
$\displaystyle=\frac{1}{K}\left[\begin{array}[]{cccc}\mathbf{W}\mathbf{W}^{\dagger}&\mathbf{W}\mathbf{S}_{1}\mathbf{S}_{2}^{T}\mathbf{W}^{\dagger}&\cdots&\mathbf{W}\mathbf{S}_{1}\mathbf{S}_{T}^{T}\mathbf{W}^{\dagger}\\\
\mathbf{W}\mathbf{S}_{2}\mathbf{S}_{1}^{T}\mathbf{W}^{\dagger}&\mathbf{W}\mathbf{W}^{\dagger}&\cdots&\mathbf{W}\mathbf{S}_{2}\mathbf{S}_{T}^{T}\mathbf{W}^{\dagger}\\\
\vdots&&\ddots&\vdots\\\
\mathbf{W}\mathbf{S}_{T}\mathbf{S}_{1}^{T}\mathbf{W}^{\dagger}&\mathbf{W}\mathbf{S}_{T}\mathbf{S}_{2}^{T}\mathbf{W}^{\dagger}&\cdots&\mathbf{W}\mathbf{W}^{\dagger}\end{array}\right],$
(54)

where $\mathbf{S}_{t}\mathbf{S}_{t}^{T}=\mathbf{I}$, and

$\displaystyle\left[\mathbf{W}\mathbf{W}^{\dagger}\right]_{q,r}$
$\displaystyle=\frac{1}{Q}\sum_{l=1}^{L-1}e^{j\frac{2\pi}{Q}l(q-r)}$
$\displaystyle=\frac{1}{Q}e^{j\frac{\pi(L-1)}{Q}(q-r)}\frac{\sin\pi\frac{L}{Q}(q-r)}{\sin\pi\frac{1}{Q}(q-r)},$
(55)

so that the diagonal elements of $\mathbf{R}$ are

$R_{m,m}=\frac{L}{KQ}.$ (56)

To evaluate the upper bound for $E[|C_{n,n}|^{2}]$, we need to evaluate
$\rho_{m,p}^{2}$ for the $K^{2}$ blocks of $\mathbf{R}$ that have a
contribution from $x_{n}$, which includes $K$ diagonal blocks and $K(K-1)$
off-diagonal blocks. To do this, we’ll need to normalize each block by
$(\frac{L}{KQ})^{2}$, then sum the squared magnitudes for those normalized
values.

To better understand each of the off-diagonal blocks for some value for $n$,
suppose that $n=2$, $L=5$, and the combination corresponding to
$\mathbf{S}_{1}$ contains $x_{1},x_{2},x_{3},x_{8}$ and $x_{9}$, and the
combination corresponding to $\mathbf{S}_{5}$ contains
$x_{2},x_{4},x_{7},x_{9}$ and $x_{10}$. The two combinations share the modes
$x_{2}$ and $x_{9}$, so that

$\displaystyle\frac{1}{K}\left[\mathbf{W}\mathbf{S}_{1}\mathbf{S}_{5}^{T}\mathbf{W}^{\dagger}\right]_{q,r}$
$\displaystyle=\frac{1}{K}\left[\mathbf{w}_{2}\mathbf{w}_{1}^{\dagger}+\mathbf{w}_{5}\mathbf{w}_{4}^{\dagger}\right]_{q,r}$
$\displaystyle=\frac{1}{QK}\left[e^{j\frac{2\pi}{Q}2q}e^{-j\frac{2\pi}{Q}1r}+e^{j\frac{2\pi}{Q}5q}e^{-j\frac{2\pi}{Q}4r}\right],$
(57)

where $\mathbf{w}_{l}$ is the $l$th column of $\mathbf{W}$. We get this result
because $x_{2}$ is in the second position in the first combination, first
position in the second combination; and $x_{9}$ is in the fifth position in
the first combination, fourth position in the second combination. Now, because

$\displaystyle\left|\frac{1}{K}\left[\mathbf{W}\mathbf{S}_{1}\mathbf{S}_{5}^{T}\mathbf{W}^{\dagger}\right]_{q,r}\right|^{2}$
$\displaystyle=\left(\frac{1}{QK}\right)^{2}\left|e^{j\frac{2\pi}{Q}2q}e^{-j\frac{2\pi}{Q}1r}+e^{j\frac{2\pi}{Q}5q}e^{-j\frac{2\pi}{Q}4r}\right|^{2}$
$\displaystyle\leq\left(\frac{1}{QK}\right)^{2}2^{2},$ (58)

we can infer the general result that

$\left|\frac{1}{K}\left[\mathbf{W}\mathbf{S}_{t}\mathbf{S}_{\tau}^{T}\mathbf{W}^{\dagger}\right]_{q,r}\right|^{2}\leq\left(\frac{1}{QK}\right)^{2}L_{t,\tau}^{2},$
(59)

where $L_{t,\tau}$ is the number of modes in common between the $t$th and
$\tau$th combinations. If we divide by $(\frac{L}{KQ})^{2}$ to obtain the
corresponding element of $\rho^{2}$, then sum over each block, we obtain the
result that the contribution from each block that contains $x_{n}$ to the sum
that defines $E[|C_{n,n}|^{2}]$ is less than or equal to
$Q^{2}L^{2}_{t,\tau}/L^{2}$. The total contribution for the $K^{2}$ blocks is
then

$\displaystyle\sum_{m=1}^{M}\sum_{p=1}^{M}\alpha_{m,n}\alpha_{p,n}\rho^{2}_{m,p}\leq\frac{Q^{2}}{L^{2}}\sum_{(t,\tau)\in{\cal
T}_{n}}L^{2}_{t,\tau},$ (60)

where ${\cal T}_{n}$ is the $K^{2}$ pairs of $L$-element combinations
$(t,\tau)$ that include $x_{n}$.

The number of modes in common between the $t$th and $\tau$th combinations,
$L_{t,\tau}$, must be less than or equal to $L$ and greater than or equal to
$1$. When $N$ is much greater than $L$, though, nearly all of the values for
$L_{t,\tau}$ are equal to one. This follows from the fact that the probability
of no common elements in any two randomly-selected combinations of $L-1$
elements from $N-1$ is equal to

$\frac{\displaystyle\binom{N-L}{L-1}}{\displaystyle\binom{N-1}{L-1}}\approx\frac{\left(1-\frac{L}{N}\right)^{2N-2L}}{\left(1-\frac{2L}{N}\right)^{N-2L}}\xrightarrow{N\gg
L}1.$ (61)

where we use Stirling’s approximation to expand the factorials. In the limit
of $N\gg L$, then, we can place an upper bound on $E[|C_{n,n}|^{2}]$:

$\displaystyle E\left[|C_{n,n}|^{2}\right]$
$\displaystyle=\sum_{m=1}^{M}\sum_{p=1}^{M}a_{m,n}^{2}(a_{p,n}^{*})^{2}\Gamma_{m,p}$
$\displaystyle\leq\frac{1}{(QK)^{2}}\sum_{m=1}^{M}\sum_{p=1}^{M}\alpha_{m,n}\alpha_{p,n}\rho_{m,p}^{2}$
$\displaystyle\leq\frac{1}{L^{2}}.$ (62)

Therefore, by selecting $L$ large enough, a phase-retrieval system can get
arbitrarily close to the quantum-limit for field estimation in the limit of a
large number of modes. Of course when $N$ and $L$ are large, the number of
combinations $T$ is generally impractical, but this result shows that a phase-
retrieval system can encode all the information needed to attain the quantum-
limit and achieve a mean-square error equal to the total number of modes, or a
per mode error equal to 1.

## Example Measurement Systems

This section presents simulations of phase retrieval systems derived from the
strategy of the previous section. We consider cases from $L=2$ to $L=5$.
Recognizing that measurement of $\binom{N}{L}$ groups as suggested in the last
section is impractical, we instead consider randomly selected subsets of the
groups proposed above. Here we show that $Ln\log{}n$ measurements are
sufficient to approach the quantum limit and that such measurements allow
computationally tractable phase retrieval. As suggested in the previous
section, increasing $L$ reduces MSE.

For $L=2$, we assume random connections which require that 1) each mode in x
is connected pairwise to $p=2\log_{2}n$ randomly selected modes, and 2) there
does not exist a subset $\bf\tilde{x}\subset x$ such that no connections exist
between modes in $\bf\tilde{x}$ and modes in $\bf x\backslash\tilde{x}$. This
is easily understood as a connected graph with $n$ vertices, and each vertex
has $p$ adjacent vertices. Each group of modes is measured using $\mathbf{W}$
as defined in Eqn. 31 with $Q=3$. In total, we consider $n\log_{2}n$ pairs and
$3n\log_{2}n$ measurements. For $L=3,4,5$ we select $Q=L$. The connection
strategy for $L>2$ is the same as for $L=2$, each connection randomly draws
$L$ modes. For $n=2^{k}$, each mode is measured/selected $L\log n$ times. The
number of total measurements is $Ln\log n$.

We solve for $x$ by first-order gradient-based optimization using the Adam
optimizer implemented in TensorFlow [37, 38]. MSE and computation time as a
function of $n$ are summarized in Figure 1. In these simulations we
independently draw the value of each mode in the input field from complex
standard normal distribution $CN(0,1)$. Measurements are drawn from a Poisson
distribution with $\bar{K}=10^{4}$ per mode. We repeated the simulation
several times for each value of $n$ and averaged the results. As shown in the
figure and expected from the results of the previous section, MSE falls as $L$
and $n$ increase. The slight rise in MSE at high dimension for $L=5$ is
associated with quantization error in dividing $10^{4}$ photons into 100
measurements, increasing the number of photons per mode decreases the error.
Figure 1(b) shows the computation time per mode, i.e., total runtime divided
by the signal dimension. For the operating system and other hardware supports,
we use Ubuntu 18.04, Intel Core i7-7700 @ 3.60GHz and 16 GB of RAM. The
computer code used for this paper is available online [39].

(a) MSE

(b) TIME

Figure 1: The reconstruction error and the computing time per mode for
different $L$. The reconstruction error of a holographic system is also shown
for comparison. For $\log n<8$ the time to launch the kernel dominates the
total runtime.

$n$ is ultimately limited by computing resources. For example, to process a
$1024\times 1024$ image ($2^{20}$ modes) requires $>10^{7}$ measurements and
an optimization engine with at least 10GB RAM. As $n$ increases additional
challenges arise in interferometer design, physical implementation and
numerical estimation. One approach to scalable design may be to build the
interferometer network for a large field on top of blocks from smaller fields.
This system is illustrated in Figure 2. The field x is divided into several
blocks, called rank 0 blocks. Within each rank 0 block, random connections and
phase retrieval can be applied as described above. The two rank 0 blocks are
called connected when each mode in one block is randomly connected to one or
more modes in another block. The relative phase between the two connected
blocks can be computed from the cross-block connections. Figure 2(a) shows an
example of two connected rank 0 blocks: each mode is connected to 4 modes
(mode A, B and C) within the block and 2 modes across the block (mode A and
D).

(a) Connections between modes

(b) The multiscale connection strategy

Figure 2: Illustration of the multiscale system. (a) Connections between
modes. The bounding boxes represents the rank 0 blocks, the blue dots
represent the mode to measure, and the arrows represent the connections
between modes within (purple) and across (red) the block. (b) The multiscale
connection strategy. The labelled squares represent rank 0 blocks, the arrows
represent the connections between blocks, and the blocks are put together when
the relative phases between them are retrieved.

We find the global phase by successively building higher rank blocks. This
process is illustrated in Figure 2(b) with 8 rank 0 blocks. First we build
rank 1 blocks by connecting pairs of neighboring rank 0 blocks and applying
proper phase shift to one of the block. Then each rank 1 block consists of 2
rank 0 blocks. Next we build rank 1 blocks by pairing rank 1 blocks. Each rank
0 block is connected to the corresponding rank 0 block in the other rank 1
block, and the relative phase between the the rank 1 blocks can be estimated
from the all the cross-block connections. For example, in Figure 2(b), block 1
is connected to block 3, block 2 is connected to block 4, and the relative
phase between the two rank 1 blocks can be estimated by averaging the relative
phases between the rank 0 blocks. After phase shift, each rank 2 block
consists of 4 rank 0 blocks. This process continues until all the relative
phases are recovered.

We test this multiscale system with simulations. We assume $L=2$, $n=2^{k}$
and $p=2\log n$. There are $2^{q}$ modes in a rank 0 block, each connected to
another $2q$ modes in the same block, and each mode is connected to 2 more
modes in another rank 0 block when a higher rank block is to be built. Because
the cross-block connections are only used to find the relative phases between
blocks, only a short exposure time is necessary. In the simulation, only five
percent of the signal energy is used for cross-block connections. The
reconstruction error is shown in Figure 3. For $q=12$ reconstruction time is
dominated by the intrablock phase retrieval, scaling to $n=2^{20}$ requires
phase retrieval over 256 blocks of 4096 elements. The processing time per
pixel remains approximately equal to the minimum associated with $n=2^{12}$ at
the expense of a 0.7 dB increase in MSE.

Figure 3: Simulation results on the multiscale system. Four different rank 0
blocks ($q=5,8,10,12$) are implemented, along with the global optimization and
holography. In the global optimization method, when $n>12$, the random mapping
is scaled up from $n=12$.

We have attempted to compare results presented here with previous phase
retrieval studies, which have shown that $m\geq 4n-4$ measurements are
sufficient to uniquely recover $x$ [40, 41, 42]. The number of measurements
made under our approach, $Ln\log n$, is much larger than in these studies, but
the computational complexity of phase retrieval under our approach is much
lower than previous work using dense random matrices $A$, which require
$O(n^{2})$ operations per forward projection rather than the $O(n\log n)$
required here. Our computational complexity is similar to Fourier
measurements, but simple Fourier transform sampling does not robustly achieve
quantum limited phase retrieval. Unfortunately, direct comparison between
results is not possible because (1) previous methods fail to converge for a
significant fraction of random inputs, rendering calculation of mean MSE
impossible, and (2) higher computational complexity renders previous methods
are impractical for $n>2^{12}$. In contrast, the proposed method here always
recovers the signal. Because the proposed method uses sparse connections,
built-in operators for sparse matrices in TensorFlow, or other platforms, can
be used to speed computing and decrease memory, allowing for the
reconstruction of megapixel fields.

## Discussion

This paper considers optimization of $\mathbf{A}$ for estimation of
$\mathbf{x}$ from the intensities of $\mathbf{y}=\mathbf{A}\mathbf{x}$. In the
presence of Poisson noise we prove that there exist $\mathbf{A}$ such that
$\mathbf{x}$ can be estimated with quantum-limited precision. We further
demonstrate in simulation that there exist practical measurement strategies
that return $\mathbf{x}$ with precision near this limit. The implication of
these results is that phase retrieval systems can match the precision of
holographic measurement systems. We have not proven, but we here speculate
that phase retrieval may be capable of exceeding holographic performance when
other system parameters are considered. In particular, holography must assume
that the reference and signal are mutually coherent whereas for phase
retrieval the signal need only be self-coherent.

We recognize that the interferometers proposed here will not be easily
constructed, but we suggest that their construction forms a suitable challenge
for emerging integrated photonics technologies. In practice, of course, the
use of programmable interconnect networks in such devices could lead to
measurement matrices significantly more advanced than those proposed here.
Beyond issues of physical implementation, the results presented here raise
more technological and theoretical questions than they answer. For example,
the performance of the proposed system in analyzing partially coherent fields
and natural image data remains to be explored.

## References

  * [1] D Gabor.  A new microscopic principle.  Nature, 161(4098):777, 1948.
  * [2] Ulf Schnars, Claas Falldorf, John Watson, and Werner Jüptner.  Digital holography.  In Digital Holography and Wavefront Sensing, pages 39–68. Springer, 2015.
  * [3] RW Gerehberg and WO Saxton.  A practical algorithm for the determination of phase from image and diffraction plane pictures.  Optik, 35(2):237–246, 1972.
  * [4] RA Gonsalves.  Phase retrieval from modulus data.  JOSA, 66(9):961–964, 1976.
  * [5] James R Fienup.  Phase retrieval algorithms: a comparison.  Applied optics, 21(15):2758–2769, 1982.
  * [6] Yuxin Chen and Emmanuel J Candès.  Solving random quadratic systems of equations is nearly as easy as solving linear systems.  Communications on Pure and Applied Mathematics, 70(5):822–883, 2017\.
  * [7] William H Southwell.  Wave-front estimation from wave-front slope measurements.  JOSA, 70(8):998–1006, 1980.
  * [8] M. P. Rimmer.  Method for evaluating lateral shearing interferograms.  Appl. Opt., 13(3):623–629, Mar 1974.
  * [9] Carlton M Caves.  Quantum-mechanical noise in an interferometer.  Physical Review D, 23(8):1693, 1981.
  * [10] Yoav Shechtman, Yonina C Eldar, Oren Cohen, Henry Nicholas Chapman, Jianwei Miao, and Mordechai Segev.  Phase retrieval with application to optical imaging: a contemporary overview.  IEEE signal processing magazine, 32(3):87–109, 2015.
  * [11] Congli Wang, Xiong Dun, Qiang Fu, and Wolfgang Heidrich.  Ultra-high resolution coded wavefront sensor.  Optics express, 25(12):13736–13746, 2017.
  * [12] Ryoichi Horisaki, Yusuke Ogura, Masahiko Aino, and Jun Tanida.  Single-shot phase imaging with a coded aperture.  Optics letters, 39(22):6466–6469, 2014.
  * [13] Emmanuel J Candes, Xiaodong Li, and Mahdi Soltanolkotabi.  Phase retrieval from coded diffraction patterns.  Applied and Computational Harmonic Analysis, 39(2):277–299, 2015\.
  * [14] JM Rodenburg, AC Hurst, AG Cullis, BR Dobson, F Pfeiffer, O Bunk, C David, K Jefimovs, and I Johnson.  Hard-x-ray lensless imaging of extended objects.  Physical review letters, 98(3):034801, 2007.
  * [15] Gang Wang, Georgios B Giannakis, and Yonina C Eldar.  Solving systems of random quadratic equations via truncated amplitude flow.  IEEE Transactions on Information Theory, 64(2):773–794, 2017.
  * [16] Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch, Christoph Studer, and Tom Goldstein.  Phasepack: A phase retrieval library.  Asilomar Conference on Signals, Systems, and Computers, 2017.
  * [17] Martin Harwit and Neil Sloane.  Hadamard transform optics.  Academic Press, 1979.
  * [18] Edward E Fenimore and Thomas M Cannon.  Coded aperture imaging with uniformly redundant arrays.  Applied optics, 17(3):337–347, 1978.
  * [19] M Strojnik, G Paez, and M Mantravadi.  Lateral shear interferometers.  Optical Shop Testing, pages 122–184, 2007.
  * [20] D Malacara-Doblado and I Ghozeil.  Hartmann, hartmann-shack, and other screen tests.  Optical Shop Testing, 361397, 2007.
  * [21] Luc Gilles, Curtis R Vogel, and Brent L Ellerbroek.  Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction.  JOSA A, 19(9):1817–1822, 2002.
  * [22] Amos Talmi and Erez N Ribak.  Wavefront reconstruction from its gradients.  JOSA A, 23(2):288–297, 2006.
  * [23] Iacopo Mochi and Kenneth A Goldberg.  Modal wavefront reconstruction from its gradient.  Applied Optics, 54(12):3780–3785, 2015.
  * [24] Remy Tumbar, Ronald A Stack, and David J Brady.  Wave-front sensing with a sampling field sensor.  Applied optics, 39(1):72–84, 2000.
  * [25] Luca Pezzè, Philipp Hyllus, and Augusto Smerzi.  Phase-sensitivity bounds for two-mode interferometers.  Physical Review A, 91(3):032103, 2015.
  * [26] Emmett N Leith and Juris Upatnieks.  Reconstructed wavefronts and communication theory.  JOSA, 52(10):1123–1130, 1962.
  * [27] Jacques Carolan, Christopher Harrold, Chris Sparrow, Enrique Martín-López, Nicholas J. Russell, Joshua W. Silverstone, Peter J. Shadbolt, Nobuyuki Matsuda, Manabu Oguma, Mikitaka Itoh, Graham D. Marshall, Mark G. Thompson, Jonathan C. F. Matthews, Toshikazu Hashimoto, Jeremy L. O’Brien, and Anthony Laing.  Universal linear optics.  Science, 349(6249):711–716, 2015.
  * [28] William R. Clements, Peter C. Humphreys, Benjamin J. Metcalf, W. Steven Kolthammer, and Ian A. Walmsley.  Optimal design for universal multiport interferometers.  Optica, 3(12):1460–1465, Dec 2016.
  * [29] Tiehui Su, Guangyao Liu, Katherine E Badham, Samuel T Thurman, Richard L Kendrick, Alan Duncan, Danielle Wuchenich, Chad Ogden, Guy Chriqui, Shaoqi Feng, et al.  Interferometric imaging using si 3 n 4 photonic integrated circuits for a spider imager.  Optics express, 26(10):12801–12812, 2018.
  * [30] Joseph C. Marron and Kirk S. Schroeder.  Three-dimensional lensless imaging using laser frequency diversity.  Appl. Opt., 31(2):255–262, Jan 1992.
  * [31] Jason W. Stafford, Bradley D. Duncan, and David J. Rabb.  Phase gradient algorithm method for three-dimensional holographic ladar imaging.  Appl. Opt., 55(17):4611–4620, 2016.
  * [32] Steven M. Kay.  Fundamentals of Statistical Signal Processing, Volume 1: Estimation Theory.  Prentice Hall, 1993.
  * [33] Donald L. Snyder and Michael I. Miller.  Random point processes in time and space.  Springer-Verlag New York, 1991.
  * [34] Songcan Lai, Brian King, and Mark A Neifeld.  Wave front reconstruction by means of phase-shifting digital in-line holography.  Optics communications, 173(1-6):155–160, 2000.
  * [35] Jung-Ping Liu, Ting-Chung Poon, Gui-Syu Jhou, and Po-Jung Chen.  Comparison of two-, three-, and four-exposure quadrature phase-shifting holography.  Appl. Opt., 50(16):2443–2450, Jun 2011.
  * [36] Kenneth S. Miller.  Complex stochastic processes.  Addison-Wesley, 1974.
  * [37] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [38] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [39] Wang, C. and Brady, D. J. and Schulz, T. J. .  Multiscale phase retrieval computer code, 2020.  https://github.com/djbrady/multiscalePhaseRetrieval/.
  * [40] Gang Wang, Georgios B Giannakis, Yousef Saad, and Jie Chen.  Solving almost all systems of random quadratic equations.  arXiv preprint arXiv:1705.10407, 2017.
  * [41] Radu Balan, Pete Casazza, and Dan Edidin.  On signal reconstruction without phase.  Applied and Computational Harmonic Analysis, 20(3):345–356, 2006\.
  * [42] Aldo Conca, Dan Edidin, Milena Hering, and Cynthia Vinzant.  An algebraic characterization of injectivity in phase retrieval.  Applied and Computational Harmonic Analysis, 38(2):346–356, 2015\.


isotropy assumption is reduced to the moment condition
${\mathbb{E}}X_{i}^{2}=1$ ($1\leq i\leq n$).

In particular, it follows from (16.5) that

${\mathbb{E}}_{\theta}\,\rho(F_{\theta},\Phi)\leq{\mathbb{E}}_{\theta}\,\|F_{\theta}-\Phi\|_{\rm
TV}\leq\frac{c}{n}.$ (16.6)

##  17\. Improved Rates Under Correlation-Type Conditions

Up to a logarithmically growing term, the improved rate as in the upper bound
(16.3) can be achieved under more flexible correlation-type conditions (in
comparison with independence). For example, one may consider an optimal value
$\Lambda=\Lambda(X)$ in the relation

${\rm
Var}\bigg{(}\sum_{i,j=1}^{n}a_{ij}X_{i}X_{j}\bigg{)}\leq\Lambda\sum_{i,j=1}^{n}a_{ij}^{2}\qquad(a_{ij}\in{\mathbb{R}}),$
(17.1)

which we call that the random vector $X=(X_{1},\dots,X_{n})$ satisfies a
second order correlation condition with constant $\Lambda$. This quantity is
finite as long as the moment ${\mathbb{E}}\,|X|^{4}$ is finite.

To relate $\Lambda$ to the moment-type characteristics which we discussed
before, one may apply (17.1) with $a_{ij}=\delta_{ij}$ or (as another option)
with $a_{ij}=\theta_{i}\theta_{j}$,
$\theta=(\theta_{1},\dots,\theta_{n})\in{\mathbb{S}}^{n-1}$. This gives that

$\sigma_{4}^{2}\leq\Lambda,\quad
m_{4}^{2}\leq\sup_{\theta\in{\mathbb{S}}^{n-1}}{\mathbb{E}}S_{\theta}^{4}\leq
1+\Lambda,$

where in the last inequality we should assume that
${\mathbb{E}}S_{\theta}^{2}=1$ for all $\theta$ (i.e. $X$ is isotropic). In
the latter case, necessarily $\Lambda\geq\frac{n-1}{n}$, so that $\Lambda$ is
bounded away from zero.

If the distribution of $X$ is “regular” in some sense, one may also bound
$\Lambda$ from above. For example, this is the case when it shares a Poincaré-
type inequality

$\lambda_{1}{\rm Var}(u(X))\leq{\mathbb{E}}\,|\nabla u(X)|^{2},$ (17.2)

which is required to hold in the class of all bounded, smooth functions $u$ on
${\mathbb{R}}^{n}$ with a constant $\lambda_{1}>0$ independent of $u$ (called
the spectral gap). We then have

$\Lambda\leq\frac{4}{\lambda_{1}^{2}},\quad\Lambda\leq\frac{4}{\lambda_{1}},$
(17.3)

where in the second inequality we assume that $X$ is isotropic.

The following relation is established in [10].

Theorem 17.1. If the distribution of $X$ is isotropic and symmetric about the
origin, then

${\mathbb{E}}_{\theta}\,\rho(F_{\theta},\Phi)\leq c\Lambda\,\frac{\log n}{n}.$
(17.4)

The proof is based on the second order spherical concentration phenomenon
which was developed in [7] with the aim of applications to randomized central
limit theorems. It indicates that the deviations of any smooth function
$u(\theta)$ on ${\mathbb{S}}^{n-1}$ from the mean
${\mathbb{E}}_{\theta}u(\theta)$ are at most of the order $1/n$, provided that
$u$ is orthogonal in $L^{2}({\mathbb{R}}^{n},{\mathfrak{s}}_{n-1})$ to all
linear functions and has a “bounded” Hessian (the matrix of second order
partial derivatives). Being applied to the characteristic functions
$u(\theta)=f_{\theta}(t)$, this property yields an upper bound

${\mathbb{E}}_{\theta}\,|f_{\theta}(t)-f(t)|^{2}\,\leq\,\frac{c\Lambda
t^{4}}{n^{2}}$

on every interval $|t|\leq An^{1/5}$ with constants $c>0$ depending on the
parameter $A\geq 1$ only. This estimate can be used to bound the integrals in
(8.4) to get a similar variant of (17.4).

The symmetry hypothesis in Proposition 17.1 may be dropped, if $\Lambda$ is
replaced by $\lambda_{1}^{-1}$ which is a larger quantity according to (17.3).
In addition, one can control large deviations of the distance
$\rho(F_{\theta},\Phi)$ for most of the directions $\theta$ (rather than on
average). The corresponding assertions are obtained in [11].

Theorem 17.2. Let $X$ be an isotropic random vector in ${\mathbb{R}}^{n}$ with
mean zero and a positive Poincaré constant $\lambda_{1}$. Then

${\mathbb{E}}_{\theta}\,\rho(F_{\theta},\Phi)\,\leq\,c\lambda_{1}^{-1}\,\frac{\log
n}{n}.$ (17.5)

Moreover, for all $r>0$,

${\mathfrak{s}}_{n-1}\Big{\\{}\rho(F_{\theta},\Phi)\geq
c\lambda_{1}^{-1}\frac{\log n}{n}\,r\Big{\\}}\,\leq\,2\,e^{-\sqrt{r}}.$

The logarithmic term in (17.5) may be removed using the less sensitive
$L^{2}$-distance:

${\mathbb{E}}_{\theta}\,\omega^{2}(F_{\theta},\Phi)\leq\frac{c}{\lambda_{1}^{2}\,n^{2}}.$

There is an extensive literature devoted to bounding the spectral gap
$\lambda_{1}$ from below. In particular, it is positive for any log-concave
probability distribution on ${\mathbb{R}}^{n}$. A well-known conjecture raised
by Kannan, Lovász and Simonovits asserts that $\lambda_{1}$ is actually
bounded away from zero, as long as the random vector $X$ has an isotropic log-
concave distribution (cf. [18]). The best known dimensional lower bound up to
date is due to Klartag and Lehec [22] who showed that

$\lambda_{1}\geq\frac{c}{(\log n)^{\alpha}}$

for some absolute positive constants $c$ and $\alpha$ (one may take
$\alpha=10$). Applying this bound in Theorem 17.2, we therefore obtain:

Corollary 17.3. Let $X$ be an isotropic random vector in ${\mathbb{R}}^{n}$
with mean zero and a log-concave probability distribution. Then with some
absolute positive constants $c$ and $\alpha$

${\mathbb{E}}_{\theta}\,\rho(F_{\theta},\Phi)\,\leq\,\frac{c(\log
n)^{\alpha}}{n}.$ (17.6)

Thus, there is a certain extension of Klartag’s bound (16.6) at the expense of
a logarithmic factor to the entire class of isotropic log-concave probability
distributions on ${\mathbb{R}}^{n}$.

One may also argue in the opposite direction: upper bounds of the form

${\mathbb{E}}_{\theta}\,\rho(F_{\theta},\Phi)\leq\frac{c(\log
n)^{\beta}}{n},\quad\beta>0,$

in the class of log-concave probability distributions on ${\mathbb{R}}^{n}$
imply lower bounds $\lambda_{1}\geq c\,(\log n)^{-\beta^{\prime}}$ with some
$\beta^{\prime}>0$, cf. [10].

## References

  * [1]
  * [2] M. Anttila, K. Ball, and I. Perissinaki. The central limit problem for convex bodies. Trans. Amer. Math. Soc. 355 (2003), no. 12, 4723–4735.
  * [3] H. Bateman. Higher transcendental functions, Vol. II. McGraw-Hill Book Company, Inc., 1953, 396 pp.
  * [4] S. G. Bobkov. On concentration of distributions of random weighted sums. Ann. Probab. 31 (2003), no. 1, 195–215.
  * [5] S. G. Bobkov. Closeness of probability distributions in terms of Fourier-Stieltjes transforms. Russian Math. Surveys, vol. 71, issue 6, (2016), 1021–1079. Translated from: Uspekhi Matem. Nauk, vol. 71, issue 6 (432), (2016), 37–98.
  * [6] S. G. Bobkov. Edgeworth corrections in randomized central limit theorems. Geometric Aspects of Functional Analysis, 2256 (2020), pp. 71–97.
  * [7] S. G. Bobkov, G. P. Chistyakov, and F. Götze. Second-order concentration on the sphere. Commun. Contemp. Math. 19 (2017), no. 5, 1650058, 20 pp.
  * [8] S. G. Bobkov, G. P. Chistyakov, and F. Götze. Gaussian mixtures and normal approximation for V. N. Sudakov’s typical distributions. Zap. Nauchn. Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov. (POMI) 457 (2017), Veroyatnost i Statistika. 25, 37–52; reprinted in J. Math. Sci. (N.Y.) 238 (2019), no. 4, 366–376.
  * [9] S. G. Bobkov, G. P. Chistyakov, and F. Götze. Berry-Esseen bounds for typical weighted sums. J. Electron. Probab. 23 (2018), no. 92, 1–22.
  * [10] S. G. Bobkov, G. P. Chistyakov, and F. Götze. Normal approximation for weighted sums under a second order correlation condition. Ann. Probab. 48 (2020), no. 3, 1202–1219.
  * [11] S. G. Bobkov, G. P. Chistyakov, and F. Götze. Poincaré-type inequalities and normal approximation for weighted sums. Electron. J. Probab. 25 (2020), Paper No. 155, 31 pp.
  * [12] S. G. Bobkov, G. P. Chistyakov, and F. Götze. Concentration and Gaussian approximation for randomized sums. Monograph (in preparation).
  * [13] S. G. Bobkov, and F. Götze. Concentration inequalities and limit theorems for randomized sums. Probab. Theory Related Fields 137 (2007), no. 1–2, 49–81.
  * [14] S. Brazitikos, A. Giannopoulos, P. Valettas, and B-H. Vritsiou. Geometry of isotropic convex bodies. Mathematical Surveys and Monographs, 196. American Mathematical Society, Providence, RI, 2014. xx+594 pp.
  * [15] A. Dvoretzky, J. Kiefer, and J. Wolfowitz. Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator. Ann. Math. Statist. 27 (1956), 642–669.
  * [16] Goldstein, L.; Reinert, G. Stein’s method and the zero bias transformation with application to simple random sampling. Ann. Appl. Probab. 7 (1997), no. 4, 935–952.
  * [17] G. Gromov, and V. D. Milman. A topological application of the isoperimetric inequality. Amer. J. Math. 105 (1983), 843–854.
  * [18] R. Kannan, L. Lovász, and M. Simonovits. Isoperimetric problems for convex bodies and a localization lemma. Discrete and Comput. Geom. 13 (1995), 541–559.
  * [19] B. S. Kashin, and A. A Saakyan. Orthogonal series. Translated from the Russian by Ralph P. Boas. Translation edited by Ben Silver. Translations of Mathematical Monographs, 75. American Mathematical Society, Providence, RI, 1989, xii+451 pp.
  * [20] S. Kaczmarz, and G. Steinhaus. Theory of orthogonal series. Warszawa, Lwow, 1935; Russian ed.: Izdat. Fiz.-Mat. Lit., Moscow, 1958, 507 pp.
  * [21] B. Klartag. A Berry-Esseen type inequality for convex bodies with an unconditional basis. Probab. Theory Related Fields 145 (2009), no. 1–2, 1–33.
  * [22] B. Klartag, and J. Lehec. Bourgain’s slicing problem and KLS isoperimetry up to polylog. arXiv:2203.15551v2, 10 April 2022.
  * [23] B. Klartag, and S. Sodin. Variations on the Berry-Esseen theorem. Teor. Veroyatn. Primen. 56 (2011), no. 3, 514–533; reprinted in: Theory Probab. Appl. 56 (2012), no. 3, 403–419.
  * [24] M. Ledoux. On Talagrand’s deviation inequalities for product measures. ESAIM Probab. Statist. 1 (1995/97), 63–87.
  * [25] M. Ledoux. Concentration of measure and logarithmic Sobolev inequalities. Séminaire de Probabilités XXXIII. Lect. Notes in Math. 1709 (1999), 120–216, Springer.
  * [26] P. Massart. The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality. Ann. Probab. 18 (1990), no. 3, 1269–1283.
  * [27] V. V. Petrov. Sums of independent random variables. Translated from the Russian by A. A. Brown. Ergebnisse der Mathematik und ihrer Grenzgebiete, Band 82. Springer–Verlag, New York–Heidelberg, 1975. x+346 pp.
  * [28] V. V. Petrov. Limit theorems for sums of independent random variables (Russian), Nauka, Moscow, 1987. 318 pp.
  * [29] R. Salem, and A. Zygmund. On lacunary trigonometric systems. Proc. Nat. Acad. Sci. USA, 33 (1947), 333–338.
  * [30] R. Salem, and A. Zygmund. On lacunary trigonometric series. II. Proc. Nat. Acad. Sci. U. S. A. 34 (1948), 54–62.
  * [31] V. N. Sudakov. Typical distributions of linear functionals in finite-dimensional spaces of high dimension. (Russian) Soviet Math. Dokl. 19 (1978), 1578–1582; translation in: Dokl. Akad. Nauk SSSR, 243 (1978), no. 6, 1402–1405.