# Label-free quantum super-resolution imaging using entangled multi-mode
squeezed light

Daniel Soh<EMAIL_ADDRESS>Sandia National Laboratories, Livermore,
California 94550  Eric Chatterjee<EMAIL_ADDRESS>Sandia National
Laboratories, Livermore, California 94550

###### Abstract

In this study, we explore the theoretical application of entangled multi-mode
squeezed light for label-free optical super-resolution imaging. By generating
massively entangled multi-mode squeezed light through an array of balanced
beam splitters, using a single-mode squeezed light input, we create a multi-
mode quantum light state with exceptional entanglement and noise suppression
below the shot noise level. This significantly reduces imaging measurement
errors compared to classical coherent state light imaging when the same number
of photons are used on the imaging sample. We demonstrate how to optimize the
imaging system’s parameters to achieve the Heisenberg imaging error limit,
taking into account the number of entangled modes and photons used. We also
examine the effects of optical losses in the imaging system, necessitating
adjustments to the optimized parameters based on the degree of optical loss.
In practical applications, this new quantum imaging approach reduces the
number of photons needed to achieve the same image quality by two orders of
magnitude compared to classical imaging methods that use non-entangled, non-
squeezed coherent state light.

##  I Introduction

Super-resolution imaging is a technique showing an object’s visual features
smaller than the so-called Rayleigh diffraction limit, which is approximately
on the order of probing light’s wavelength. Typical classical (as opposed to
quantum) super-resolution imaging relies on a labeling technique where
optically sensitive materials (dyes) are injected into samples, and either
active or passive fluorescence imaging techniques are used to localize the
optically responding labels, revealing nanoscale features. This kind of
technique includes photoactivated localization microscopy (PALM), stochastic
optical reconstruction microscopy (STORM), stimulated emission depletion
(STED) microscopy, and structured illumination microscopy (SIM) (for a
comprehensive review on these techniques, see [1, 2]). However, these
classical techniques have shortfalls, including limited target biomarkers
(only a few at a time) and enormously long imaging acquisition time. Most of
all, sample preparation (labeling) is time-consuming, and the toxicity of dyes
on the benign imaging samples is not well understood.

There are other non-optical imaging methods that allow even a sub-nanometer
resolution. These include electron microscopy (EM) [3] and soft x-ray
tomography (SXT) [4]. However, these use high energy particles as imaging
probes, which may result in undesired side-effects: the high-energy particles
may disturb the sample’s quantum states (e.g., ionization, elevation to higher
energy states) or even destroy fragile samples such as sensitive chemical
molecules or biosamples. For these sensitive samples, optical imaging provides
an excellent alternative method due to the low energy optical photons leading
to less-interfering imaging, which however is limited in imaging resolution
due to the finite size of photon wavepacket.

Many attempts were made for label-free optical imaging. Recognizing that the
ultimate imaging resolution is limited by the quantum fluctuation of light
(i.e., shot noise) [5], multiple research results focused on how to reduce the
effect of quantum noise in the probing light [6, 7]. These included the usage
of quantum entanglement of single photons [8, 9] and squeezed light [10, 11,
12, 13]. More recently, quantum sensing in general observed a broader use of
the combination of quantum entanglement and squeezing in order to enhance the
signal-to-noise ratio of measurement under the condition of low photon numbers
[14, 15]. In this context, a single-mode squeezed light is used to generate
massively entangled multi-mode squeezed light (MEMSL) that bears the excellent
feature of reducing measurement errors in the Heisenberg limit with respect to
the number of entangled modes [16]. So far, the related research has focused
on measurements of lumped elements. Inspired by this new development of the
combination of quantum entanglement and squeezing, we present for the first
time a new quantum imaging technique that is label-free utililzing the MEMSL
as probe light to enhance the image resolution. MEMSL imaging has a superior
performance with a significantly reduced imaging error in the Heisenberg limit
where the imaging error scales as $1/MN$ ($M$: number of entangled modes in
MEMSL, $N$: number of photons in each mode), rather than the standard limit
$1/\sqrt{MN}$ of classical imaging utilizing non-entangled coherent light
imaging repeating $M$-times. While one can achieve super-resolution imaging
even with a classical method using a coherent state of light, provided that
sufficient averaging reduces the measurement noise, the new MEMSL imaging
accomplishes the same level of imaging resolution with much reduced
requirement for averaging, which therefore accomplishes significantly shorter
imaging acquisition time with a much smaller number of photons used.

There are many cases where the number of photons impinging on the imaging
sample must be limited: quantum state tomography of molecules requires only a
small number of probing photons in order not to disturb the quantum state to
be measured [17, 18, 19]. Bioimaging, in general, requires only a small number
of photons in order not to disturb the bioprocesses or overheat the fragile
biomolecules [20]. Furthermore, chemical surface reactive intermediate
molecules are so fragile that excessive probe light may dissociate or ionize
them, making them hard to detect [21]. All these cases demand limiting the
number of photons impinging on the sensitive and fragile imaging samples.
Hence, it is imperative to study how one can accomplish the best imaging
resolution under such a strict limitation on number of photons impinging on
the samples. Our new MEMSL imaging addresses this important question by
providing a breakthrough quantum imaging utilizing entanglement and squeezing.

In this paper, we start with introducing MEMSL and how it is generated, which
is followed by the theory of quantum MEMSL imaging that is compared to two
other cases, namely, non-entangled single-mode squeezed light imaging
repeating multiple times and classical imaging using coherent state light
repeating multiple times. Then, we discuss the obtainable image resolution
with a limited number of photons. We also discuss the impact of optical losses
in the imaging system. We then provide numerical case studies and discuss how
effective MEMSL imaging is in realistic situations.

##  II Quantum Imaging Theory

Figure 1: A method to generate a 4-mode MEMSL using an array of balanced beam
splitters with an input from a single-mode squeezed light.

In this section, we briefly introduce MEMSL that will play an important role
in our scheme and discuss how one can utilize MEMSL as a probe for quantum
imaging. In order to create the MEMSL light state, one starts with a single-
mode squeezed light. A network of beam-splitters produces a multi-mode output
when only one input is seeded by the single-mode squeezed light while all the
other inputs are in vacuum state [14, 15]. For example, let us consider a
4-mode MEMSL as shown in Fig. 1, which uses three 50-50 beam splitters. It is
well known that the multi-mode output is quantum mechanically entangled,
resulting in a relation [16]

$b_{e}=\frac{1}{2}(b^{(1)}+b^{(2)}+b^{(3)}+b^{(4)}),$ (1)

where $b^{(i)}$ represents the annihilation operator of the $i$th output port
mode and $b_{e}$ is the annihilation operator of the input single-mode
squeezed light. The above equation has a semi-classical interpretation that
the field injected in the squeezed light input port is divided into the four
output modes equally. Quantum mechanically, the photons in the input squeezed
mode are divided into four modes with equal probabilities. In general, it can
be shown that the following holds for MEMSL with $M$ modes [16]:

$b_{e}=\frac{1}{\sqrt{M}}\sum_{m=1}^{M}b^{(m)}.$ (2)

Figure 2: A time-delayed feeding scheme for MEMSL’s individual modes to the
imaging station.

The multiple modes in MEMSL with M modes are used in the imaging scheme one
mode at a time (Fig. 2). In this scheme, each output mode of $M$-mode MEMSL is
selected sequentially using a switching array with differential time delays.
Hence, the $m$th mode will be fed into the imaging station after $m\Delta T$
time delay where $\Delta T$ is appropriately set between two adjacent
measurements in the imaging station. The measurement is taken for each of the
$m$ modes, and later processed together.

Figure 3: A one-dimensional 4-$f$ imaging system.

Following the treatment in Kolobov and Fabre [5], we will consider a one-
dimensional imaging system shown in Fig. 3 where an object with size $Y$ is
placed in the object plane (coordinate $y$). Generalization of MEMSL imaging
to higher dimensions is rather straightforward (see Section VI. Discussions
and conclusions for more details). A pair of lenses with focal length $f$ and
pupil size $d$ are used to image the object on the image plane (coordinate
$x$). Since the lenses have only a finite pupil size, higher spatial-frequency
components in the pupil plane are lost. As a result, the image on the image
plane has smooth boundaries comprising of only lower spatial-frequency
components [5]. Let us normalize the coordinates such that
$s=2x/Y,s^{\prime}=2y/Y$. This normalization allows the generalization of the
theory to other two-lens imaging systems with different focal lengths where
$s^{\prime}=2y/(Yf^{\prime}/f)$ with $f^{\prime}$ representing the focal
length of the second (imaging plane side) lens. We will also define a space-
bandwidth product

$c=\frac{\pi dY}{2\lambda f},$ (3)

where $\lambda$ is the wavelength of the light.

In the following, we will introduce three imaging methods, namely, (1) MEMSL
imaging, (2) single-mode squeezed light imaging repeating $M$ times, and (3)
classical imaging using coherent-state light repeating $M$ times, in order to
demonstrate the superior performance of MEMSL imaging over the two others.

###  II.1 MEMSL imaging method

We emphasize that we extend the quantum imaging much further by introducing
the optimization of squeezing and displacement parameters in the context of
applying the MEMSL technique, while the Kolobov and Fabre’s results described
merely the advantage of using squeezed light in imaging [5]. Our imaging
method accomplishes a Heisenberg imaging error reduction with respect to the
number of photons impinging on the imaging sample and the number of modes in
MEMSL, which is a significant milestone for general quantum imaging.

Let us assume a plane wave for each of the MEMSL output modes. Our imaging
scheme measures any combination of light phase and amplitude using a homodyne
or a heterodyne measurement (or any superposition between the two) after the
probing light interacts with the object. We consider, without losing
generality of analysis, the case where the object’s spatial information is
imprinted on the optical phase of the interacting light. The case where the
object’s spatial information is imprinted in the amplitude of light (or any
combination of amplitude and phase) can be treated similarly, just by rotating
the squeezing angle of the squeezed light in the phase plane appropriately for
minimizing the quantum noise in the desired light quadrature angle.

We apply a MEMSL treatment and finally reconstruct the objective in a new way
that fully utilizes the quantum entanglement in MEMSL to achieve the minimal
error in the image reconstruction (see the general imaging treatment and the
definitions of prolate spheroidal functions $\varphi_{j}(s^{\prime})$,
$\chi_{j}(s^{\prime})$, and $\psi(s)$ in Appendix A). For this, let’s define
the estimate object field (in an orthogonal decomposition form)
$\tilde{a}^{\mathrm{avg}}(s^{\prime})=\sum_{j=0}^{Q}\tilde{a}_{j}^{\mathrm{avg}}\varphi_{j}(s^{\prime})$
for $|s^{\prime}|\leq 1$. Here, we limit the summation to $j=Q$, which is the
highest order of useful prolate spheroidal wave functions. The orders higher
than $Q$ do more harm than good because the quantities
$\langle\tilde{a}_{j>Q}^{\mathrm{avg}}\rangle$ are noisy. The maximum order
number $Q$ is determined by both quantum and measurement noise. It will be
shown later that $Q$ is directly related to the obtainable resolution of the
reconstructed object field because higher orders of the prolate spheroidal
wave functions provide faster spatial variation (oscillation) of the object
field features (i.e., a better imaging resolution of the reconstructed object
field is obtained if higher orders of prolate spheroidal functions are added).

We now show that an average of $e_{j}^{(m)}$ over all values of the mode index
$m$ in MEMSL will accomplish the minimum imaging error. Eq. (73) results in a
simple relation for the estimate
$\tilde{a}_{j}^{\mathrm{avg}}=(1/\sqrt{\lambda_{j}})e_{j}^{\mathrm{avg}}$
where $e_{j}^{\mathrm{avg}}=(1/M)\sum_{m=1}^{M}e_{j}^{(m)}$: the estimate
object field for $|s^{\prime}|\leq 1$ is

$\displaystyle\tilde{a}^{\mathrm{avg}}(s^{\prime})=\sum_{j=0}^{Q}\tilde{a}_{j}^{\mathrm{avg}}\varphi_{j}(s^{\prime})=\sum_{j=1}^{Q}\frac{e_{j}^{\mathrm{avg}}}{\lambda_{j}}\psi_{j}(s^{\prime}).$
(4)

Again, the ultimate goal of imaging is to obtain $\phi(s^{\prime})$, which is
the spatially-dependent phase distribution (Eq. (57)). We note the following
crucial relation between $\tilde{a}^{\mathrm{avg}}(s^{\prime})$ and
$\phi(s^{\prime})$, which is derived using the characteristic MEMSL relation
in Eq. (2):

$\tilde{a}^{\mathrm{avg}}(s^{\prime})=\frac{1}{M}\sum_{m=1}^{M}b^{(m)}e^{i\phi(s^{\prime})}=\frac{1}{\sqrt{M}}b_{e}e^{i\phi(s^{\prime})}.$
(5)

Combed with Eq. (4), this relation extracts $\phi(s^{\prime})$ from the
measured $e_{j}^{\mathrm{avg}}$.

Let us now assume that $\phi(s^{\prime})$ is sufficiently small for all
$|s^{\prime}|\leq 1$. We introduce the quadrature operators for $b_{e}$ as

$\displaystyle b_{e}$ $\displaystyle=b_{e1}+ib_{e2},$ $\displaystyle b_{e1}$
$\displaystyle=(b_{e}+b_{e}^{\dagger})/2,$ $\displaystyle b_{e2}$
$\displaystyle=(b_{e}-b_{e}^{\dagger})/(2i).$ (6)

The quadrature operators are Hermitian. We also introduce the quadrature
operators for $b^{(m)}$ similarly: $b_{1}^{(m)}=(b^{(m)}+b^{(m)\dagger})/2$
and $b_{2}^{(m)}=(b^{(m)}-b^{(m)\dagger})/(2i)$. We assume a displaced
squeezed state for $b_{e}$-field such that

$\displaystyle\langle b_{e1}\rangle=\sqrt{2}\alpha_{e},\quad\langle
b_{e2}\rangle=0,$ $\displaystyle\langle(\Delta
b_{e1})^{2}\rangle=\frac{e^{+2r_{e}}}{4},\quad\langle(\Delta
b_{e2})^{2}\rangle=\frac{e^{-2r_{e}}}{4},$ (7)

where $\alpha_{e}$ is the displacement and $r_{e}$ is the squeezing parameter.
We note that this state is a “phase squeezed” state, suitable for phase
estimate. Then, one obtains from Eq. (5)

$\langle\tilde{a}_{2}^{\mathrm{avg}}(s^{\prime})\rangle=\frac{\sqrt{2}\alpha_{e}}{\sqrt{M}}\sin\phi(s^{\prime})\approx\alpha_{e}\phi(s^{\prime})\sqrt{\frac{2}{M}},$
(8)

where
$\tilde{a}_{2}^{\mathrm{avg}}(s^{\prime})=(\tilde{a}^{\mathrm{avg}}(s^{\prime})-\tilde{a}^{\mathrm{avg\dagger}}(s^{\prime}))/(2i)$
is a quadrature operator. Based on these and Eq. (4), we define the estimate
$\tilde{\phi}(s^{\prime})$, which is a Hermitian quantum field operator, as
following:

$\displaystyle\tilde{\phi}(s^{\prime})$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\sum_{j=0}^{Q}\frac{1}{\alpha_{e}\lambda_{j}}\sqrt{\frac{M}{2}}\psi_{j}(s^{\prime})e^{(m)}_{j2}$
$\displaystyle=\frac{1}{\alpha_{e}\sqrt{2M}}\sum_{j=0}^{Q}\frac{\left(\sum_{m=1}^{M}e^{(m)}_{j2}\right)}{\lambda_{j}}\psi_{j}(s^{\prime}),$
(9)

where $e_{j2}^{(m)}=(e^{(m)}_{j}-e^{(m)\dagger}_{j})/(2i)$, representing the
imaginary part of the image component created by the light in mode $m$ and
corresponding to the $j$th spheroidal function.

The goal of imaging is to obtain the best estimate
$\langle\tilde{\phi}(s^{\prime})\rangle$ under the influence of the quantum
noise embedded in $e_{j2}^{(m)}$. This noise in turn relates to the noise
embedded in the source photon field $b_{2}^{(m)}$ through the following (see
Appendix E for the derivation):

$e^{(m)}_{j2}=b_{2}^{(m)}\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s)+b_{1}^{(m)}\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\psi_{j}(s^{\prime}),$
(10)

leading to the following for the composite imaginary image component for the
$j$th spheroidal function:

$\displaystyle\sum_{m=1}^{M}e^{(m)}_{j2}=$ $\displaystyle
b_{e2}\sqrt{M}\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s^{\prime})+b_{e1}\sqrt{M}\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\psi_{j}(s^{\prime}).$
(11)

The noise in the reconstructed object field, resulting from the noise in the
image field, becomes the following (see Appendix E):

$\displaystyle\begin{split}\sigma&=\sqrt{\int_{-\infty}^{\infty}ds^{\prime}~{}\langle(\Delta\tilde{\phi}(s^{\prime}))^{2}\rangle}\\\
&\approx\frac{e^{-r_{e}}}{2\sqrt{2}\alpha_{e}}\sqrt{\sum_{j=0}^{Q}\frac{A_{j}^{2}}{\lambda_{j}^{2}}},\end{split}$
(12)

which is the root-mean-squared sum of all noises in the reconstructued
spatially-dependent phase. The coefficients $A_{j}$ are defined as follows:

$A_{j}=\left|\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s^{\prime})\right|,$
(13)

Note that the reconstructed-object-field noise associated with a given basis
function $j$ scales inversely with $\lambda_{j}$. This is because for a lower
$\lambda_{j}$, the object-to-image attenuation is more pronounced. As such,
when deriving the reconstructed object field from the image field, a greater
amplification is required, which also causes significant amplification of the
noise from the image field. Based on the fact that $\lambda_{j}$ drops rapidly
with increasing $j$, we can make the approximation
$A_{Q}^{2}/\lambda_{Q}^{2}\gg\sum_{j=1}^{Q-1}A_{j}^{2}/\lambda_{j}^{2}$, which
is numerically verifiable. This leads to

$\sigma\approx\frac{e^{-r_{e}}A_{Q}}{2\sqrt{2}\alpha_{e}\lambda_{Q}}.$ (14)

We now consider the fact that the number of photons impinging on the sample
must be limited. The average number of photons hitting the sample in each mode
of MEMSL is given as

$N_{e}=\frac{1}{M}(\alpha_{e}^{2}+\sinh^{2}r_{e}),$ (15)

since the photons in the single-mode squeezed light are equally divided into
the $M$ modes of MEMSL.

Now, we introduce the optimization for obtaining the minimum imaging error
$\sigma$ by varying the squeezing parameter $r_{e}$ and the displacement
$\alpha_{e}$ while fixing $N_{e}=N$, for a given number of modes $M$ in the
MEMSL. Substituting Eq. (15) into Eq. (14), we can express the amplitude-
squared of the noise, i.e. $\sigma^{2}$, in terms of the squeezing parameter
$r_{e}$ as follows:

$\sigma^{2}\approx\frac{A_{Q}^{2}}{8\lambda_{Q}^{2}}\bigg{(}\frac{e^{-2r_{e}}}{MN-\sinh^{2}{r_{e}}}\bigg{)}.$
(16)

Taking the derivative of $\sigma^{2}$ with respect to $e^{-2r_{e}}$ and
setting it to zero (corresponding to the minimum value for $\sigma^{2}$), the
optimal solution set is easily obtained to be

$e^{-r_{e}}=\frac{1}{\sqrt{1+2MN}},\quad\alpha_{e}=\frac{\sqrt{MN(1+MN)}}{\sqrt{1+2MN}}.$
(17)

The tradeoff can be conceptualized as follows: The phase uncertainty per
photon (i.e., the angle subtended by the second-quadrature uncertainty in the
phase space distribution) approximately equals the ratio between the second-
quadrature uncertainty and the amplitude. As the squeezing parameter $r_{e}$
is increased, the second-quadrature uncertainty per photon decreases as
desired, but the amplitude per photon also decreases. Consequently, if the
light is too strongly squeezed, the phase uncertainty per photon (i.e., the
phase noise) will start to increase, requiring optimization. With the above
solution, the optimal (minimal) imaging error $\sigma$ is

$\sigma^{\mathrm{opt}}=\left(\frac{A_{Q}}{2\sqrt{2}\lambda_{Q}}\right)\frac{1}{\sqrt{MN(1+MN)}}.$
(18)

We note that, for a large $MN$, the error scales as $\sim 1/MN$, which is a
Heisenberg limit with respect to both $M$ and $N$.

###  II.2 Non-entangled squeezed light imaging with $M$-time repeating
measurements

To compare the MEMSL imaging technique above with other alternatives, we
introduce non-entangled single-mode squeezed light imaging. In this method, a
single squeezed light is used for $M$ times and the result is averaged.

Let us designate the annihilation operator for each non-entangled single-mode
squeezed light as $b_{e}^{(m)}=b_{e1}^{(m)}+ib_{e2}^{(m)}$ for
$m=1,2,\cdots,M$ with $b_{e1}^{(m)}=(b_{e}^{(m)}+b_{e}^{(m)\dagger})/2$ and
$b_{e2}^{(m)}=(b_{e}^{(m)}-b_{e}^{(m)\dagger})/(2i)$. The statistics of the
quadrature operators is

$\displaystyle\langle
b_{e1}^{(m)}\rangle=\sqrt{2}\alpha^{\prime}_{e},\quad\langle
b_{e2}^{(m)}\rangle=0,$ $\displaystyle\langle(\Delta
b_{e1}^{(m)})^{2}\rangle=\frac{e^{+2r^{\prime}_{e}}}{4},\quad\langle(\Delta
b_{e2}^{(m)})^{2}\rangle=\frac{e^{-2r^{\prime}_{e}}}{4},$ (19)

It is straightforward to derive

$\tilde{a}_{2}^{\mathrm{avg}}(s^{\prime})\approx\sqrt{2}\alpha^{\prime}_{e}\tilde{\phi}(s^{\prime}),$
(20)

which leads to

$\tilde{\phi}(s^{\prime})=\frac{1}{\alpha_{e}M\sqrt{2}}\sum_{j=0}^{Q}\frac{\left(\sum_{m=1}^{M}e_{j2}^{(m)}\right)}{\lambda_{j}}\psi_{j}(s^{\prime}),$
(21) $\displaystyle\sum_{m=1}^{M}e_{j2}^{(m)}=$
$\displaystyle\left(\sum_{m=1}^{M}b_{e2}^{(m)}\right)\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s^{\prime})$
$\displaystyle+\left(\sum_{m=1}^{M}b_{e1}^{(m)}\right)\int_{-1}^{1}ds^{\prime}\tilde{\phi}(s^{\prime})\psi_{j}(s^{\prime}),$
(22)

and finally, one obtains

$\sigma\approx\frac{e^{-r^{\prime}_{e}}A_{Q}}{2\sqrt{2M}\alpha^{\prime}_{e}\lambda_{Q}}.$
(23)

This implies that the error in this case scales as $\sim 1/\sqrt{M}$, which is
a well known result for an average over multiple non-entangled measurement
trials. Also, we note that the photon number in each trial amounts to
$N_{e}={\alpha^{\prime}_{e}}^{2}+\sinh^{2}r^{\prime}_{e}$. Substituting this
into Eq. (23), $\sigma^{2}$ becomes

$\sigma^{2}\approx\frac{A_{Q}^{2}}{8M\lambda_{Q}^{2}}\bigg{(}\frac{e^{-2r^{\prime}_{e}}}{N-\sinh^{2}{r^{\prime}_{e}}}\bigg{)}$
(24)

The optimal solutions of $\alpha^{\prime}_{e},r^{\prime}_{e}$ that minimize
$\sigma^{2}$ are

$e^{-r^{\prime}_{e}}=\frac{1}{\sqrt{1+2N}},\quad\alpha^{\prime}_{e}=\frac{\sqrt{N(1+N)}}{\sqrt{1+2N}},$
(25)

which leads to the optimal (minimal) imaging error

$\sigma^{\mathrm{opt}}=\left(\frac{A_{Q}}{2\sqrt{2}\lambda_{Q}}\right)\frac{1}{\sqrt{M}\sqrt{N(1+N)}}.$
(26)

We note that for large $N$, the error scales as $\sim 1/(N\sqrt{M})$, which is
a standard limit with respect to $M$ but a Heisenberg limit with respect to
$N$. This is due to the fact that in the absence of entanglement, squeezing
degrades the amplitude per photon more rapidly, with the rate of degradation
with respect to squeezing increasing with $M$. Consequently, without
entanglement, the noise decreases more slowly with the number of modes $M$.
This result reveals that the entanglement in different modes of MEMSL is
essential to make the $M$-dependence be a Heisenberg limit.

###  II.3 Non-entangled coherent state light imaging with $M$-time repeating
measurements

Next, we calculate the imaging error when classical non-entangled non-squeezed
coherent state light is used for $M$ times and the result is averaged.

This case is quite similar to the previous case of non-entangled squeezed
light imaging except $r^{\prime}_{e}=0$. We designate the annihilation
operator of the individual coherent-state light as
$b_{e}^{(m)}=b_{e1}^{(m)}+ib_{e2}^{(m)}$ where

$\displaystyle\langle b_{e1}^{(m)}\rangle=\sqrt{2}\alpha^{\prime}_{e},\langle
b_{e2}^{(m)}\rangle=0,$ $\displaystyle\langle(\Delta
b_{e1}^{(m)})^{2}\rangle=\langle(\Delta b_{e2}^{(m)})^{2}\rangle=\frac{1}{4},$
(27)

it is straightforward to obtain, replacing $r^{\prime}_{e}=0$ in Eq. (23):

$\sigma\approx\frac{A_{Q}}{2\sqrt{2M}\alpha^{\prime}_{e}\lambda_{Q}}.$ (28)

Then, with the fact that the number of photons in individual imaging trials is
$N_{e}={\alpha^{\prime}_{e}}^{2}=N$, the above turns into

$\sigma\approx\left(\frac{A_{Q}}{2\sqrt{2}\lambda_{Q}}\right)\frac{1}{\sqrt{MN}},$
(29)

which implies that the error scales as the standard limit for both $M$ and
$N$. Note that there is no need to optimize $\alpha^{\prime}_{e}$ here because
it is set by the number of photons used in each measurement. It is also worth
noting that the same result is obtained even for entangled coherent state
light. This reveals that squeezing and optimization of the squeezing parameter
$r_{e}$ are necessary in order to make the $N$-dependence be a Heisenberg
limit.

##  III Obtainable image resolution

In the following, we will perform an estimate of imaging resolution. We show
that the obtainable image resolution is a direct consequence of the integer
$Q$, which is the highest-order index of the prolate spheroidal wave functions
that are used to reconstruct the object field in the previous section. Recall
that the minimum image noise from the MEMSL imaging technique is given in Eq.
(18), which is inversely proportional to $\lambda_{Q}$. We note that, if the
index $Q$ increases by one, $\lambda_{Q}$ decreases significantly (orders of
magnitude for a sufficiently large $Q$), which demands to limit $Q$. To
determine $Q$, we cut the highest order index $Q$ where the accumulated
imaging noise with the summation of the prolate spheroidal functions up to
$(Q-1)$st order equals the total number of photons used (i.e., the signal-to-
noise ratio is unity). It is due to the observation that higher $Q$ will
suffer more from measurement noise since higher-order prolate spheroidal
functions have faster spatial oscillations, causing significantly sharper
object-to-image attenuation (hence requiring greater amplification when
reconstructing the object from the image, which also increases the
amplification of the noise). We also exploit the fact that the noise is zero
for for odd spheroidal functions (see the definition of $A_{j}$ in equation
(13) and the parity of the odd orders of prolate spheroidal functions
(Appendix B)). As a result, by using an odd value of $Q$ as the cut-off, we
can obtain a resolution corresponding to that odd value of $Q$ while the noise
corresponds only to the even $(Q-1)$st-order spheroidal function. Therefore,
we limit $Q$ such that

$\sigma_{Q-1}^{2}\approx MN,$ (30)

where we used the fact that the noise energy equals $\sigma^{2}$. Then, we
obtain the condition to determine $Q$ in MEMSL imaging such that

$\frac{A_{Q-1}^{2}}{\lambda_{Q-1}^{2}}\lessapprox 8M^{2}N^{2}(1+MN).$ (31)

One chooses the maximum integer $Q$ that satisfies this inequality.

Once the value of $Q$ is determined using the above, the obtainable resolution
$D$ is determined by considering the fact that the prolate spheroidal function
of order $Q$ ($\psi_{Q}(s^{\prime})$) has approximately $Q$ number of zeros in
the interval $[-1,1]$ [22] and $Q+1$ number of mountains and valleys in the
same interval (see Fig. 8). Hence, the image resolution is given as

$D\simeq\frac{Y}{Q+1}=R\frac{S}{Q+1},$ (32)

where $S=2c/\pi=dY/\lambda f$ is the Shannon number, and $R=Y/S=\lambda f/d$
is the Rayleigh resolution length.

##  IV Impact of optical losses

Optical losses degrade both the noise-squeezing and the entanglement in MEMSL.
It is because optical losses couple the probing optical fields used in the
imaging to the environments, leading to leakage of information in an
irreversible way. Therefore, optical losses will degrade the imaging
performance. In this section, we analyze the impact of optical losses on the
imaging noise.

In Guo et al., the loss is modeled as a lumped single beam splitter located
between the light source and the sample [16], which models loss mechanisms in
various locations into a single lumped loss. In modeling point of view, this
approach is much more attractive than other alternatives because it allows a
tractable treatment of the optical loss. The position of the lumped loss
element is not important as one can scale the model beam splitter’s ratio
appropriately according to the actual loss amount. We follow this convention,
modeling the optical loss as the probing light undergoing a beam splitting
action:

$b^{(m)}\rightarrow\sqrt{\tau}b^{(m)}+\sqrt{1-\tau}d^{(m)}_{v},$ (33)

where $d^{(m)}_{v}$ is the annihilation operator of the external vacuum mode
that is introduced via the beam splitter. Here, $\tau$ is the beam splitter’s
splitting ratio, which models the optical loss of $1-\tau$. We note that
$\langle d^{(m)}_{v}\rangle=0$ and $\langle(\Delta
d^{(m)}_{v1,v2})^{2}\rangle=1/4$ where
$d^{(m)}_{v1}=(d^{(m)}_{v}+d^{(m)\dagger}_{v})/2$ and
$d^{(m)}_{v2}=(d^{(m)}_{v}-d^{(m)\dagger}_{v})/(2i)$.

With this, we proceed to modify the previous analysis, now incorporating the
optical loss. First, we study the impact of optical loss on MEMSL imaging: we
modify the recipe Eq. (5) into the following, after reflecting the fact that
only $\sqrt{\tau}$ ratio of the probe light will arrive at the object after
the loss:

$\tilde{a}^{\mathrm{avg}}(s^{\prime})=\frac{\sqrt{\tau}}{M}\sum_{m=1}^{M}b^{(m)}e^{i\tilde{\phi}(s^{\prime})}=\sqrt{\frac{\tau}{M}}b_{e}e^{i\tilde{\phi}(s^{\prime})}.$
(34)

This leads to

$\displaystyle\tilde{a}_{2}^{\mathrm{avg}}(s^{\prime})$
$\displaystyle\approx\alpha_{e}\sqrt{\frac{2\tau}{M}}\tilde{\phi}(s^{\prime}),$
$\displaystyle\tilde{\phi}(s^{\prime})$
$\displaystyle=\frac{1}{\alpha_{e}\sqrt{2\tau
M}}\sum_{j=0}^{Q}\frac{\left(\sum_{m=1}^{M}e^{(m)}_{j2}\right)}{\lambda_{j}}\psi_{j}(s^{\prime}).$
(35)

In addition, Eq. (10) is modified to

$\displaystyle e^{(m)}_{j2}=$
$\displaystyle\left(\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s^{\prime})\right)\left(\sqrt{\tau}b^{(m)}_{2}+\sqrt{1-\tau}d^{(m)}_{v2}\right)$
$\displaystyle+\left(\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\psi_{j}(s^{\prime})\right)\left(\sqrt{\tau}b^{(m)}_{1}+\sqrt{1-\tau}d^{(m)}_{v1}\right).$
(36)

Therefore, we calculate

$\displaystyle\int_{-1}^{1}ds^{\prime}\langle(\Delta\tilde{\phi}(s^{\prime}))^{2}\rangle$
$\displaystyle=\frac{1}{2\tau\alpha_{e}^{2}M}\sum_{j=0}^{Q}\frac{\langle\Delta(\sum_{m=1}^{M}e^{(m)}_{j2})^{2}\rangle}{\lambda_{j}^{2}}$
$\displaystyle\approx\frac{1}{2\tau\alpha_{e}^{2}M}\times$
$\displaystyle\sum_{j=0}^{Q}\frac{A_{j}^{2}\left(M\tau\langle(\Delta
b_{e2})^{2}\rangle+(1-\tau)\sum_{m=1}^{M}\langle(\Delta
d^{(m)})^{2}\rangle\right)}{\lambda_{j}^{2}}$
$\displaystyle=\frac{e^{-2r_{e}}+1/\tau-1}{8\alpha_{e}^{2}}\sum_{j=0}^{Q}\frac{A_{j}^{2}}{\lambda_{j}^{2}}$
$\displaystyle\approx\frac{A_{Q}^{2}}{8\alpha_{e}^{2}\lambda_{Q}^{2}}\left(e^{-2r_{e}}+1/\tau-1\right).$
(37)

Hence, we obtain a modified imaging error of lossy MEMSL imaging:

$\sigma\approx\frac{A_{Q}}{2\sqrt{2}\alpha_{e}\lambda_{Q}}\sqrt{e^{-2r_{e}}+1/\tau-1}.$
(38)

We note that the above reduces to the lossless MEMSl imaging error in Eq. (14)
as $\tau\rightarrow 1$ (i.e., lossless). We note that this is a monotonically
decreasing function of $\tau$, implying that the imaging error reduces
monotonically as $\tau$ increases (i.e., the loss becomes smaller).

As the loss increases, the number of photons in each mode on the sample
reduces. If we are interested in setting the same number of photons $N$ on the
sample regardless of the loss, we will have to increase either the
displacement $\alpha_{e}$ or the squeezing parameter $r_{e}$ as the number of
photons is given as Eq. (15). It is possible to re-optimize the system
parameters by re-designing the squeezing parameter $r_{e}$ and the
displacement $\alpha_{e}$ to minimize the imaging error for given optical loss
$1-\tau$ so that one can still have the same number of photons $N$ impinging
on the sample, but now with much higher initial optical power. If one designs
the system to have the same $N$, it is intuitively expected that the imaging
error may relatively stay to be similar.

With the optical loss, we modify the number of photons in each branch of MEMSL
since only $\tau$ ratio of photons will arrive at the object after the loss:

$N_{e}=\frac{\tau}{M}(\alpha_{e}^{2}+\sinh^{2}r_{e}).$ (39)

We then optimize $\alpha_{e}$ and $r_{e}$ while fixing $N_{e}=N$ as we
discussed above, yielding the following optimal solutions that minimize the
noise $\sigma$:

$\displaystyle e^{-r_{e}}=\sqrt{\frac{\Lambda+\tau}{1+4MN+\tau}},$
$\displaystyle\alpha_{e}=$
$\displaystyle\frac{1}{\sqrt{2}}\sqrt{\frac{8M^{2}N^{2}(1-\tau)-\tau(\Lambda-1)-2MN(\tau(\Lambda+2\tau-2)-1)}{(1-\tau)\tau(1+4MN+\tau)}},$
(40)

where $\Lambda=\sqrt{1+4MN(1-\tau)}$. Therefore, the optimal (minimal) imaging
error is given as

$\sigma^{\mathrm{opt}}=\frac{A_{Q}}{2\lambda_{Q}}\sqrt{\frac{(1-\tau)\left(1+\tau\Lambda\right)}{8M^{2}N^{2}(1-\tau)+2MN(1-\tau(\Lambda+2(\tau-1)))+\tau(1-\Lambda)}}.$
(41)

It is easy to verify that the above approaches the result of a lossless case
in Eq. (18) as $\tau\rightarrow 1$.

Next, we study the impact of optical loss in non-entangled squeezed light
imaging with $M$-time repeating measurements. The analysis is similar, but
having a crucial difference so that we now have the imaging error

$\sigma\approx\frac{A_{Q}}{2\sqrt{2M}\alpha_{e}^{\prime}\lambda_{Q}}\sqrt{e^{-2r_{e}^{\prime}}+1/\tau-1}.$
(42)

The optimal solutions are

$\displaystyle
e^{-r_{e}^{\prime}}=\sqrt{\frac{\Lambda^{\prime}+\tau}{1+4N+\tau}},$
$\displaystyle\alpha_{e}^{\prime}=$
$\displaystyle\frac{1}{\sqrt{2}}\sqrt{\frac{8N^{2}(1-\tau)-\tau(\Lambda^{\prime}-1)-2N(\tau(\Lambda^{\prime}+2\tau-2)-1)}{(1-\tau)\tau(1+4N+\tau)}},$
(43)

where $\Lambda^{\prime}=\sqrt{1+4N(1-\tau)}$. These lead to the optimal
imaging error of

$\sigma^{\mathrm{opt}}=\frac{A_{Q}}{2\lambda_{Q}\sqrt{M}}\sqrt{\frac{(1-\tau)\left(1+\tau\Lambda^{\prime}\right)}{8N^{2}(1-\tau)+2N(1-\tau(\Lambda^{\prime}+2(\tau-1)))+\tau(1-\Lambda^{\prime})}}.$
(44)

It is easy to verify that this solution approaches that in Eq. (26) as
$\tau\rightarrow 1$ (i.e., lossless).

Finally, for the case of classical coherent light imaging with $M$-time
repeating measurement, the number of photons used requires
$\tau{\alpha^{\prime}_{e}}^{2}=N$. Hence, with a reduced $\tau$, one must
increase $\alpha^{\prime}_{e}$ to maintain the same number of photons
impinging on the sample. However, the imaging error does not depend on $\tau$
because the number of photons measured are maintained to be identical
regardless of $\tau$. In this case, the imaging error in Eq. (29) remains the
same regardless of optical loss.

##  V Numerical examples and analysis

We now present numerical examples of our quantum super-resolution scheme. For
numerical evaluations of prolate spheroidal functions and their eigenvalues,
we use the following [23]:

$\displaystyle\lambda_{n}(c)$ $\displaystyle=\frac{2c}{\pi}R_{0n}^{2}(c,1),$
$\displaystyle\psi_{n}(c,x)$
$\displaystyle=\sqrt{\frac{\lambda_{n}(c)}{\int_{-1}^{1}dx^{\prime}S^{2}_{0n}(c,x^{\prime})}}S_{0n}(c,x),$
(45)

where $R_{mn}(c,x),S_{mn}(c,x)$ are the radial and the angular solution with
order $m,n$ of the Helmholtz wave equation of the first kind, respectively. We
explicitly made $c$ as a parameter of the functions indicating that both the
eigenvalues $\lambda_{n}$ and the prolate spheroidal functions $\psi_{n}$ are
functions of $c$. The following relation is extremely useful to numerically
evaluate $A_{j}$ in Eq. (13) [23] :

$\int_{-\infty}^{\infty}dt\psi_{n}(c,t)e^{i\omega
t}=i^{n}\sqrt{\frac{2\pi}{c}}\psi_{n}\left(c,\frac{\omega}{c}\right),$ (46)

which leads to

$A_{j}=\left|\int_{-\infty}^{\infty}ds\psi_{j}(c,s)\right|=\sqrt{\frac{2\pi}{c}}\left|\psi_{j}(c,0)\right|.$
(47)

For calculating the values of $\psi_{n}(c,x)$, we used the built-in functions
of a python package (Scipy), namely, $R_{mn}(c,x)=$
scipy.special.pro_rad1(m,n,c,x) and $S_{mn}(c,x)=$
scipy.special.pro_ang1(m,n,c,x), respectively. These numerical packages are
able to evaluate the cases for $|x|<1$. For the numerical evaluation for
$|x|\geq 1$, we utilize the Eq. (63):

$\psi_{j}(s)=\frac{1}{\lambda_{j}}\int_{-1}^{1}ds^{\prime}\,\psi_{j}(s^{\prime})\frac{\sin[c(s^{\prime}-s)]}{\pi(s^{\prime}-s)}.$
(48)

We first demonstrate the utility of the MEMSL imaging technique where we
assume lenses with a focal length $f=10$ mm, the wavelength of light
$\lambda=780$ nm, and a lens diameter $d$ of 2 inches. In this case, the
Rayleigh diffraction limit is $\lambda f/d=154$ nm. We also assume that the
object has a size $Y=300$ nm. These parameter values imply $c=\pi dY/2\lambda
f=3.07$. We then assume that the MEMSL has eight entangled modes ($M=8$), and
each mode has 6 photons ($N=6$). In this case, the optimal squeezing and
displacement of the input single-mode squeezed mode to MEMSL generator from
Eq. (17) are $e^{-r_{e}}=0.10$ and $\alpha_{e}=4.92$. Also, this example
implies that the highest order number $Q$ according to the Eq. (31) is 7
(since the noise energy corresponding to $Q-1=6$ is on the same order as the
total photon number), implying that the resolution one can achieve is 38 nm
according to Eq. (32). This example of MEMSL imaging beats the Rayleigh
diffraction limit deeply by roughly a factor of 5 with such a small number of
photons used (total 48 photons). On the other hand, in order to accomplish the
same level of imaging error, the classical imaging using non-entangled non-
squeezed coherent state light with $M$-time repeating measurement would
require 294 photons each time for averaging over eight times (total 2352
photons), which is roughly two orders-of-magnitude larger number of photons on
the sample than MEMSL imaging. Therefore, we learn that MEMSL imaging
significantly reduces the required number of photons.

In order to simulate the imaging system with quantum noise, we utilize the Eq.
(10), which is refined using Eq. (46):

$e^{(m)}_{j2}=b_{2}^{(m)}\sqrt{\frac{2\pi}{c}}\psi_{j}(c,0)+b_{1}^{(m)}\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\psi_{j}(s^{\prime}).$
(49)

The quantum noise enters via $b_{2}^{(m)}$ and $b_{1}^{(m)}$, which are
quadrature operators. A more informative quantity is the average
$(1/M)\sum_{m=1}^{M}e_{j2}^{(m)}$ since this quantity is directly used to
reconstruct $\tilde{\phi}(s^{\prime})$, the original object’s spatial phase
distribution. For the case of MEMSL imaging, we use the Eq. (11), where the
noise is given in Eq. (7) with Gaussian noise statistics. For coherent light
imaging method, we use Gaussian noise with $\langle
b_{1}^{(m)}\rangle=\sqrt{2}\alpha^{\prime}$, $\langle b_{2}^{(m)}\rangle=0$,
$\langle(\Delta b_{1,2})^{2}\rangle=1/4$. The image field’s second quadrature
is given as
$e_{2}^{(m)}(s^{\prime})=\sum_{j=1}^{\infty}e_{j2}^{(m)}\psi_{j}(s^{\prime})$,
which now contains the appropriate quantum noise for different cases of MEMSL
imaging or non-entangled squeezed/coherent light imaging with $M$-time
repeating measurement.

An alternative numerical simulation approach is to directly generate random
data for
$e_{2}^{\mathrm{avg}}(s^{\prime})=(1/M)\sum_{m=1}^{M}e_{2}^{(m)}(s^{\prime})$
for each point $s^{\prime}$. Note that, from Eq. (71), we obtain the average
field:

$\langle
e_{2}^{\mathrm{avg}}(s^{\prime})\rangle=\int_{-1}^{1}ds^{\prime}\frac{\sin[c(s^{\prime}-s)]}{\pi(s^{\prime}-s)}\left(\frac{1}{M}\sum_{m=1}^{M}\langle
a_{2}^{(m)}(s^{\prime})\rangle\right).$ (50)

With $\langle a_{2}^{(m)}(s^{\prime})\rangle\simeq\langle
b_{2}^{\mathrm{(m)}}\rangle\phi(s^{\prime})$ (see Eq. (57)), the above average
is easily calculated. For the calculation of the variance $\langle(\Delta
e_{2}^{\mathrm{avg}}(s^{\prime}))^{2}\rangle$, we separate the case of MEMSL
imaging and the case of the coherent light imaging repeating $M$-times. For
MEMSL imaging, the calculation is

$\displaystyle\langle(\Delta e_{2}^{\mathrm{avg}}(s^{\prime}))^{2}\rangle$
$\displaystyle=\frac{1}{M^{2}}\times$
$\displaystyle\sum_{m=1}^{M}\sum_{j,j^{\prime}=0}^{\infty}\left(\langle
e_{2j}^{(m)}e_{2j^{\prime}}^{(m)}\rangle-\langle e_{2j}^{(m)}\rangle\langle
e_{2j^{\prime}}^{(m)}\rangle\right)\psi_{j}(s^{\prime})\psi_{j^{\prime}}(s^{\prime})$
$\displaystyle=\frac{1}{M^{2}}\sum_{m=1}^{M}\sum_{j=0}^{\infty}\langle(\Delta
e_{2j}^{(m)})^{2}\rangle\psi_{j}^{2}(c,s^{\prime})$
$\displaystyle\simeq\frac{\langle(\Delta
b_{2e})^{2}\rangle}{M}\sum_{j=0}^{\infty}\frac{2\pi}{c}\psi_{j}^{2}(c,0)\psi_{j}^{2}(c,s^{\prime})$
$\displaystyle=\frac{G(c,s^{\prime})}{4M(1+2MN)},$ (51)

where

$G(c,s^{\prime})=\sum_{j=0}^{\infty}\frac{2\pi}{c}\psi_{2j}^{2}(c,0)\psi_{2j}^{2}(c,s^{\prime}).$
(52)

Here, the second equation in (51) used the fact that $e_{2j}$ and
$e_{2j^{\prime}}$ are independent random variables unless $j=j^{\prime}$. The
third equation used the MEMSL relation in (2) and the fact that the noise
contribution of the second term in Eq. (49) is negligible. The last line used
the optimal solution for MEMSL in Eq. (17) and the fact that $\psi_{j}(c,0)=0$
for odd $j$s (see Fig. 8 (b)). It is known that the following sum converges
(see Eq. (78) of Moore and Cada [23]):

$\sum_{j=0}^{\infty}\frac{2\pi}{c}\psi_{2j}^{2}(c,0)=E^{c},$ (53)

where $E^{c}$ is a constant. Hence, $G(c,s^{\prime})$ also converges because
$\psi_{2j}^{2}(c,s^{\prime})\leq\psi_{c,0}^{2}(0)$ for all $j\geq 1$. In fact,
the summation quickly converges numerically (the numerical accuracy is
sufficient even if we sum only up to $j=20$). Numerical evaluations of the
function $G(c,s^{\prime})$ are shown in Fig. 4. Higher $c$ values concentrate
more towards the center. Incidentally, the domain of $s^{\prime}$ where
$G(c,s^{\prime})$ is large coincides with region of dense information. Since
both the expected value and the variance of the Gaussian random variable
$e_{2}^{\mathrm{avg}}(s^{\prime})$ are known for each $s^{\prime}$, we can
numerically simulate the entire noisy $e_{2}^{\mathrm{avg}}(s^{\prime})$ for
all $s^{\prime}$.

Figure 4: Numerical evaluation of $G(c,s^{\prime})$ for various $c$ values.

For coherent light imaging repeating $M$ times, the following holds:

$\displaystyle\langle(\Delta e_{2}^{\mathrm{avg}}(s^{\prime}))^{2}\rangle$
$\displaystyle\simeq\frac{1}{M^{2}}\sum_{m=1}^{M}\sum_{j=0}^{\infty}\langle(\Delta
b_{2}^{(m)})^{2}\rangle\frac{2\pi}{c}\psi_{j}^{2}(c,0)\psi_{j}^{2}(c,s^{\prime})$
$\displaystyle=\frac{G(c,s^{\prime})}{4M}.$ (54)

Likewise, because we know both the expected value and the variance of
$e_{2}^{\mathrm{avg}}(s^{\prime})$ for all points $s^{\prime}$ for the
classical imaging using coherent state light, we can numerically simulate the
entire noisy $e_{2}^{\mathrm{avg}}(s^{\prime})$ for all $s^{\prime}$.

Let us compare the above two results in Eq. (51) and (54): MEMSL imaging
significantly reduces the error in each measurement data
$e_{2}^{\mathrm{avg}}(s^{\prime})$ at all points $s^{\prime}$ because the
noise variance is inversely proportional to roughly $M^{2}$, which is due to
the entanglement in the MEMSL branch modes, and also to the number of photons
$N$ in each MEMSL branch mode. In contrast, the coherent light imaging lacks
this feature, and the error in $e_{2}^{\mathrm{avg}}(s^{\prime})$ is inversely
proportional to only $M$, due to lack of entanglement in the independent
$M$-time measurement. It is remarkable to observe that the number of photons
$N$ does not reduce the error of $e_{2}^{\mathrm{avg}}(s^{\prime})$ at each
individual point $s^{\prime}$ for the classical imaging using coherent state
light, while it reduces the accumulated error $\sigma$ (the phase distribution
error, see Eq. (29)). This can be understood because the error in phase
distribution estimate follows $\Delta\tilde{\phi}(s)\propto\Delta
a^{\mathrm{avg}}(s)/\sqrt{N}$ (c.f., Eq. (8)). So, while the phase
distribution estimate error reduces as $1/\sqrt{N}$ in the coherent light
imaging repeating $M$-times, the measured value $a^{\mathrm{avg}}(s)$, which
is obtained from $e^{\mathrm{avg}}(s^{\prime})$, is not affected by $N$.

We note that achieving high-level squeezing is hard. For example, the best
measured noise squeezing result so far is 15 dB [24]. In imaging, the actual
measured level of squeezing matters because only imaging performance depends
on measured data. According to the optimal squeezing solution in Eq. (17), the
optimal squeezed noise power ($e^{-2r_{e}}$) is desired to be $1/(1+2MN)$. The
15 dB squeezing implies only, for example, four photons ($N=4$) in each mode
for eight modes ($M=8$). The hardship in accomplishing large squeezing forces
averaging over multiple independent measurements, rather than accomplishing an
excessively large level of squeezing. In contrast, the coherent light imaging
repeating $M$-times does not see any difference in the roles of $M$ and $N$
(see Eq. (29)), which makes sense physically since repeating more times using
coherent light is exactly same as adding more photons in a single measurement
as the independent photons reduce the noise as $\sim 1/\sqrt{N}$. For
numerical simulations, we assume a fictitious one-dimensional object plane’s
phase distribution as shown in Fig. 5 (dashed-line), which is designed
intentionally to be asymmetric to accommodate both the even and the odd orders
of prolate spheroidal functions. Then, we performed numerical simulation of
the image plane fields for both MEMSL imaging with $M=8$ and $N=6$
(corresponding to 20-dB-squeezed input light (i.e.,
$10\log_{10}(e^{-2r_{e}})=-20$) to the array of balanced beam splitters to
generate the MEMSL light) (which used the same $M=8$ and $N=6$), both
averaging over 50000 independent measurements. Fig. 5 shows numerically
simulated image plane quadrature distribution $e_{2}(s)$ with added quantum
noise for MEMSL imaging and coherent light imaging. It is clearly shown that
the MEMSL imaging’s (simulated) measured field data shows significantly
reduced noise compared to the coherent light imaging.

((a))

((b))

Figure 5: Image field quadrature values (solid lines) with added quantum noise
resulting from the object’s phase distribution (dashed lines). The measurement
is repeated for $50000$ times with $M=4$ and $N=12$ for both (a) MEMSL imaging
and (b) coherent light imaging. The imaging system’s space-bandwidth is
$c=3.07$.

We follow the reconstruction recipe in Eq. (9) for estimating
$\tilde{\phi}(s^{\prime})$. For the maximum index $Q$, we follow the criterion
in Eq. (30) for both MEMSL and classical imaging, which are good for a single-
shot object field estimate. For averaging $N^{\mathrm{avg}}$-time independent
repeating measurements, the equation must slightly change to
$\sigma\rightarrow\sigma/\sqrt{N^{\mathrm{avg}}}$ as the repeating
measurements are independent (i.e., unentangled). Also, the total number of
photons used must change $MN\rightarrow MNN^{\mathrm{avg}}$. These change the
maximum $Q$ criterion to be

$\frac{A_{Q-1}^{2}}{\lambda_{Q-1}^{2}}\leq
8M^{2}N^{2}N^{\mathrm{avg2}}(1+MN).$ (55)

According to this, the estimated maximum $Q$ is 7 for MEMSL. On the other
hand, for the coherent light imaging repeating $M$-times, we use the Eq. (29)
for applying to Eq. (30)

$\frac{A_{Q-1}^{2}}{\lambda_{Q-1}^{2}}\leq 8M^{2}N^{2}N^{\mathrm{avg2}},$ (56)

which produces the estimated maximum $Q=5$. We found that adding a higher-
order prolate spheroidal functions made the estimate
$\tilde{\phi}(s^{\prime})$ quickly diverging due to the denominator
$\sqrt{\lambda_{j}}$ in Eq. (73), which decreases fast by multiple order-of-
magnitudes with a slightly increased index $j$. Fig. 6 shows the reconstructed
object plane’s optical phase distribution $\tilde{\phi}(s^{\prime})$ for both
MEMSL imaging and coherent light imaging. While the reconstructed optical
phase distribution from MEMSL imaging accurately shows three lobes with
somewhat inaccurate heights, that of the classical light imaging completely
missed the salient features of the object’s optical phase distribution (i.e.,
wrong number of lobes, completely wrong heights).

We note that our numerical example demonstrated the efficacy of our MEMSL
quantum imaging for probing nanoscale structures: revealing such a complex and
fine structure with size far below the Rayleigh diffraction limit demands the
imaging performance to beat the Rayleigh diffraction limit deeply. The
accomplished imaging resolution from MEMSL imaging in the numerical example is
roughly 38 nm ($=300/(7+1)$ with $Q=7$) whereas the Rayleigh diffraction limit
in this case was 154 nm. We note that the obtainable resolution can be
improved by several means. Having more entangled modes ($M$) with increased
number of photons ($N$) improves the optical resolution at the cost of large
squeezing in the input single-mode squeezed light for MEMSL generator. In
practice, increasing the number of independent measurement trials (i.e.,
increasing $N^{\mathrm{avg}}$) is a viable solution.

We emphasize that, with the same number of modes $(M)$ and the same number of
photons $(N)$ for a single measurement, MEMSL imaging is always superior to
the classical light imaging due to the combined benefit from the squeezing and
the entanglement in MEMSL. On the contrary, the effect of increasing
$N^{\mathrm{avg}}$ is basically identical between MEMSL and classical light
imaging since the measurements are independent from each other for every
measurement trial.

Figure 6: Reconstructed object plane’s optical phase distributions
$\tilde{\phi}(s^{\prime})$ from MEMSL imaging (solid line) and classical light
imaging (dotted line), along with the original accurate test object plane’s
optical phase distribution $\phi(s^{\prime})$ (dashed line).

((a))

((b))

Figure 7: (a) Comparison of optimized imaging error $\sigma$ from classical
coherent light imaging (dashed), non-entangled squeezed light imaging
(dotted), and MEMSL imaging (solid) as functions of light transmission $\tau$.
(b) Optimal values of displacement $\alpha_{e}^{2}$ and initial squeezing
($10\log_{10}(e^{-2r_{e}})$) of MEMSL imaging and non-entangled squeezed light
imaging as a function of light transmission $\tau$.

Next, we compare the optimized imaging errors $\sigma$ of MEMSL imaging, non-
entangled single-mode squeezed light imaging repeating multiple times, and
classical imaging using non-entangled non-squeezed coherent light imaging for
various values of $\tau$ (the optical loss is given as ($1-\tau$)) (see Fig.
7). It is clear that the optimized imaging error of MEMSL imaging is
approximately an order-of-magnitude smaller than the classical coherent light
imaging for various values of $\tau$. The imaging error of MEMSL imaging is
also consistently better than that of non-entangled repeating squeezed light
imaging as well. While the imaging error values stay relatively similar for
different $\tau$ values as expected, it is somewhat counter-intuitive to
observe that lower transmission values (higher optical losses) lead to smaller
optimized imaging error for both MEMSL imaging and non-entangled squeezed
light imaging repeating multiple times.

We also compare the optimal configuration of MEMSL imaging and non-entangled
repeating squeezed light imaging (Fig. 7 (b)). We note that the displacement
$\alpha_{e}$ and the level of squeezing $e^{-2r_{e}}$ of MEMSL imaging
represent the properties of the single-mode squeezed light used as an input
for MEMSL generation while $\alpha^{\prime}_{e}$ and $e^{-2r^{\prime}_{e}}$ of
non-entangled squeezed light imaging repeating multiple times are properties
of the individual light source for each measurement, which is somewhat similar
to the individual mode of MEMSL. Therefore, it is difficult to give a direct
fair comparison between $\alpha_{e}$ and $\alpha^{\prime}_{e}$, and also
$e^{-2r_{e}}$ and $e^{-2r^{\prime}_{e}}$. With this caveat, we find that
$\alpha_{e}^{2}$ is nearly the same as $M{\alpha^{\prime}}_{e}^{2}$ for all
values of $\tau$, which makes good physical sense, since dividing the single-
mode input squeezed light into $M$ modes corresponds to the individual
squeezed light in the non-entangled repeating squeezed light imaging.
Naturally, the optimal configuration for $alpha_{e}^{2}$ and
${\alpha^{\prime}_{e}}^{2}$ increases as $\tau$ decreases (i.e., optical loss
increases). In contrast, the optimal configuration for squeezing is reverse,
indicating that the role of squeezing becomes less important as the loss
increases, and the imaging performance is rapidly compensated by large
displacement values $(\alpha_{e}^{2},{\alpha^{\prime}_{e}}^{2})$, which is in
consistence with the expectation that optical loss will quickly degrade the
benefit coming from squeezing. It is also important to note that the required
squeezing of MEMSL is much more demanding than that of non-entangled squeezed
light imaging repeating multiple times as shown clearly in the graph in Fig. 7
(b). The reason for this is that the squeezing in the single-mode squeezed
input for generating MEMSL is a resource to create the entanglement among the
$M$ modes of MEMSL, which will accomplish much reduced imaging error in MEMSL
imaging than non-entangled repeating squeezed light imaging. Therefore, the
extra degree of squeezing requirement of MEMSL is the price to pay to
accomplish a better imaging performance.

##  VI Discussions and conclusion

We presented a new label-free quantum imaging technique using MEMSL as probe
light. It is clearly shown that MEMSL’s imaging performance is superior
compared to other methods presented here, namely, non-entangled single-mode
squeezed light imaging repeating multiple times and classical imaging using
coherent state light repeating multiple times. The MEMSL imaging collects the
measurement data for each mode in MEMSL and combines them later to reconstruct
the original object information. Due to quantum entanglement and squeezing,
the imaging error in MEMSL image accomplishes the Heisenberg limit where the
error scales as $1/MN$ ($M$: number of entanglement modes in MEMSL and $N$:
number of photons impinging on the sample in each MEMSL mode), while that of
non-entangled squeezed light imaging repeating multiple times scales as
$1/N\sqrt{M}$, and that of classical imaging using coherent state light
repeating multiple times scales as the standard limit $1/\sqrt{NM}$. The
Heisenberg limit $1/M$ is ascribed to the quantum entanglement among $M$-modes
in MEMSL whereas the Heisenberg limit $1/N$ is ascribed to the squeezing.

The performance of MEMSL imaging in real world is limited by how much
squeezing one can accomplish. The achievable squeezing limits the product of
$MN$ in the optimal imaging configurations (see Eq. (17)). With the practical
squeezing level of $15-20$ dB, we demonstrated that the imaging error of MEMSL
imaging is still approximately an order of magnitude better than that of using
coherent state light repeating multiple times. This reduced imaging noise led
to more faithful reconstruction of the object information: MEMSL imaging
accurately reconstructed quite complex object features with a resolution that
deeply beats the Rayleigh diffraction limit while classical imaging using
coherent state light repeating multiple times fails to extract such complex
object features, when both imaging methods used the same number of modes ($M$)
and the same number of photons ($N$) impinging upon the imaging sample. In the
case of MEMSL imaging with $c=3.07$, $M=8$, $N=6$, classical imaging using
coherent state light repeating multiple times will need approximately two
orders-of-magnitude larger number of photons in order to accomplish the same
imaging error as MEMSL imaging.

We also provided a comprehensive optimization strategy when imaging system
suffers from optical losses. The price to pay in lossy imaging system is to
use a squeezed light input to MEMSL generator with a larger displacement
(i.e., higher light power) with somewhat less degree of squeezing than the
optimal lossless case. We showed that such an arrangement accomplishes the
same image quality under the same limit for the number of photons impinging on
the sample as the lossless case. While a lossless imaging system sees equal
contributions in the number of photons (i.e.,
$N_{e}=\alpha_{e}^{2}+\sinh^{2}r_{e}$) from the displacement
($\alpha_{e}^{2}$) and squeezing ($\sinh^{2}r_{e}$) to accomplish a good
imaging quality, lossy imaging systems rely more on the displacement
($\alpha_{e}^{2}$) than squeezing ($\sinh^{2}r_{e}$) since squeezing quickly
degrades with respect to optical losses.

While we derived the theory based on a one-dimensional imaging using one-
dimensional prolate spheroidal (Slepian) functions, we note that it is
straightforward to adopt higher-dimension prolate spheroidal functions [25,
26, 27] to expand our theory to higher-dimensional imaging. Since two-
dimensional prolate spheroidal functions also satisfy the crucial integral
relation similar to Eq. (69), which guarantees that the two-dimensional
prolate spheroidal functions are still the eigenfunctions after the two-
dimensional integrals, all the theory stands roughly the same, except changing
the one-dimensional integrals to two-dimensional integrals. With this caveat,
our one-dimensional analysis can be appropriately adapted for a two-
dimensional imaging application. Furthermore, it is important to mention that,
because we normalized the coordidnates $s,s^{\prime}$ in the object and the
image plane, our theory also can be applied to different two-lens imaging
systems that have disparate focal lengths, after appropriately scaling
normalized coordinates $s,s^{\prime}$.

## Acknowledgements

We thank Eric Chatterjee and Philiip Chrostoski for various valuable
discussions.

The work was supported by the Laboratory Directed Research and Development
program at Sandia National Laboratories, a multimission laboratory managed and
operated by National Technology and Engineering Solutions of Sandia, LLC., a
wholly owned subsidiary of Honeywell International, Inc., for the U.S.
Department of Energy’s National Nuclear Security Administration under Contract
No. DE-NA-003525. This paper describes objective technical results and
analysis. Any subjective views or opinions that might be expressed in the
paper do not necessarily represent the views of the U.S. Department of Energy
or the United States Government.

##  Appendix A Prolate Spheroidal Functions and Recovery of Objects

In this section, we briefly present how one recovers the object from the
imaging plane measurements, largely based on the results in [5] (see Fig. 3).

Passing through the object, the optical field acquires spatially-dependent
phases $\phi(s^{\prime})$: for example, the $m$-th mode in MEMSL becomes

$a^{(m)}(s^{\prime})=b^{(m)}e^{i\phi(s^{\prime})},\quad|s^{\prime}|\leq 1.$
(57)

The goal of imaging is to obtain $\phi(s^{\prime})$ accurately, allowing for
the reconstruction of the object at every point $s^{\prime}$. The object-to-
image transformation is given as follows (see Appendix C for a derivation):

$L[\psi](s)=\int_{-\infty}^{\infty}ds^{\prime}\psi(s^{\prime})\frac{\sin{(c(s^{\prime}-s))}}{\pi(s^{\prime}-s)},$
(58)

where $c$ is the previously defined space bandwidth product. Note that this
accords with the expression given by Kolobov and Fabre [5]. In order to
reconstruct the object from the image information, we decompose the object
field into basis fields. We start with the in-pupil eigenfields $\varphi_{j}$
of the object-to-image transformation $L$, which satisfy the following
relationship for $|s|\leq 1$:

$L[\varphi_{j}](s)=\lambda_{j}\varphi_{j}(s).$ (59)

In the small-object limit, however, the image size extends well beyond the
object size, as previously discussed. For each in-pupil field $\varphi_{j}$,
we define the out-pupil field $\chi_{j}$ to represent the part of the image in
the region $|s|>1$, less a proportionality constant that is a function of
$\lambda_{j}$:

$\chi_{j}(s)=\frac{L[\varphi_{j}](s)}{\sqrt{\lambda_{j}(1-\lambda_{j})}}.$
(60)

For a classical light source, the object could simply be decomposed into the
basis of the in-pupil series $\varphi_{j}$, since no scattering point exists
for $|s^{\prime}|>1$. However, due to quantum fluctuation, it is necessary to
include a superposition of the out-pupil series $\chi_{j}$ as well, with the
corresponding coefficients featuring an expectation value of zero. Therefore,
the object plane field is decomposed as follows:

$a^{(m)}(s^{\prime})=\sum_{j=0}^{\infty}\left(a_{j}^{(m)}\varphi_{j}(s^{\prime})+c_{j}^{(m)}\chi_{j}(s^{\prime})\right),-\infty<s^{\prime}<\infty.$
(61)

Here, $a_{j}^{(m)}$ and $c_{j}^{(m)}$ are superposing weights for $j$th order
orthogonal functions $\varphi_{j}(s^{\prime})$ and $\chi_{j}(s^{\prime})$,
respectively. Note that $\langle c_{j}^{(m)}\rangle=0$ for all values of $j$.
In general, Eq. (59) is solved by the prolate spheroidal wave functions
(a.k.a., Slepian functions) $\psi_{j}(s^{\prime})$, which satisfy the
following conditions [5]:

$\displaystyle L[\psi](s)=\psi(s),$ (62)
$\displaystyle\int_{-1}^{1}ds^{\prime}\psi_{j}(s^{\prime})\frac{\sin{(c(s^{\prime}-s))}}{\pi(s^{\prime}-s)}=\lambda_{j}\psi_{j}(s).$
(63)

The orthogonal functions $\varphi_{j}(s^{\prime})$ and $\chi_{j}(s^{\prime})$
relate to $\psi_{j}(s^{\prime})$ as follows:

$\displaystyle\varphi_{j}(s^{\prime})=\left\\{\begin{array}[]{ll}\frac{\psi_{j}(s^{\prime})}{\sqrt{\lambda_{j}}},&|s^{\prime}|\leq
1,\\\ 0,&|s^{\prime}|>1.\end{array}\right.$ (66)
$\displaystyle\chi_{j}(s^{\prime})=\left\\{\begin{array}[]{ll}0,&|s^{\prime}|\leq
1,\\\
\frac{\psi_{j}(s^{\prime})}{\sqrt{1-\lambda_{j}}},&|s^{\prime}|>1\end{array}\right.$
(69)

where $c=\pi dY/2\lambda f$ is the space-bandwidth product in Eq. (3). Here,
$\lambda_{j}$’s are the eigenvalues of $j$th order prolate spheroidal wave
functions under the integral defined in Eq. (63) with the property
$1\geq\lambda_{0}>\lambda_{1}>\cdots>0$. Appendix B shows the eigenvalues
$\lambda_{j}$ and both the odd and the even orders of the prolate spheroidal
functions, which are all functions of the value $c$.

For quantum optics treatment, we elevate the superposing weights $a_{j}^{(m)}$
and $c_{j}^{(m)}$ to bosonic annihilation operators, satisfying commutation
relations $[a_{j}^{(m)},a_{k}^{(m)\dagger}]=\delta_{j,k}$ and
$[c_{j}^{(m)},c_{k}^{(m)\dagger}]=\delta_{j,k}$ as is well justified in
Kobolov and Fabre [5]. Then, $a^{(m)}(s^{\prime})$ becomes a quantum field
operator. We also note that the image plane field is decomposed as

$e^{(m)}(s)=\sum_{j=0}^{\infty}e_{j}^{(m)}\psi_{j}(s),\quad-\infty<s<\infty.$
(70)

We note the relation between the object plane field and the image plane field
([5], where the integral domain is not explicitly written after the authors
added the quantum treatment including the domain $|s^{\prime}|>1$, but we make
it explicit and now add the outside region in the integral):

$e^{(m)}(s)=L[a^{(m)}](s).$ (71)

The image superposition coefficients (which we elevate to bosonic annihilation
operators) thus take the following form in terms of the object coefficients
(see the derivation in Appendix D):

$e^{(m)}_{j}=\sqrt{\lambda_{j}}a^{(m)}_{j}+\sqrt{1-\lambda_{j}}c^{(m)}_{j}.$
(72)

One has to estimate the average object field
$a^{\mathrm{avg}}(s^{\prime})=(1/M)\sum_{m=1}^{M}a^{(m)}(s^{\prime})$ after
one takes the measurement of image field $e^{(m)}(s)$, for all
$m=1,2,\cdots,M$. The recipe to obtain $a^{(m)}_{j}$ for each of
$a^{(m)}(s^{\prime})$ using Eq. (72) is, following the treatment of Kolobov
and Fabre [5]:

$\tilde{a}_{j}^{(m)}=\frac{\tilde{e}_{j}^{(m)}}{\sqrt{\lambda_{j}}}.$ (73)

The above equation is based on the fact that the expected values of Eq. (72)
lead to $\langle e_{j}^{(m)}\rangle=\sqrt{\lambda_{j}}\langle
a_{j}^{(m)}\rangle$ since $\langle c_{j}^{(m)}\rangle=0$ for all $j$ because
values of $a^{(m)}(s^{\prime})$ outside the domain $|s^{\prime}|\leq 1$ is
zero (i.e., the object only exists in the range $|s^{\prime}|\leq 1$, and the
coefficients $c_{j}^{(m)}$ represent quantum noise and thus cancel out when
averaged).

##  Appendix B Prolate Spheroidal Functions and Eigenvalues

Figure 8 shows the even and odd prolate spheroidal functions, along with the
eigenvalues corresponding to the object-to-image transformation for each of
the functions.

((a))

((b))

((c))

Figure 8: Prolate spheroidal functions. (a) Eigenvalues, (b) even orders
($c=6$), and (c) odd orders ($c=6$).

##  Appendix C Derivation of Object-to-Image Transformation

We now heuristically derive the object-to-image transformation $L$ given the
setup in Fig. 3, starting with an object with the profile $\psi(s^{\prime})$.
In general, for an object size much smaller than the pupil size yet much
larger than the light wavelength (i.e., $\lambda\ll Y\ll d$, assuming that the
focal length $f$ and pupil size $d$ are on the same order), the size of the
image would approximately equal that of the object, with each point on the
object focusing onto the corresponding point on the image screen. This is due
to the fact that the diffraction-induced spread of the image (on the order of
$\lambda f/d$) is neligible compared to the size of the object $Y$.
Furthermore, the large size of the object relative to the wavelength implies a
narrow beam angle, ensuring that all of the propagating beams are transmitted
through both lenses and onto the screen. Therefore, the intensity of the beams
arriving on the screen would match that of the beams leaving the object.
Consequently, the object-to-image transformation reduces to identity, i.e.
$L[\psi](s)=\psi(s)$.

For the case of a very small object (below the wavelength scale), however, the
pupil size becomes a limiting factor. Specifically, the object size is now
much smaller than the diffraction-induced spread $\lambda f/d$ on the screen,
resulting in the image size being significantly larger than the object size
(i.e., $L[\psi](s)\neq 0$ for $|s|>1$). By the same principle, if we consider
very closely-spaced points on a given object, then the diffraction patterns
corresponding to the different points will overlap and interfere on the image
screen. The one-to-one mapping between object and image is thus lost, and a
more intricate analysis is required to recreate the object from the image
information.

We analyze the image by splitting the object into infinitessimal segments of
length $ds^{\prime}$ and summing the contributions of the beams from each
segment onto the point $s$ on the screen. Since a generic segment at point
$s^{\prime}$ on the object axis gives rise to a field of amplitude
$\psi(s^{\prime})ds^{\prime}$, the overall amplitude of the light focused on
the point $s$ is calculated by summing the diffraction patterns using the
following integral:

$L[\psi](s)\propto\int_{-\infty}^{\infty}ds^{\prime}\psi(s^{\prime})\frac{\sin{\Big{(}\frac{d\pi}{\lambda}\sin{(\theta^{\prime}-\theta)}\Big{)}}}{\frac{d\pi}{\lambda}\sin{(\theta^{\prime}-\theta)}},$
(74)

where $\theta^{\prime}$ represents the angle at which the beams originating
from the object point $s^{\prime}$ traveled prior to reaching the pupil, and
$\theta$ represents the angle at which the beams reaching the screen point $s$
traveled after diffracting through the pupil. Note from Fig. 3 that any beam
traveling through the origin (i.e., through the pupil on the optical axis) at
a generic angle $\theta^{\prime\prime}$ is bent by the second lens such that
it becomes parallel to the optical axis afterward. Therefore, given an object
size much smaller than the pupil size (i.e., $Y\ll d$), we apply the small-
angle approximation to obtain the relationships $\theta=Ys/(2f)$ and
$\theta^{\prime}=Ys^{\prime}/(2f)$. Substituting these into the above
expression for the object-to-image transformation $L$, we find the following:

$L[\psi](s)\propto\int_{-\infty}^{\infty}ds^{\prime}\psi(s^{\prime})\frac{\sin{(c(s^{\prime}-s))}}{c(s^{\prime}-s)},$
(75)

where $c=\pi dY/(2\lambda f)$ is the previously defined space-bandwidth
product.

It is important to note that the fact that the object-to-image transformation
is proportional to the above expression rather than equal to it. This is due
to the fact that in the small-object limit ($dY\ll\lambda f$), the wide
scattering angle of the light from the object and the limited pupil/lens size
ensure that only a fraction of the beams actually enter the lens and make it
through the aperture and the second lens to the screen. Heuristically, we know
that the tangent of the scattering half-angle of the light is on the order of
$\lambda/(2\pi Y)$, since the photon momentum along the propagation axis is
$h/\lambda$ while the limitation on the transverse object size ensures a
transverse momentum uncertainty with a range of about $h/(\pi Y)$. On the
other hand, the tangent of the half-angle subtended by the lens is $d/(2f)$.
Consequently, the fraction of the light amplitude that passes through the lens
is on the order of $\sqrt{\pi dY/(\lambda f)}$ for a one-dimensional setup.
Extrapolating to the case of a two-dimensional setup with circular symmetry,
the amplitude attenuation factor thus becomes about $\pi dY/(\lambda f)$. This
is on the same order as the result found by Kolobov and Fabre [5], i.e.,
$dY/(2\lambda f)=c/\pi$. Multiplying the expression in Eq. (75) by this factor
yields the following equality for the object-to-image transformation:

$L[\psi](s)=\int_{-\infty}^{\infty}ds^{\prime}\psi(s^{\prime})\frac{\sin{(c(s^{\prime}-s))}}{\pi(s^{\prime}-s)},$
(76)

where the term $c$ in the denominator of Eq. (75) has been replaced with
$\pi$.

##  Appendix D Calculating the Image’s Superposition Coefficients

The superposing weight $e_{j}^{(m)}$ is obtained through

$e_{j}^{(m)}=\frac{\int_{-\infty}^{\infty}ds\,e^{(m)}(s)\psi_{j}(s)}{\int_{-\infty}^{\infty}ds\,\psi_{j}^{2}(s)},$
(77)

because $\psi_{j}(s)$ are real functions.

Now, we introduce a linear transformation $L$ such that

$L[f](s)=\int_{-\infty}^{\infty}ds^{\prime}\frac{\sin[c(s-s^{\prime})]}{\pi(s-s^{\prime})}f(s^{\prime}),\quad-\infty<s<\infty.$
(78)

Then, it easily follows that

$L[\varphi_{j}](s)=\sqrt{\lambda_{j}}\psi_{j}(s),\quad
L[\chi_{j}](s^{\prime})=\sqrt{1-\lambda_{j}}\psi_{j}(s^{\prime}).$ (79)

From these, we obtain the following:

$\displaystyle L[a^{(m)}](s)$
$\displaystyle=\sum_{j=0}^{\infty}\left(a^{(m)}_{j}L[\varphi_{j}](s)+c^{(m)}_{j}L[\chi_{j}](s)\right)$
$\displaystyle=\sum_{j=0}^{\infty}\left(\sqrt{\lambda_{j}}a^{(m)}_{j}+\sqrt{1-\lambda_{j}}c^{(m)}_{j}\right)\psi_{j}(s).$
(80)

From Eq. (70), Eq. (71), and the above expressions, we obtain the following
relation:

$e^{(m)}_{j}=\sqrt{\lambda_{j}}a^{(m)}_{j}+\sqrt{1-\lambda_{j}}c^{(m)}_{j},$
(81)

and we can easily show that $[e^{(m)}_{j},e^{(m)\dagger}_{k}]=\delta_{jk}$,
satisfying the quantization condition of the bosonic annihilation operators
$e^{(m)}_{j}$.

##  Appendix E Derivation of Quantum Noise

Here, we analyze the noise in the reconstructed object field. We define the
noise as

$\sigma=\sqrt{\int_{-\infty}^{\infty}ds^{\prime}~{}\langle(\Delta\tilde{\phi}(s^{\prime}))^{2}\rangle},$
(82)

which is the root-mean-squared sum of all noises in the reconstructued
spatially-dependent phase. It is straightforward to derive

$\displaystyle e^{(m)}_{j2}$
$\displaystyle=\mathrm{Im}\left[\sqrt{\lambda_{j}}a_{j}^{(m)}+\sqrt{1-\lambda_{j}}c_{j}^{(m)}\right]$
$\displaystyle=\mathrm{Im}\left[\sqrt{\lambda_{j}}\int_{-1}^{1}ds^{\prime}b^{(m)}e^{i\phi(s^{\prime})}\varphi_{j}(s^{\prime})+\sqrt{1-\lambda_{j}}\int_{|s^{\prime}|>1}ds^{\prime}b^{(m)}\chi_{j}(s^{\prime})\right]$
$\displaystyle\approx\sqrt{\lambda_{j}}\left(b_{2}^{(m)}\int_{-1}^{1}ds^{\prime}\varphi_{j}(s^{\prime})+b_{1}^{(m)}\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\varphi_{j}(s^{\prime})\right)+\sqrt{1-\lambda_{j}}\int_{|s^{\prime}|>1}ds^{\prime}b_{2}^{(m)}\chi_{j}(s^{\prime})$
$\displaystyle=b_{2}^{(m)}\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s)+b_{1}^{(m)}\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\psi_{j}(s^{\prime})$
(83)

Using the MEMSL’s crucial identity in Eq. (2), we obtain

$\displaystyle\sum_{m=1}^{M}e^{(m)}_{j2}=$ $\displaystyle
b_{e2}\sqrt{M}\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s^{\prime})+b_{e1}\sqrt{M}\int_{-1}^{1}ds^{\prime}\phi(s^{\prime})\psi_{j}(s^{\prime}).$
(84)

Now, we apply the Parseval’s theorem to calculate the variance of the noise:

$\displaystyle\int_{-\infty}^{\infty}ds^{\prime}\langle(\Delta\tilde{\phi}(s^{\prime}))^{2}\rangle$
$\displaystyle=\frac{1}{2\alpha_{e}^{2}M}\sum_{j=0}^{Q}\frac{\langle\Delta(\sum_{m=1}^{M}e^{(m)}_{j2})^{2}\rangle}{\lambda_{j}^{2}}\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}^{2}(s^{\prime})$
$\displaystyle=\frac{1}{2\alpha_{e}^{2}M}\sum_{j=0}^{Q}\frac{\langle\Delta(\sum_{m=1}^{M}e^{(m)}_{j2})^{2}\rangle}{\lambda_{j}^{2}}$
$\displaystyle\approx\frac{1}{2\alpha_{e}^{2}M}\sum_{j=0}^{Q}\frac{MA_{j}^{2}\langle(\Delta
b_{e2})^{2}\rangle}{\lambda_{j}^{2}}$
$\displaystyle=\frac{e^{-2r_{e}}}{8\alpha_{e}^{2}}\sum_{j=0}^{Q}\frac{A_{j}^{2}}{\lambda_{j}^{2}},$
(85)

where we defined

$A_{j}=\left|\int_{-\infty}^{\infty}ds^{\prime}\psi_{j}(s^{\prime})\right|.$
(86)

Here, the fourth line comes from the conventional approximation to ignore the
contribution of the second term in Eq. (84) for the variance calculation since
$|\phi(s^{\prime})|\ll 1$ (for example, see Guo et al. [16]). More precisely,
the approximation becomes well justified with a conservative sufficient
condition

$\left|\int_{-1}^{1}ds^{\prime}~{}\phi(s^{\prime})\right|<\frac{e^{-r_{e}}}{\left|\int_{-1}^{1}ds^{\prime}\psi_{0}(s^{\prime})\right|}.$
(87)

## References

  * Cremer  _et al._ [2017] C. Cremer, A. Szczurek, F. Schock, A. Gourram, and U. Birk, Super-resolution microscopy approaches to nuclear nanostructure imaging, Methods 123, 11 (2017).
  * Poole and Mostaço-Guidolin [2021] J. J. Poole and L. B. Mostaço-Guidolin, Optical microscopy and the extracellular matrix structure: a review, Cells 10, 1760 (2021).
  * De Jonge and Ross [2011] N. De Jonge and F. M. Ross, Electron microscopy of specimens in liquid, Nature nanotechnology 6, 695 (2011).
  * Parkinson  _et al._ [2008] D. Y. Parkinson, G. McDermott, L. D. Etkin, M. A. Le Gros, and C. A. Larabell, Quantitative 3-d imaging of eukaryotic cells using soft x-ray tomography, Journal of structural biology 162, 380 (2008).
  * Kolobov and Fabre [2000] M. I. Kolobov and C. Fabre, Quantum limits on optical resolution, Physical review letters 85, 3789 (2000).
  * Beskrovnyy and Kolobov [2005] V. N. Beskrovnyy and M. I. Kolobov, Quantum limits of super-resolution in reconstruction of optical objects, Physical Review A 71, 043802 (2005).
  * Wang  _et al._ [2012] H. Wang, S. Han, and M. I. Kolobov, Quantum limits of super-resolution of optical sparse objects via sparsity constraint, Optics express 20, 23235 (2012).
  * Giovannetti  _et al._ [2004] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum-enhanced measurements: beating the standard quantum limit, Science 306, 1330 (2004).
  * Sajia and Qian [2021] A. Sajia and X.-F. Qian, Super-resolution of two unbalanced point sources assisted by the entangled partner, arXiv preprint arXiv:2112.01725 (2021).
  * Treps  _et al._ [2005] N. Treps, V. Delaubert, A. Maître, J.-M. Courty, and C. Fabre, Quantum noise in multipixel image processing, Physical Review A 71, 013820 (2005).
  * Embrey  _et al._ [2015] C. Embrey, M. Turnbull, P. Petrov, and V. Boyer, Observation of localized multi-spatial-mode quadrature squeezing, Physical Review X 5, 031004 (2015).
  * Treps  _et al._ [2004] N. Treps, N. Grosse, W. P. Bowen, M. Hsu, A. Maître, C. Fabre, H. Bachor, and P. K. Lam, Nano-displacement measurements using spatially multimode squeezed light, Journal of Optics B: Quantum and Semiclassical Optics 6, S664 (2004).
  * Feng  _et al._ [2013] F. Feng, S. wen Bi, B. zhu Lu, and M. hua Kang, Long-term stable bright amplitude-squeezed state of light at 1064 nm for quantum imaging, Optik 124, 1070 (2013).
  * Zhuang  _et al._ [2018] Q. Zhuang, Z. Zhang, and J. H. Shapiro, Distributed quantum sensing using continuous-variable multipartite entanglement, Physical Review A 97, 032329 (2018).
  * Zhang and Zhuang [2021] Z. Zhang and Q. Zhuang, Distributed quantum sensing, Quantum Science and Technology 6, 043001 (2021).
  * Guo  _et al._ [2020] X. Guo, C. R. Breum, J. Borregaard, S. Izumi, M. V. Larsen, T. Gehring, M. Christandl, J. S. Neergaard-Nielsen, and U. L. Andersen, Distributed quantum sensing in a continuous-variable entangled network, Nature Physics 16, 281 (2020).
  * Zhang  _et al._ [2021] M. Zhang, S. Zhang, Y. Xiong, H. Zhang, A. A. Ischenko, O. Vendrell, X. Dong, X. Mu, M. Centurion, H. Xu, _et al._ , Quantum state tomography of molecules by ultrafast diffraction, Nature communications 12, 1 (2021).
  * Mouritzen and Mølmer [2006] A. S. Mouritzen and K. Mølmer, Quantum state tomography of molecular rotation, The Journal of chemical physics 124, 244311 (2006).
  * Skovsen  _et al._ [2003] E. Skovsen, H. Stapelfeldt, S. Juhl, and K. Mølmer, Quantum state tomography of dissociating molecules, Physical review letters 91, 090406 (2003).
  * Voronin and Zheltikov [2010] A. Voronin and A. Zheltikov, Ionization penalty in nonlinear optical bioimaging, Physical Review E 81, 051918 (2010).
  * Cremer  _et al._ [1996] P. S. Cremer, X. Su, Y. R. Shen, and G. A. Somorjai, Ethylene hydrogenation on pt (111) monitored in situ at high pressures using sum frequency generation, Journal of the American Chemical Society 118, 2942 (1996).
  * Slepian and Sonnenblick [1965] D. Slepian and E. Sonnenblick, Eigenvalues associated with prolate spheroidal wave functions of zero order, The Bell System Technical Journal 44, 1745 (1965).
  * Moore and Cada [2004] I. C. Moore and M. Cada, Prolate spheroidal wave functions, an introduction to the slepian series and its properties, Applied and Computational Harmonic Analysis 16, 208 (2004).
  * Vahlbruch  _et al._ [2016] H. Vahlbruch, M. Mehmet, K. Danzmann, and R. Schnabel, Detection of 15 db squeezed states of light and their application for the absolute calibration of photoelectric quantum efficiency, Physical review letters 117, 110801 (2016).
  * Slepian [1964] D. Slepian, Prolate spheroidal wave functions, fourier analysis and uncertainty—iv: extensions to many dimensions; generalized prolate spheroidal functions, Bell System Technical Journal 43, 3009 (1964).
  * Shkolnisky [2007] Y. Shkolnisky, Prolate spheroidal wave functions on a disc—integration and approximation of two-dimensional bandlimited functions, Applied and Computational Harmonic Analysis 22, 235 (2007).
  * Yang  _et al._ [2002] Q. X. Yang, M. A. Lindquist, L. Shepp, C.-H. Zhang, J. Wang, and M. B. Smith, Two dimensional prolate spheroidal wave functions for mri, Journal of Magnetic Resonance 158, 43 (2002).
  * Boyer  _et al._ [2008] V. Boyer, A. M. Marino, R. C. Pooser, and P. D. Lett, Entangled images from four-wave mixing, Science 321, 544 (2008).


Navigating data-centric artificial intelligence with DC-Check: Advances, challenges, and opportunities.
IEEE Transactions on Artificial Intelligence, 2023b.

[Shwartz et al., 2022]
Daniel Shwartz, Uri Stern, and Daphna Weinshall.
The dynamic of consensus in deep networks and the identification of noisy labels.
arXiv preprint arXiv:2210.00583, 2022.

[Snoek et al., 2018]
Jasper Snoek, Alex Wiltschko, and Ali Rahimi.
Winner's curse? On pace, progress, and empirical rigor.
In International Conference on Learning Representations - Workshop, 2018.

[Song et al., 2022]
Hwanjun Song, Minseok Kim, Dongmin Park, Yooju Shin, and Jae-Gil Lee.
Learning from noisy labels with deep neural networks: A survey.
IEEE Transactions on Neural Networks and Learning Systems, 2022.

[Sorscher et al., 2022]
Ben Sorscher, Robert Geirhos, Shashank Shekhar, Surya Ganguli, and Ari S Morcos.
Beyond neural scaling laws: Beating power law scaling via data pruning.
In Advances in Neural Information Processing Systems, 2022.

[Sukhbaatar et al., 2015]
Sainbayar Sukhbaatar, Joan Bruna, Manohar Paluri, Lubomir Bourdev, and Rob Fergus.
Training convolutional networks with noisy labels.
In 3rd International Conference on Learning Representations, ICLR 2015, 2015.

[Sun et al., 2023]
Hao Sun, Alex J Chan, Nabeel Seedat, Alihan Hüyük, and Mihaela van der Schaar.
When is off-policy evaluation useful? a data-centric perspective.
arXiv preprint arXiv:2311.14110, 2023a.

[Sun et al., 2023]
Hao Sun, Boris van Breugel, Jonathan Crabbé, Nabeel Seedat, and Mihaela van der Schaar.
What is flagged in uncertainty quantification? latent density models for uncertainty categorization.
Advances in Neural Information Processing Systems, 36, 2023b.

[Swayamdipta et al., 2020]
Swabha Swayamdipta, Roy Schwartz, Nicholas Lourie, Yizhong Wang, Hannaneh Hajishirzi, Noah A Smith, and Yejin Choi.
Dataset cartography: Mapping and diagnosing datasets with training dynamics.
In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pp.  9275–9293, 2020.

[Tian et al., 2021]
Junjiao Tian, Yen-Chang Hsu, Yilin Shen, Hongxia Jin, and Zsolt Kira.
Exploring covariate and concept shift for detection and confidence calibration of out-of-distribution data.
arXiv preprint arXiv:2110.15231, 2021.

[Toneva et al., 2019]
Mariya Toneva, Alessandro Sordoni, Remi Tachet des Combes, Adam Trischler, Yoshua Bengio, and Geoffrey J Gordon.
An empirical study of example forgetting during deep neural network learning.
In International Conference on Learning Representations, 2019.

[Vanschoren et al., 2014]
Joaquin Vanschoren, Jan N Van Rijn, Bernd Bischl, and Luis Torgo.
OpenML: Networked science in machine learning.
ACM SIGKDD Explorations Newsletter, 150 (2):0 49–60, 2014.

[Wang & Jia, 2023]
Jiachen T Wang and Ruoxi Jia.
Data banzhaf: A robust data valuation framework for machine learning.
In International Conference on Artificial Intelligence and Statistics, pp.  6388–6421. PMLR, 2023.

[West, 2020]
Amanda West.
Machine Learning Models for the Future: The rise of the data scientist.
Refinitiv Labs, 2020.

[Whang et al., 2023]
Steven Euijong Whang, Yuji Roh, Hwanjun Song, and Jae-Gil Lee.
Data collection and quality challenges in deep learning: A data-centric AI perspective.
The VLDB Journal, pp.  1–23, 2023.

[Winkens et al., 2020]
Jim Winkens, Rudy Bunel, Abhijit Guha Roy, Robert Stanforth, Vivek Natarajan, Joseph R Ledsam, Patricia MacWilliams, Pushmeet Kohli, Alan Karthikesalingam, Simon Kohl, et al.
Contrastive training for improved out-of-distribution detection.
arXiv preprint arXiv:2007.05566, 2020.

[Xia et al., 2021]
Xiaobo Xia, Tongliang Liu, Bo Han, Mingming Gong, Jun Yu, Gang Niu, and Masashi Sugiyama.
Sample selection with uncertainty of losses for learning with noisy labels.
In International Conference on Learning Representations, 2021.

[Yang et al., 2023]
Jingkang Yang, Kaiyang Zhou, and Ziwei Liu.
Full-spectrum out-of-distribution detection.
International Journal of Computer Vision, pp.  1–16, 2023.

[Yuksekgonul et al., 2023]
Mert Yuksekgonul, Linjun Zhang, James Zou, and Carlos Guestrin.
Beyond confidence: Reliable models should also consider atypicality.
arXiv preprint arXiv:2305.18262, 2023.

[Zhang et al., 2022]
Shujian Zhang, Chengyue Gong, Xingchao Liu, Pengcheng He, Weizhu Chen, and Mingyuan Zhou.
ALLSH: Active learning guided by local sensitivity and hardness.
In Findings of the Association for Computational Linguistics: NAACL 2022, pp.  1328–1342, 2022.

[Zhang et al., 2021]
Yivan Zhang, Gang Niu, and Masashi Sugiyama.
Learning noise transition matrix from only noisy labels via total variation regularization.
In International Conference on Machine Learning, pp.  12501–12512. PMLR, 2021.

[Zhu et al., 2022]
Zhaowei Zhu, Jialu Wang, and Yang Liu.
Beyond images: Label noise transition matrix estimation for tasks with lower-quality features.
In International Conference on Machine Learning, pp.  27633–27653. PMLR, 2022.