# Scattering-lens based quantum imaging beyond shot noise

Dong Li  Microsystems and Terahertz Research Center, China Academy of
Engineering Physics, Chengdu Sichuan, 610200, China  Institute of Electronic
Engineering, China Academy of Engineering Physics, Mianyang Sichuan, 621999,
China  Yao Yao  Microsystems and Terahertz Research Center, China Academy of
Engineering Physics, Chengdu Sichuan, 610200, China  Institute of Electronic
Engineering, China Academy of Engineering Physics, Mianyang Sichuan, 621999,
China<EMAIL_ADDRESS>
###### Abstract

The scheme of optical imaging using scattering lens can provide a resolution
beyond the classical optical diffraction limit with a coherent-state input.
Nevertheless, due to the shot noise of the coherent state, the corresponding
signal-to-noise ratio and resolution are both still shot-noise-limited. In
order to circumvent this problem, we theoretically propose an alternative
scheme where the squeezed state (with a sub-shot noise) is considered as input
and the quantum noise is then suppressed below the shot-noise level.
Consequently, when comparing with the previous imaging scheme (using
combination of coherent state and scattering lens), our proposal is able to
achieve an enhanced signal-to-noise ratio for a given scattering lens.
Meanwhile, it is demonstrated that the resolution is also improved. We believe
that this method may afford a new way of using squeezed states and enable a
higher performance than that of using coherent state and scattering lens.

## Introduction

A significant fraction of the research activities in the field of high-
resolution imaging involves the scattering lens [1, 2, 3, 4, 5, 6, 7, 8, 9,
10]. This is due to the fact that the optical imaging system using scattering
lens could provide a better resolution than that of the conventional lens.
Different from the conventional lens consisting of an ordered structure, the
scattering lens is indeed a disordered medium [11, 12] which generally
comprises randomly distributed small particles for light scattering.
Previously, the disordered medium was deemed not suitable for the optical
imaging. Since after a beam propagates through a disordered medium, it
generally produces a speckle pattern owing to the multiple scattering, which
was supposed to deteriorate the original information carried by the incident
light. However, in recent years, it has been found that the disordered medium
has the ability to overcome the classical diffraction limit [1, 2] and can be
used for the optical imaging [3, 4], even with a better performance than that
of the conventional lens. In order to promote a more profound imagery, the
disordered medium, utilized in the high-resolution imaging system, is usually
called a scattering lens.

In a traditional imaging configuration, it is mainly formed by a set of
optical elements, in order from the object side: an input object, two
conventional lenses, and a CCD camera. In contrast to the traditional scheme,
many other methods with an extra scattering lens [1, 2, 3, 4] have been
proposed for the high-resolution imaging in the past decades. According to the
location of the scattering lens in the optical circuit, these schemes can be
roughly divided into two categories: (I) behind the object [2, 3] (i.e.
between the object and the conventional lens, corresponding to the wide-field
imaging) and (II) in front of the object [1] (in respect of the narrow-field
imaging). It is worth pointing out that in case (I), the light behind the
object would transport through the scattering lens while in case (II), the
light illuminating the object is generated from the scattering lens. Although
these two kinds of schemes show optical circuits with the different
structures, both of them can realize an image with the enhanced resolution.
This is because the scattering lens made of randomly distributed nanoparticles
could increase the effective numerical aperture [1, 2].

Figure 1: Sketch of optical imaging scheme utilizing a scattering lens with
the squeezed-state input. A scattering lens focuses a wavefront-shaped beam
[modulated by the spatial light modulator (SLM)] on a small volume. The focus
precisely illuminates a small object in the object plane. By scanning the
focus and collecting the light behind the object, the image of the object can
be obtained.

Particularly, we will concentrate on the case (II) in this work. In 2011, a
high-resolution scheme of optical imaging using a scattering lens [1] is
proposed and experimentally realized as depicted in Fig. 1. In this scheme,
the setup can be briefly described as follows: one wavefront-shaped coherent-
state beam transports through a scattering lens and then produces a nanosized
focus to illuminate a small object in the object plane. By scanning the
nanosized focus and collecting the corresponding light behind the small
object, the image of the object can be obtained. Importantly, this experiment
yields some fantastic outcomes that the scattering lens with a coherent-state
input can achieve a sub-$100$ nm resolution at visible wavelength, better than
that of the conventional lens (an optimal resolution of order of around $200$
nm at visible wavelength [1]).

Although this method provides a high resolution with a coherent-state input,
the performance of imaging is still shot-noise-limited due to the shot noise
of the coherent state (in fact, the performance of imaging is related to the
quantum noise of the focused beam illuminating the object. When the coherent
states are considered as input, the focused beam is still a coherent state
[13]). To circumvent this problem, one possible way is to reduce the shot
noise. Fortunately, the shot noise is actually a classical limit which can be
beaten by the quantum technique [14, 15, 16].

As a typical nonclassical state, the squeezed state [17, 18, 19] is of
significant importance since it possesses a quantum noise which can be below
the shot noise level (note that a coherent state corresponds to the shot
noise) [20, 21, 22]. As a consequence, the squeezed state can improve the
signal-to-noise ratio (SNR) [23, 24, 25] and has been utilized in numerous
applications, such as, quantum-enhanced magnetometer [26, 27], gravitational
wave detection [28, 29, 30], and quantum imaging using only conventional
lenses (without any scattering lens) [31, 32, 33, 34, 35, 36, 37].

In order to achieve high performance, we propose an alternative scheme where
the input is the squeezed state instead of the coherent state, as shown in
Fig. 1. In the presence of the squeezed-state input, the SNR of the optical
imaging system is analyzed and the related resolution is also investigated. In
addition, the comparison is performed between the squeezed and coherent
states. It is found that the squeezed-state input leads to both an improved
SNR and an enhanced resolution due to the suppressed quantum noise.

## Results

### Propagation of the quantized light through a scattering lens

Figure 2: Quantized light propagating through a scattering lens (a) in the
absence of wavefront shaping, (b) in the presence of wavefront shaping.
$\hat{a}_{a^{\prime}}^{\rm{in}}$ ($\hat{a}_{b^{\prime}}^{\rm{in}}$) represents
the annihilation operator of the input mode and $\hat{a}_{a}$ ($\hat{a}_{b}$)
the output mode. When the beams are injected, without the wavefront shaping in
(a), the scattering lens separates the light into different optical channels
randomly. As a result, the output presents a speckle pattern. In (b), with the
wavefront shaping, the scattering lens couples the beams into the desired
optical paths. Hence the output presents an ordered pattern. The wavefront
shaping, performed by a spatial light modulator in (b), controls the phase of
incident light.

Unlike the case of the conventional lens, the light focusing via a scattering
lens [39, 40, 41, 38, 42, 43] actually contains two indispensable processes:
(I) multiple scattering of light inside the scattering lens and (II) shaping
the wavefront of light before the light transporting through the scattering
lens. As a matter of fact, wavefront shaping is an emerging technology for
optical imaging and focusing through disordered media [42, 43, 44], by
modulating the incident wavefront, which paves a way for manipulating the
scattered light in an expected pattern. Generally, wavefront shaping can be
performed by a spatial light modulator in experiment as shown in Fig. 1. The
spatial light modulator acting as a reprogrammable matrix of pixels imprints
desired phase values on the coherent wavefront.

Let us first review the process (I). Fig. 2(a) depicts the propagation of
quantized light through a scattering lens [45, 49, 46, 47, 48], which
comprises randomly distributed small particles for light scattering. To
characterize the scattering lens, two primary factors are introduced: the
transport mean free path $l$ and the thickness $L$. If $l\ll L$, the multiple
scattering events would occur and result in a speckle pattern [50]. Hereafter
we define $s\equiv L/l$ which determines the degree of disorder.

After multiple scattering, the scattered mode $b$ [49] can be written as

$\hat{a}_{b}=\sum_{a^{\prime}}{t_{a^{\prime}b}\hat{a}_{a^{\prime}}^{\rm{in}}}+\sum_{b^{\prime}}{r_{b^{\prime}b}\hat{a}_{b^{\prime}}^{\rm{in}}},$
(1)

where $\hat{a}^{\rm{in}}_{a^{\prime}}$ ($\hat{a}^{\rm{in}}_{b^{\prime}}$)
denote the annihilation operator of the incident modes $a^{\prime}$
($b^{\prime}$) and obey the commutation relation
$[\hat{O},\hat{O}^{\dagger}]=1$
($\hat{O}=\hat{a}^{\rm{in}}_{a^{\prime}},\hat{a}^{\rm{in}}_{b^{\prime}}$). The
transmission and reflection coefficients $t_{a^{\prime}b}$ and
$r_{b^{\prime}b}$, subject to a constraint [13]
$\sum_{a^{\prime}}|t_{a^{\prime}b}|^{2}+\sum_{b^{\prime}}|r_{b^{\prime}b}|^{2}=1$,
can be approximately regarded as complex Gaussian random variables [51].
Accordingly, $t_{a^{\prime}b}=\sqrt{T_{a^{\prime}b}}e^{i\phi_{a^{\prime}b}}$
and $r_{b^{\prime}b}=\sqrt{R_{b^{\prime}b}}e^{i\phi_{b^{\prime}b}}$ where
$\phi_{a^{\prime}b}$ ($\phi_{b^{\prime}b}$) is indeed uniformly distributed in
the interval [$0,2\pi$] while $T_{a^{\prime}b}$ and $R_{b^{\prime}b}$ are the
variables of Rayleigh distribution. In addition, the ensemble-averaged
transmission and reflection coefficients are given by
$\overline{T_{a^{\prime}b}}=1/(Ms)$ and $\overline{R_{b^{\prime}b}}=(1-1/s)/M$
[52, 53], where $M$ represents the number of transmission channels and the
overline indicates the average over ensembles. It is easily seen that as the
disorder strength $s$ increases, the average transmission coefficient
$\overline{{T_{a^{\prime}b}}}$ decreases. It is worth noting that Eq. (1)
quantifies the very general input-output relation. The specific
characteristics of the multiple scattering disordered medium are represented
by the reflection and transmission coefficients. For instance,
$t_{a^{\prime}b}$ describes the coupling between the output mode $b$ and the
input mode $a^{\prime}$.

Second, consider process (II) wavefront shaping as depicted in Fig. 2(b). In
the presence of wavefront shaping, the scattered light can be directed towards
a focus in any desired output mode with a scattering lens [38].
Mathematically, the related input-output relation [13] can be characterized as

$\hat{a}_{b}^{\rm{w}}=\sum_{a^{\prime}}{|t_{a^{\prime}b}|\hat{a}_{a^{\prime}}^{\rm{in}}}+\sum_{b^{\prime}}{r_{b^{\prime}b}\hat{a}_{b^{\prime}}^{\rm{in}}},$
(2)

where the superscript $w$ denotes wavefront shaping. In contrast to Eq. (1),
the complex transmission coefficient $t_{a^{\prime}b}$ is replaced by
$|t_{a^{\prime}b}|$ in Eq. (2), which results from the fact that the phase
modulator exactly compensates the phase retardation of each transmission
channel in the scattering lens.

In experiment, it can be realized for focusing light with a scattering lens
via wavefront shaping under the current condition in laboratory nowadays,
since the setting of phase modulation has been intensively investigated in
theory and experiments over recent decades [54, 55, 56, 57, 58]. However,
comparing with previous works concentrating mainly on the enhanced intensity
of the focused beam, we will focus on the quantum fluctuation of the focused
beam.

According to Eq. (2), the expectation value of photon number operator of the
focused beam correspondingly arrives at

$\displaystyle\langle\hat{n}_{b}^{\rm{w}}\rangle$ $\displaystyle=$
$\displaystyle\langle\hat{a}^{\rm{w}\dagger}_{b}\hat{a}^{\rm{w}}_{b}\rangle$
(3) $\displaystyle=$
$\displaystyle\sum_{a^{\prime}a^{\prime\prime}}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|\langle\hat{a}^{\rm{in}\dagger}_{a^{\prime}}\hat{a}^{\rm{in}}_{a^{\prime\prime}}\rangle+\sum_{b^{\prime}b^{\prime\prime}}r_{b^{\prime}b}^{\ast}r_{b^{\prime\prime}b}\langle\hat{a}^{\rm{in}\dagger}_{b^{\prime}}\hat{a}^{\rm{in}}_{b^{\prime\prime}}\rangle$
$\displaystyle+\sum_{a^{\prime}b^{\prime}}(|t_{a^{\prime}b}|r_{b^{\prime}b}\langle\hat{a}^{\rm{in}\dagger}_{a^{\prime}}\hat{a}^{\rm{in}}_{b^{\prime}}\rangle+h.c.),$

where
$\langle\hat{n}_{k}^{\rm{in}}\rangle=\langle\hat{a}_{k}^{\rm{in}\dagger}\hat{a}_{k}^{\rm{in}}\rangle$
($k=a^{\prime},b^{\prime}$) and this expectation value is universal for any
input state.

### Signal-to-noise ratio

#### Variance of the photon number of the focused beam

Consider the squeezed states as input,
$|\Psi^{\rm{in}}\rangle=[\hat{D}(\alpha)\hat{S}(\zeta)|0\rangle]^{\otimes N}$,
with $N$ being the number of input modes,
$\hat{D}(\alpha)=e^{\alpha\hat{a}^{\dagger}-\alpha^{\ast}\hat{a}}$ the
displacement operator, and $\hat{S}(\zeta)=e^{(-\zeta\hat{a}^{\dagger
2}+\zeta^{\ast}\hat{a}^{2})/2}$ the squeezing operator (the complex number
$\alpha=|\alpha|e^{i\phi_{\alpha}}$ and the complex number
$\zeta=ge^{i\phi_{s}}$ with the real number $g$ denoting the squeezing
strength). For simplicity, we assume that the number of transmission channels
is equal to the number of input modes $M=N$. Note that in our scheme, the
input beam on the left-hand side of the scattering lens is the squeezed state
whereas the one on the right-hand side is actually the vacuum state (i.e.
$\langle\hat{n}^{\rm{in}}_{b^{\prime}}\rangle=0$).

The variance of operator $\hat{O}$ is defined as

$\displaystyle\langle(\Delta\hat{O})^{2}\rangle\equiv\langle\hat{O}^{2}\rangle-\langle\hat{O}\rangle^{2},$
(4)

where $\hat{O}=\hat{n}_{b}^{\rm{w}}$. That is to say, to obtain the variance,
it requires to compute $\langle\hat{n}_{b}^{\rm{w}}\rangle$ and
$\langle(\hat{n}_{b}^{\rm{w}})^{2}\rangle$.

The expectation value $\langle\hat{n}_{b}^{\rm{w}}\rangle$, according to Eq.
(3), can be obtained

$\displaystyle\langle\hat{n}_{b}^{\rm{w}}\rangle=\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}\sinh^{2}g+\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}||\alpha|^{2},$
(5)

where we present the derivation in Section Methods. Consider that
$|\alpha|^{2}\gg\sinh^{2}g$, the second term in Eq. (5) dominates (To the best
of our knowledge, the maximum achievable value for the squeezing parameter is
around $g\approx 1.5$ in experiment [59] and correspondingly
$\sinh^{2}g\approx 4.53$. In contrast, the order of magnitude of
$|\alpha|^{2}$ can be easily greater than that of $\sinh^{2}g$ for a bright
squeezed state in experiment [60]. Therefore, it is reasonable to assume that
$|\alpha|^{2}\gg\sinh^{2}g$.). As a result, Eq. (5) is roughly equal to

$\displaystyle\langle\hat{n}_{b}^{\rm{w}}\rangle\simeq$
$\displaystyle|\alpha|^{2}(\sum_{a^{\prime}a^{\prime\prime}}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|).$
(6)

According to the definition of variance, one can write the variance of
operator $\hat{n}^{\rm{w}}_{b}$ as

$\displaystyle\langle(\Delta\hat{n}^{\rm{w}}_{b})^{2}\rangle\equiv\langle(\hat{n}^{\rm{w}}_{b})^{2}\rangle-\langle\hat{n}^{\rm{w}}_{b}\rangle^{2}.$
(7)

The corresponding variance of photon number is found to be

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle$ $\displaystyle=$
$\displaystyle\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}T_{a^{\prime}b}T_{a^{\prime\prime}b}2\sinh^{2}g\cosh^{2}g+\sum_{a^{\prime}=1}^{M}\sum_{b^{\prime}=1}^{M}T_{a^{\prime}b}R_{b^{\prime}b}\sinh^{2}g$
(8)
$\displaystyle+|\alpha|^{2}(\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|)[1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}(1-e^{-2g})],$

where we have set $\phi_{\alpha}=\phi_{s}=0$ and the detailed derivation is
shown in Section Methods. Particularly, if $|\alpha|^{2}$ is sufficiently
large for the third term of Eq. (8) to dominate, the variance can be
simplified to

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle\simeq\langle\hat{n}_{b}^{\rm{w}}\rangle[1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}(1-e^{-2g})].$
(9)

From Eq. (9), it is easy to find that the variance of photon number depends
upon the value of the squeezing parameter $g$. Thus it pays to consider the
two opposing limits $g=0$ and $g\approx 1.5>0$ (in experiment, $g\approx 1.5$
has been reported [59]). First, consider $g=0$ (i.e. the coherent-state
input), one can rewrite Eq. (9) as

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle=\langle\hat{n}_{b}^{\rm{w}}\rangle,$
(10)

which is exactly as expected. Note that when the input is a coherent state,
the focused mode is still a coherent state due to the linear optical process
in the scattering lens (the variance of photon number of a coherent state is
equal to its mean photon number). Next, we consider the opposing limit
$g\approx 1.5$ (i.e. the squeezed-state input and $e^{-2g}\to 0$). In this
situation, Eq. (9) can be then approximately reduced to

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle\simeq\langle\hat{n}_{b}^{\rm{w}}\rangle[1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}].$
(11)

Comparing Eqs. (10) and (11), one can find that the squeezed-state input leads
to the reduction of quantum fluctuation of the focused beam. Moreover, from
Eq. (11), it is easily seen that the degree of this reduction is related to a
specific set of transmission coefficients of the scattering lens.

For convenience, the Fano factor is introduced and defined as

$\displaystyle
F\equiv\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle/\langle\hat{n}_{b}^{\rm{w}}\rangle.$
(12)

When the input is a coherent state (i.e. $g=0$), by plugging Eq. (11) into Eq.
(12), the Fano factor is found to be $F=1$. In contrast to the coherent-state
input, we consider that the squeezed state is injected (with a large $g$, i.e.
$e^{-2g}\to 0$). By inserting Eq. (11) into Eq. (12), it is easy to verify
that the Fano factor of the focused beam is roughly equal to

$\displaystyle F\simeq 1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b},$ (13)

which is smaller than that of the coherent state. This result reveals that the
squeezed-state input can achieve a focused beam with the lower quantum noise
than that of the coherent state.

By averaging over all disorder ensembles, according to Eq. (13), the average
Fano factor is found to be

$\displaystyle\overline{F}\simeq 1-\frac{1}{s},$ (14)

where the overline means the average and
$\overline{\sum_{a^{\prime}}T_{a^{\prime}b}}=N\overline{T_{a^{\prime}b}}=1/s$
is used [49]. From Eq. (14), one can see that with the increase of $s$, the
average Fano factor increases monotonically which indicates that the output
quantum noise increases. It is worth pointing out that this conclusion is
suitable for the regime of the strong squeezing strength ($e^{-2g}\to 0$).

Figure 3: The ratio of the average SNR $\overline{R}$ to the mean photon
number $\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}$ as a function of (a)
the squeezing parameter $g$ and (b) the disorder strength $s$ of the
scattering lens. The input is the product state of $N$ squeezed states.
Parameters: $|\alpha|^{2}=10000$.

#### Signal-to-noise ratio

In our scheme, The light illuminating the object for imaging is the focused
beam as shown in Fig. 1. Accordingly, the focused beam determines the SNR of
the imaging system [36] which can be defined as

$\displaystyle
R\equiv\frac{\langle\hat{n}^{\rm{w}}_{b}\rangle^{2}}{\langle(\Delta\hat{n}^{\rm{w}}_{b})^{2}\rangle}=\frac{\langle\hat{n}^{\rm{w}}_{b}\rangle}{F}=\frac{\langle\hat{n}^{\rm{w}}_{b}\rangle}{1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}(1-e^{-2g})},$
(15)

where $\hat{n}^{\rm{w}}_{b}$ represents the photon number operator of the
focused mode, $F$ denotes the corresponding Fano factor and $T_{a^{\prime}b}$
characterizes the coupling between the output mode $b$ and the input mode
$a^{\prime}$. When $g=0$ (i.e. the coherent-state input and $F=1$), the SNR is
naturally given by

$\displaystyle R=\langle\hat{n}^{\rm{w}}_{b}\rangle,$ (16)

and the corresponding average SNR can be written as
$\overline{R}=\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}.$ When $e^{-2g}\to
0$ (i.e. the squeezed-state input), the SNR can be written as

$\displaystyle
R=\frac{\langle\hat{n}^{\rm{w}}_{b}\rangle}{F}\simeq\frac{\langle\hat{n}^{\rm{w}}_{b}\rangle}{1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}}.$
(17)

By averaging over all the disorder ensembles, the ensemble-averaged SNR is
obtained

$\displaystyle\overline{R}=\frac{\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}}{\overline{F}}\simeq\frac{\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}}{1-1/s},$
(18)

where $\overline{\sum_{a^{\prime}}T_{a^{\prime}b}}=1/s$ has been used and $s$
denotes the disorder parameter. From Eq. (18), it is easy to check that the
average SNR is improved as the decrease of disorder strength $s$ for a given
$\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}$.

Figure 4: Numerical simulation of Fano factor of the focused mode. Each point
corresponds to one Fano factor. The numerical experiments have been repeated
for 1000 times. In each simulation circle, the scattering matrix is generated
randomly. The input beam is the photon-number squeezed state
$|\Psi^{\rm{in}}\rangle=[\hat{D}(\alpha)\hat{S}(\zeta)|0\rangle]^{\otimes N}$.
Parameters: $|\alpha|^{2}=10000$, $N=50$.

Based on Eq. (15), one can obtain the universal average SNR as
$\overline{R}=\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}/[1-(1-e^{-2g})/s]$.
Figures 3(a) and 3(b) plot the average SNRs as a function of the squeezing
parameter $g$ and the disorder strength $s$, respectively. It is obvious that
the squeezed state always has a better average SNR than that of the coherent
state. In other words, the squeezed state can improve the SNR. In Fig. 3(a),
with the increase of $g$, the average SNR increases which yields that the
larger the squeezing parameter is, the better the average SNR is. However, in
Fig. 3(b), with the increase of $s$, this advantage resulting from the
squeezed states would weaken due to the increased quantum noise induced from
reflected modes (the vacuum state with the shot noise). It is worth noting
that although the average SNR becomes worse, it is still better than that of
the coherent state. To be more clear, we plot the numerical results of Fano
factor for four cases (($g=1.5,s=2$), ($g=1.5,s=6$), ($g=1,s=2$), and
($g=1,s=6$)) in Fig. 4. It is easily seen that in Fig. 4 the Fano factors are
always smaller than one, which indicates the focused modes with a sub-shot
noise under these four situations.

### Resolution

In addition to the SNR, the resolution is another significant criterion of
performance evaluation for an optical imaging system. In this section, we
would like to investigate how the squeezed states affect the resolution of
imaging system. As a matter of fact, in 2005, Beskrovnyy and Kolobov [36]
developed a model that allows for analysis of the effects of squeezed state on
the resolution. This model could be adopted in our circumstance.

In essence, the analysis model is composed by two core parts: (i) the first
step is to express the point spread function (PSF) [36] in terms of prolate
spheroidal functions [61] within a quantum mechanical framework; (ii) the
second one is to work out the “cutoff” PSF limited by the quantum noise (under
the critical condition of the SNR decaying to one). Once the “cutoff” PSF is
obtained, it is easy to find the corresponding resolution [36] since the
resolution is determined by the half-width of the “cutoff” PSF.

For the sake of clarity, we will not present the analysis model in this
manuscript, since the detailed derivation is of extreme complexity involving a
huge number of mathematical expressions (i.e. prolate spheroidal functions
[61]). Readers who are interested in more details about this model may refer
to the Beskrovnyy and Kolobov paper [36]. Based on this analysis model, the
quantity that we consider is the super-resolution factor $J$ [36] which will
be briefly introduced in the following section.

#### Definition of the super-resolution factor

In this section, we briefly introduce the definition of the super-resolution
factor [36]. The super-resolution factor, proposed by Beskrovnyy and Kolobov
[36], is the ratio of the half width of PSF in quantum optics theory to the
one in classical optics theory.

In classical optics theory, the PSF of the traditional imaging system [36] is
given by

$\displaystyle
h(z-z^{\prime})=\frac{\sin[c(z-z^{\prime})]}{\pi(z-z^{\prime})},$ (19)

where $c$ denotes the spatial transmission bandwidth [62] and $z$
($z^{\prime}$) is the spatial coordinate.

By contrast, within the quantum mechanical framework, the PSF of the modified
imaging system [36] is found to be

$\displaystyle
h^{(r)}(z,z^{\prime})=\sum_{k=0}^{Q-1}\phi_{k}(z)\phi_{k}(z^{\prime}),$ (20)

where the detailed derivation is shown in Section Methods, the real number $Q$
is determined by the SNR $R$ and $\phi_{k}(z)$ is the prolate spheroidal
function [36]. In fact, the set of functions $\phi_{k}(z)$ [$k=1,2,3...$]
constitutes a complete orthonormal basis [61].

To evaluate the resolution, the super-resolution factor $J$ [36] is introduced
and defined as

$\displaystyle J=\frac{W}{W_{Q}},$ (21)

where we present the numerical analysis method to obtain $J$ in Section
Methods, for brevity. $W$ [$W_{Q}$] is the half-width of the PSF
$h(z-z^{\prime})$ [$h^{(r)}(z,z^{\prime})$].

#### Resolution

Figure 5: Super-resolution factor $J$ as a function of the mean number of
photon in the object plane. The red-solid, green-dashed, blue-dashed, orange-
dotted curves denote the squeezed state ($s=2$), squeezed state ($s=4$),
squeezed state ($s=6$), squeezed state ($s=8$), and coherent state,
respectively. Parameters: $g=1.5$.

Let us consider a typical condition for estimating the value of mean photon
number in the focused beam illuminating the object in experiment. For
instance, if it is given an incident beam with $\lambda=694$ nm, optical power
of $1$ mW and observation time of $1$ ms, the focused beam then has an optical
power of $0.01$ mW (where it has been assumed that there is $1/100$ of total
incident power in the focused mode [38]). Finally, the mean photon number of
the focused beam can be worked out as
$\langle\hat{n}^{\rm{w}}_{b}\rangle\simeq 3.47\times 10^{10}$.

Based on the order of value of $\langle\hat{n}^{\rm{w}}_{b}\rangle$ discussed
above, we plot the super-resolution factors $J$ as a function of mean photon
number $\langle\hat{n}^{\rm{w}}_{b}\rangle$ as shown in Fig. 5. In particular,
the comparison between the squeezed-state input and the coherent-state one is
performed. From Fig. 5, one can easily see that the squeezed state has a
better resolution than that of the coherent state for a given intensity of the
focused mode. In other words, the squeezed states could enhance the resolution
owing to the sub-shot noise. We also investigate the effects of disorder
strength $s$ on the resolution in Fig. 5. It can be seen that with the
decrease of $s$, the resolution with squeezed state is improved for a fixed
intensity $\overline{\langle\hat{n}^{\rm{w}}_{b}\rangle}$. This is due to the
fact that the decrease of $s$ leads to the decrease of the average noise of
the focused mode (Eq. (14)).

It is worth noting that the increase of the super-resolution factor is not
dramatic as depicted in Fig. 5. Nevertheless, the squeezed state is still
deserved to consider as input since any tiny improvement in high-precision
measurement is very valuable.

## Discussion

### Comparison between our proposal and the previous imaging scheme

Our proposal is inspired by the scheme using coherent state with a scattering
lens in the previous work [1]. Therefore it is important to compare our
proposal with the one in the reference [1]. To be fair, we would like to
consider the same situation (such as, same disorder parameter $s$). Before
comparing these two schemes, let us estimate the value of $s$ in the former
work [1]. From the reference [1], the length of disordered medium is given by
$L=2.88$ ${\rm{\mu m}}$. The mean free path is found to be $l=0.47\pm 0.05$
${\rm{\mu m}}$ according to the previous paper [63] of A. Lagendijk (one of
the authors of the reference [1]). As a result, $s=L/l\approx 6$ in the former
work [1]. Accordingly, we also consider $s=6$ in our imaging scheme (see Figs.
3 and 5).

In Fig. 3, it is easily found that the scheme using squeezed state (dashed-
blue curve, $s=6$) has an enhanced SNR than that with coherent state (solid-
gray curve), which results in an improved resolution as depicted in Fig. 5.
That is to say, our scheme has a resolution beyond that in the reference [1]
when $s=6$ are same. Actually, this conclusion can be extended to the case of
almost any value of $s$ (including not only weak scattering lens but also
strong scattering lens).

It is worth noting that our claim on the improvement of imaging resolution is
valid when comparing our proposal and the one in the reference [1] under the
same conditions (same parameter $s$). In fact, our analysis model does not
involve the effect of $s$ on the increase of angular bandwidth of scattering
lens $or$ the spatial size of the focused spot when analyzing the super-
resolution factor $J$. This is still an open question. We would like to
investigate this effect in the future.

### SNR in the case of $N<M$

Since the number of input squeezed-state modes $N$ may not be exactly equal to
the number of transmission channels $M$ in experiment, it is worth discussing
how this affects the SNR in imaging. In the situation of $N<M$, one can
calculate the mean photon number

$\displaystyle\langle\hat{n}_{b}^{\rm{w}}\rangle_{N<M}$
$\displaystyle=\langle\hat{a}^{\rm{w}\dagger}_{b}\hat{a}^{\rm{w}}_{b}\rangle$
$\displaystyle=\sum_{a^{\prime}=1}^{N}T_{a^{\prime}b}\sinh^{2}g+\sum_{a^{\prime}=1}^{N}\sum_{a^{\prime\prime}=1}^{N}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}||\alpha|^{2}.$
(22)

Note that Eqs. (5) and (22) are different where $N<M$. Similar to Eq. (8), the
corresponding variance arrives at

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle_{N<M}=$
$\displaystyle\sum_{a^{\prime}=1}^{N}\sum_{a^{\prime\prime}=1}^{N}T_{a^{\prime}b}T_{a^{\prime\prime}b}2\cosh^{2}g\sinh^{2}g+\sum_{a^{\prime}=1}^{N}\sum_{b^{\prime}=1}^{M}T_{a^{\prime}b}R_{b^{\prime}b}\sinh^{2}g+\sum_{a^{\prime}=1}^{N}\sum_{a^{\prime\prime}=N+1}^{M}T_{a^{\prime}b}T_{a^{\prime\prime}b}\sinh^{2}g$
$\displaystyle+|\alpha|^{2}\sum_{a^{\prime}=1}^{N}\sum_{a^{\prime\prime}=1}^{N}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|[1-\sum_{a^{\prime}=1}^{N}T_{a^{\prime}b}(1-e^{-2g}),$
(23)

where we have set $\phi_{\alpha}=\phi_{s}=0$. Compared to Eq. (8), Eq. (23)
possesses an additional term of
$\sum_{a^{\prime}=1}^{N}\sum_{a^{\prime\prime}=N+1}^{M}T_{a^{\prime}b}T_{a^{\prime\prime}b}\sinh^{2}g$.
This term results from the interference between the squeezed state and the
vacuum state from the empty input port.

Figure 6: The ratio of the average signal-to-noise ratio to the mean photon
number $\overline{R_{N<M}}/\overline{\langle\hat{n}_{b}^{\rm{w}}\rangle}$ as a
function of $N/M$. $N$: number of input modes, $M$: number of transmission
channels. It is found that the increasing of $N/M$ leads to the increase in
the ratio $\overline{R_{N<M}}/\overline{\langle\hat{n}_{b}^{\rm{w}}\rangle}$.
Parameters: $|\alpha|^{2}=10000,M=50$.

By combining Eqs. (15), (22), and (23), one can obtain the signal-to-noise
ratio $R_{N<M}$. The ratio
$\overline{R_{N<M}}/\overline{\langle\hat{n}_{b}^{\rm{w}}\rangle}$ could be
numerically analyzed. Fig. (6) plots the
$\overline{R_{N<M}}/\overline{\langle\hat{n}_{b}^{\rm{w}}\rangle}$ as a
function of $N/M$. From Fig. (6), it is easily found that the ratio
$\overline{R_{N<M}}/\overline{\langle\hat{n}_{b}^{\rm{w}}\rangle}$ increases
as the increasing of $N/M$. The maximum
$\overline{R_{N<M}}/\overline{\langle\hat{n}_{b}^{\rm{w}}\rangle}$ can be
achieved if $N/M=1$.

### The output noise in the presence of $|\alpha|^{2}$ in various regimes

In contrast to the regime of large $|\alpha|^{2}$ in the previous section, we
take into account a full range of $|\alpha|^{2}$. From Eqs. (5), (8) and (12),
it is easy to compute the universal Fano factor $F^{\rm{univ}}$. Accordingly,
the average $\overline{F^{\rm{univ}}}$ could be obtained by numerical
simulation. Fig. 7 presents the Fano factor $\overline{F^{\rm{univ}}}$ as a
function of $|\alpha|^{2}/(|\alpha|^{2}+\sinh^{2}g)$. It can be seen that with
the increase of $|\alpha|^{2}/(|\alpha|^{2}+\sinh^{2}g)$, the Fano factor
$\overline{F^{\rm{univ}}}$ decreases. One can obtain the optimal
$\overline{F^{\rm{univ}}}$ when $|\alpha|^{2}/(|\alpha|^{2}+\sinh^{2}g)\simeq
1$.

It is worth pointing out that when $|\alpha|^{2}=450$, the ratio
$|\alpha|^{2}/(|\alpha|^{2}+\sinh^{2}g)\simeq 0.99$ with $g=1.5$
($\sinh^{2}g\simeq 4.53$) which corresponds roughly to the optimal
$\overline{F^{\rm{univ}}}$. Therefore, $|\alpha|^{2}=450$ could be considered
as a strong coherent intensity when $g=1.5$, although it is not a strong
intensity in experiment. Since $|\alpha|^{2}>450$ is easily achieved in
experiment, it is reasonable to assume that $|\alpha|^{2}$ is in the regime of
large value.

Figure 7: The average Fano factor $\overline{F^{\rm{univ}}}$ of the focused
mode as a function of $|\alpha|^{2}/(|\alpha|^{2}+\sinh^{2}g)$, where
$|\alpha|^{2}$ means the “coherent” intensity of input beam and
$|\alpha|^{2}+\sinh^{2}g$ indicates the corresponding total intensity.

### Effect of photon loss on the SNR

Figure 8: The effect of photon loss on the average signal-to-noise ratio.
Parameters: $|\alpha|^{2}=10000$, (a) $s=2$, (b) $s=4$, (c) $s=6$, (d) $s=8$.

In practical experiment, the photon loss is inevitable. Therefore, we consider
the effect of photon loss of focused mode on the performance of imaging (SNR,
particularly). To mimic the photon loss, we assume that the focused mode
passes through a fictional two-port beam splitter with vacuum state injected
in the other port. After propagating though the fictional beam splitter, the
final output state can be characterized by

$\displaystyle\hat{a}^{\rm{w}\dagger}_{b,L}$
$\displaystyle=p^{\ast}\hat{a}^{\rm{w}\dagger}_{b}+q^{\ast}\hat{a}^{\rm{V}\dagger},$
$\displaystyle\hat{a}^{\rm{w}}_{b,L}$
$\displaystyle=p\hat{a}^{\rm{w}}_{b}+q\hat{a}^{\rm{V}},$ (24)

where complex number $p$ and $q$ denotes transmission and reflection
coefficients ($|q|^{2}$ indicates the photon-loss rate) and $\hat{a}^{\rm{V}}$
represents the annihilation operator of vacuum state. The mean photon number
of focus state with loss is found to be

$\displaystyle\langle\hat{n}^{\rm{w}}_{b,L}\rangle=|p|^{2}\langle\hat{n}^{\rm{w}}_{b}\rangle.$
(25)

The average SNR is given by

$\displaystyle\overline{R_{L}}=\overline{\langle\hat{n}_{b,L}^{\rm{w}}\rangle}/\overline{F_{L}},$
(26)

where $\overline{F_{L}}$ can be expressed as

$\displaystyle\overline{F_{L}}=|p|^{2}\overline{F}+|q|^{2},$ (27)

with $\overline{F}$ being the Fano factor without loss (For clarity, the
detailed derivation is present in Section Methods). Fig. 8 depicts the average
SNR as a function of loss $|q|^{2}$ with various $g$. It is shown that
although the averaged SNR degrades gradually as the increase of loss, the
average SNR is still beyond the shot-noise level with small or moderate loss.

## Conclusions

Prior studies have shown that a high-resolution imaging can be achieved with a
scattering lens via coherent states. This is owing to the fact that the
scattering lens has a large numerical aperture and is able to focus the
coherent-state light to a tighter spot than the diffraction limit of that
conventional lens. Nevertheless, this is a pure classical technique without
the help of any quantum technique. In the point of view of quantum optics, the
squeezed state, as a nonclassical state, can enable a high-precision
measurement. Therefore, in the current literature, we propose an alternative
imaging scheme in which the squeezed states are considered as input instead of
the coherent states.

Compared with the previous scheme, this new method takes advantage of not only
the large numerical aperture of the scattering lens but also the sub-shot
noise of the squeezed states. Consequently, our scheme with squeezed states
has a better performance than that of the traditional imaging scheme. On the
other hand, our scheme establishes a superiority that maintaining the same
performance (i.e. resolution) requires a reduced photon number of the input
states. This can eliminate the disadvantages due to the strong coherent-state
input.

In summary, the effect of the squeezed-state input on the performance of
imaging in the optical system with a scattering lens is investigated. It is
clarified that the squeezed state leads to an improved signal-to-noise ratio
in imaging compared with the coherent state due to the suppressed quantum
noise. Moreover, it is found that the squeezed state also achieves an enhanced
resolution in contrast to the coherent state. Therefore, our results may pave
a new way to realize an image with a high resolution using scattering lens
with the squeezed-state input.

## Methods

### Variance of photon number of the focused mode

Consider the squeezed states as input,
$|\Psi^{\rm{in}}\rangle=[\hat{D}(\alpha)\hat{S}(\zeta)|0\rangle]^{\otimes N}$,
with $N$ being the number of input modes,
$\hat{D}(\alpha)=e^{\alpha\hat{a}^{\dagger}-\alpha^{\ast}\hat{a}}$ the
displacement operator, and $\hat{S}(\zeta)=e^{(-\zeta\hat{a}^{\dagger
2}+\zeta^{\ast}\hat{a}^{2})/2}$ the squeezing operator (the complex number
$\alpha=|\alpha|e^{i\phi_{\alpha}}$ is the amplitude of the beam and the
complex number $\zeta=ge^{i\phi_{s}}$ denotes the squeezing parameter, real
number $g$ being the squeezing strength). The number of transmission channels
is denoted by $M$. For simplicity, we assume that $M=N$.

For the convenience of calculation, one can rewrite the input-output relation
of a disordered medium as

$\displaystyle\hat{a}_{b}^{\rm{w}}=\sum_{a^{\prime}}|t_{a^{\prime}b}|(\hat{a}^{\rm{SV}}_{a^{\prime}}+\alpha)+\sum_{b^{\prime}}r_{b^{\prime}b}\hat{a}_{b^{\prime}}^{\rm{in}},$
(28)

where the complex number $\alpha$ denotes the coherent amplitude of the input
field and the operator $\hat{a}^{\rm{SV}}_{a^{\prime}}$ accounts for the
quantum fluctuation ($\langle\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle=0$,
$\langle
0|\hat{S}^{\dagger}(\alpha)\hat{D}^{\dagger}(\alpha)\hat{a}^{\rm{in}}_{a^{\prime}}\hat{D}(\alpha)\hat{S}(\zeta)|0\rangle\to\langle
0|\hat{S}^{\dagger}(\alpha)(\hat{a}^{\rm{SV}}_{a^{\prime}}+\alpha)\hat{S}(\zeta)|0\rangle$).

From Eq. (28), it is easy to obtain the photon number operator

$\displaystyle\hat{n}_{b}^{\rm{w}}=$
$\displaystyle\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}+\sum_{b^{\prime}=1}^{M}R_{a^{\prime}b}\hat{a}^{\rm{in}\dagger}_{b^{\prime}}\hat{a}_{b^{\prime}}^{\rm{in}}+[\sum_{a^{\prime}=1}^{M}\sum_{b^{\prime}=1}^{M}|t_{a^{\prime}b}|r_{b^{\prime}b}\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}_{b^{\prime}}^{\rm{in}}+H.c.]$
(29) $\displaystyle+\sum_{a^{\prime}=1}^{M}\sum_{(a^{\prime\prime}\neq
a^{\prime},a^{\prime\prime}=1)}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}+\sum_{b^{\prime}=1}^{M}\sum_{(b^{\prime\prime}\neq
b^{\prime},b^{\prime\prime}=1)}^{M}r^{\ast}_{b^{\prime}b}r_{b^{\prime\prime}b}\hat{a}^{\rm{in}\dagger}_{b^{\prime}}\hat{a}_{b^{\prime\prime}}^{\rm{in}}$
$\displaystyle+[\sum_{a^{\prime}=1}^{M}|t_{a^{\prime}b}|\alpha^{\ast}(\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime\prime}b}|\hat{a}^{\rm{SV}}_{a^{\prime\prime}}+\sum_{b^{\prime}=1}^{M}r_{b^{\prime}b}\hat{a}_{b^{\prime}}^{\rm{in}})+H.c.]+|\alpha|^{2}\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|.$

The expectation value of $\hat{n}_{b}^{\rm{w}}$ is then obtained

$\displaystyle\langle\hat{n}_{b}^{\rm{w}}\rangle=$
$\displaystyle\sum_{a^{\prime}=1}^{M}|t_{a^{\prime}b}|^{2}\langle\hat{n}^{\rm{SV}}_{a^{\prime}}\rangle+|\alpha|^{2}\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|,$
(30)

where
$\langle\hat{n}^{\rm{SV}}_{a^{\prime}}\rangle=\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle=\sinh^{2}g$
and we have used $\langle\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle=0$ and
$\langle\hat{a}_{b^{\prime}}\rangle=0$. Consider that
$|\alpha|^{2}\gg\sinh^{2}g$, the second term in Eq. (30) dominates. As a
result, Eq. (30) is roughly equal to

$\displaystyle\langle\hat{n}_{b}^{\rm{w}}\rangle\simeq$
$\displaystyle|\alpha|^{2}\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|.$
(31)

According to the definition of variance
$\langle(\Delta\hat{n}^{\rm{w}}_{b})^{2}\rangle\equiv\langle(\hat{n}^{\rm{w}}_{b})^{2}\rangle-\langle\hat{n}^{\rm{w}}_{b}\rangle^{2}$,
one can obtain the variance of photon number

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle=$
$\displaystyle\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}^{2}\langle(\Delta\hat{n}_{a^{\prime}}^{\rm{SV}})^{2}\rangle+\sum_{a^{\prime}=1}^{M}\sum_{b^{\prime}=1}^{M}T_{a^{\prime}b}R_{b^{\prime}b}\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle$
(32) $\displaystyle+\sum_{a^{\prime}=1}^{M}\sum_{(a^{\prime\prime}\neq
a^{\prime},a^{\prime\prime}=1)}^{M}T_{a^{\prime}b}T_{a^{\prime\prime}b}(\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\hat{a}^{\rm{SV}\dagger}_{a^{\prime\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}+\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle+\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}\rangle)$
$\displaystyle+|\alpha|^{2}\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|[\sum_{a^{\prime\prime\prime}=1}^{M}T_{a^{\prime\prime\prime}b}(2\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime\prime\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime\prime}}\rangle+1)$
$\displaystyle+\sum_{b^{\prime}=1}^{M}R_{b^{\prime}b}+(\sum_{a^{\prime\prime\prime}=1}^{M}|T_{a^{\prime\prime\prime}b}|\langle(\hat{a}^{\rm{SV}}_{a^{\prime\prime\prime}})^{2}\rangle+h.c.)],$

where
$\langle(\Delta\hat{n}_{a^{\prime}}^{\rm{SV}})^{2}\rangle=\langle(\hat{n}^{\rm{SV}}_{a^{\prime}})^{2}\rangle-\langle\hat{n}^{\rm{SV}}_{a^{\prime}}\rangle^{2}$
and we have set $\phi_{\alpha}=0$. Eq. (32) could be further simplified to

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle=$
$\displaystyle\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}T_{a^{\prime}b}T_{a^{\prime\prime}b}2\cosh^{2}g\sinh^{2}g+\sum_{a^{\prime}=1}^{M}\sum_{b^{\prime}=1}^{M}T_{a^{\prime}b}R_{b^{\prime}b}\sinh^{2}g$
(33)
$\displaystyle+|\alpha|^{2}\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|[1-\sum_{a^{\prime\prime\prime}=1}^{M}T_{a^{\prime\prime\prime}b}(1-e^{-2g})],$

where we have utilized
$\sum_{b^{\prime}}R_{b^{\prime}b}=1-\sum_{a^{\prime}}T_{a^{\prime}b}$,
$\langle\hat{n}_{a^{\prime}}^{\rm{SV}}\rangle=\sinh^{2}g$,
$\langle(\Delta\hat{n}_{a^{\prime}}^{\rm{SV}})^{2}\rangle=2\cosh^{2}g\sinh^{2}g$,
$\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\hat{a}^{\rm{SV}\dagger}_{a^{\prime\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}+\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle=\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}\hat{a}^{\rm{SV}}_{a^{\prime\prime}}\rangle=\cosh^{2}g\sinh^{2}g$,
$2\langle\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}\hat{a}^{\rm{SV}}_{a^{\prime}}\rangle+1+\langle(\hat{a}^{\rm{SV}}_{a^{\prime}})^{2}\rangle+\langle(\hat{a}^{\rm{SV}\dagger}_{a^{\prime}})^{2}\rangle=e^{-2g}$,
and set $\phi_{s}=0$. Particularly, if $|\alpha|^{2}$ is sufficiently large
for the fourth term of Eq. (33) to dominate, the variance could be reduced to

$\displaystyle\langle(\Delta\hat{n}_{b}^{\rm{w}})^{2}\rangle\simeq\langle\hat{n}_{b}^{\rm{w}}\rangle[1-\sum_{a^{\prime}=1}^{M}T_{a^{\prime}b}(1-e^{-2g})],$
(34)

where
$\langle\hat{n}_{b}^{\rm{w}}\rangle\simeq|\alpha|^{2}\sum_{a^{\prime}=1}^{M}\sum_{a^{\prime\prime}=1}^{M}|t_{a^{\prime}b}||t_{a^{\prime\prime}b}|$
is used.

### Quantum theory of optical imaging

#### Input-output relation corresponding to a spatial Fourier transform

Recall that in our imaging scheme, with the help of the squeezed states, the
focused beam has a sub-shot noise. To uncover the role of the suppressed
quantum fluctuation in the resolution of Fourier microscopy, it requires to
analyze the configuration by using the quantum mechanical language.

Figure 9: The traditional imaging system, including the object field $a(z)$,
the first lens ($L_{1}$), the field in the pupil plane $f(\xi)$, the second
lens ($L_{2}$), and the field in the image plane $e(z)$. $X$ denotes the width
of the object field and $d$ represents the finite size of pupil. For
convenience, the dimensionless spatial coordinates $z=2x/X$ in the object and
image planes and $\xi=2y/d$ in the pupil plane are introduced.

In quantum optics, the object field in Fig. 9 is described by the operator
$\hat{a}(z)$ [$\hat{a}^{\dagger}(z)$] instead of the complex amplitude $a(z)$
in classical optics. Similarly, the field in the pupil plane is characterized
by the operator $\hat{f}(\xi)$. The operators satisfy the standard commutation
relation, $[\hat{a}(z),\hat{a}^{\dagger}(z^{\prime})]=\delta(z-z^{\prime})$
and $[\hat{f}(\xi),\hat{f}^{\dagger}(\xi^{\prime})]=\delta(\xi-\xi^{\prime})$.
Since the lens $L_{1}$ performs a spatial Fourier transform between the object
and pupil fields, the relation between $\hat{a}(z)$ and $\hat{f}(\xi)$ reads
as follows [36]

$\displaystyle\hat{f}(\xi)=\sqrt{\frac{c}{2\pi}}\int^{\infty}_{-\infty}\hat{a}(z)e^{-icz\xi}dz,$
(35)

where $c$ denotes the spatial transmission bandwidth of the optical system.

#### Input-output relation in terms of the prolate spheroidal functions

In terms of the prolate spheroidal functions, the operators $\hat{a}(z)$ and
$\hat{f}(\xi)$ can be expressed as

$\displaystyle\hat{a}(z)=\sum_{k=0}^{\infty}\hat{a}_{k}\phi_{k}(z)+\sum_{k=0}^{\infty}\hat{b}_{k}\chi_{k}(z),$
(36)
$\displaystyle\hat{f}(\xi)=\sum_{k=0}^{\infty}\hat{f}_{k}\phi_{k}(\xi)+\sum_{k=0}^{\infty}\hat{g}_{k}\chi_{k}(\xi),$
(37)

where

$\displaystyle~\phi_{k}(z)=~\begin{cases}~~\frac{1}{\sqrt{\lambda_{k}}}\psi_{k}(z),&|z|\leq
1,\\\ ~~0,&|z|>1,\end{cases}\quad\quad~\chi_{k}(z)=~\begin{cases}~~0,&|z|\leq
1,\\\ ~~\frac{1}{\sqrt{1-\lambda_{k}}}\psi_{k}(z),&|z|>1,\end{cases}$ (38)

$\hat{a}_{k}$ ($\hat{b}_{k}$) and $\hat{f}_{k}$ ($\hat{g}_{k}$) denote the
annihilation operators of the prolate mode $\phi_{k}$ ($\chi_{k}$) and can be
obtained by

$\displaystyle\hat{a}_{k}=\int_{-\infty}^{\infty}\hat{a}(z)\phi_{k}(z)dz,$
(39)
$\displaystyle\hat{b}_{k}=\int_{-\infty}^{\infty}\hat{a}(z)\chi_{k}(z)dz,$
(40)
$\displaystyle\hat{f}_{k}=\int_{-\infty}^{\infty}\hat{f}(\xi)\phi_{k}(\xi)dz,$
(41)
$\displaystyle\hat{g}_{k}=\int_{-\infty}^{\infty}\hat{f}(\xi)\chi_{k}(\xi)dz.$
(42)

The operators ($\hat{a}_{k}$, $\hat{b}_{k}$, $\hat{f}_{k}$, and $\hat{g}_{k}$)
obey the standard commutation relation,
$[\hat{O}_{k},\hat{O}_{k^{\prime}}^{\dagger}]=\delta(k-k^{\prime})$
($O=a,b,g,k$). It is worthy noting that $\phi_{k}(z)$ lies within the region
of $|z|<1$ while $\chi_{k}(z)$ is distributed over the area of $|z|>1$.

Based on the properties of the prolate spheroidal functions [64], one has

$\displaystyle\int^{1}_{-1}\phi_{k}(z)e^{-icz\xi}dz=(-i)^{k}\sqrt{\frac{2\pi}{c}}\psi_{k}(\xi),$
(43)
$\displaystyle\int^{\infty}_{-\infty}\psi_{k}(z)e^{-icz\xi}dz=(-i)^{k}\sqrt{\frac{2\pi}{c}}\phi_{k}(\xi).$
(44)

By substituting Eqs. (36), (37), (43), and (44) into (35), it is easy to
obtain the relations between the photon annihilation operators of the prolate
modes in the object and pupil planes

$\displaystyle\hat{f}_{k}=(-i)^{k}(\sqrt{\lambda_{k}}\hat{a}_{k}+\sqrt{1-\lambda_{k}})\hat{b}_{k},$
(45)
$\displaystyle\hat{g}_{k}=(-i)^{k}(\sqrt{1-\lambda_{k}}\hat{a}_{k}-\sqrt{\lambda_{k}})\hat{b}_{k},$
(46)

where Eqs. (45) and (46) build the connection between the input and output
beams in the basis of the prolate spheroidal functions. Interestingly, this
input-output relation is very similar to the case of a two-port beam splitter.

#### Reconstructed field operators and modified point-spread function

From Eq. (45), the operator-valued coefficients $\hat{a}_{k}^{(r)}$ of the
reconstructed object [36] is found to be

$\displaystyle\hat{a}^{(r)}_{k}=\frac{\hat{f}_{k}}{(-i)^{k}\sqrt{\lambda_{k}}}=\hat{a}_{k}+\sqrt{\frac{1-\lambda_{k}}{\lambda_{k}}}\hat{b}_{k}.$
(47)

According to Eqs. (39) and (47), one can obtain the relation between the
reconstructed field operator $\hat{a}^{(r)}(z)$ and the object field operator
$\hat{a}(z)$ [36]

$\displaystyle\hat{a}^{(r)}(z)=\int^{1}_{-1}h^{(r)}(z^{\prime},z)\hat{a}(z^{\prime})dz^{\prime}+\sum^{Q-1}_{k=0}\sqrt{\frac{1-\lambda_{k}}{\lambda_{k}}}\hat{b}_{k}\phi_{k}(z),$
(48)

where $h^{(r)}(z^{\prime},z)$ denotes the reconstruction PSF and is given by

$\displaystyle
h^{(r)}(z^{\prime},z)=\sum_{k=0}^{Q-1}\phi_{k}(z^{\prime})\phi_{k}(z).$ (49)

From Eq. (49), it is easily found that the modified imaging system has a
reconstruction PSF which is related to the number $Q$. It is worthy pointing
out that as $Q$ increases, the accuracy of PSF in Eq. (49) is improved [36].
Particularly, in the limit of $Q\to\infty$, the PSF is given by
$h^{(r)}(z^{\prime},z)=\lim\limits_{Q\to\infty}\sum_{k=0}^{Q-1}\phi_{k}(z^{\prime})\phi_{k}(z)=\delta(z^{\prime}-z)$
[36], which reveals that the image is a complete replication of the object. In
other words, the PSF is totally accurate.

### Super-resolution factor

As discussed above, the model of input-output relation of the optical imaging
system has been reviewed. According this relation, we will deal with the
situation in our scheme particularly. Consider a point-like object placed at
the origin $z=0$ in the object plane in Fig. 1. Correspondingly, we assume the
focused beam in the object plane with a spatial distribution

$\displaystyle~\langle\hat{a}^{\dagger}(z)\hat{a}(z)\rangle=~\begin{cases}~~\frac{\langle\hat{n}^{\rm{w}}_{b}\rangle}{\varepsilon},&|z|\leq\varepsilon/2,\\\
~~0,&|z|>\varepsilon/2,\end{cases}$ (50)

where the width of the focused beam is very small $\varepsilon\sim 0$ and
$\langle\hat{n}^{\rm{w}}_{b}\rangle$ is actually the mean photon number of the
focused beam. It is easy to check that the total mean photon number of the
beam illuminating the object is given by
$\int_{-1}^{1}\langle\hat{a}^{\dagger}(z)\hat{a}(z)\rangle
dz=\langle\hat{n}^{\rm{w}}_{b}\rangle$ and the corresponding reconstruction
PSF is roughly equal to $h^{(r)}(0,z)$.

Assume that the object field (i.e. the focused beam) is a coherent state. In
this situation, according to Eq. (16), the SNR is exactly equivalent to the
mean photon number of the focused beam,

$\displaystyle
R=\langle\hat{n}^{\rm{w}}_{b}\rangle=\int_{-1}^{1}\langle\hat{a}^{\dagger}(z)\hat{a}(z)\rangle
dz.$ (51)

On the contrary, by combining Eqs. (15), (48), and (49), the SNR of
reconstructed object $R^{(r)}$ can be obtained,

$\displaystyle
R^{(r)}\equiv\frac{\langle\hat{n}^{(r)}\rangle^{2}}{\langle(\Delta\hat{n}^{(r)})^{2}\rangle}=\left(\sum_{k=0}^{Q-1}|a_{k}|^{2}\right)^{2}\Biggm{/}\left(\sum_{k=0}^{Q-1}|a_{k}|^{2}/\lambda_{k}\right),$
(52)

where the photon number operator of the reconstructed object
$\hat{n}^{(r)}=\int_{-1}^{1}\hat{a}^{(r)\dagger}(z)\hat{a}^{(r)}(z)dz$ and
$a_{k}=\int_{-1}^{1}\langle\hat{a}^{(r)}(z)\rangle\phi_{k}(z)dz$ represents
the coefficients of decomposition of $\hat{a}^{(r)}(z)$ over the prolate
function $\phi_{k}(z)$.

Note that the number $Q$ in Eq. (49) determines the PSF in the super-
resolving-Fourier-microscopy imaging system. It is easily checked that
increasing $Q$ improves the accuracy of the PSF in Eq. (49) which is related
to the resolution. The larger the number $Q$ is, the higher the super-
resolution of the reconstructed object achieves. Nevertheless, from Eq. (52),
it is easy to verify that with the increase of $Q$, the SNR in the
reconstructed object degrades. That is to say, the number $Q$ could not be
arbitrary large due to constraint from the decay of the SNR in the
reconstructed object. Without loss of generality, one can presume that the SNR
in the reconstructed object, no less than unity, can deliver the
reconstruction of the object.

To describe the superiority of this scheme in resolution, the comparison
between the traditional and modified schemes is performed. Correspondingly,
the super-resolution factor is introduced and is defined as the ratio of the
width of the diffraction-limited imaging PSF ($W$) in Eq. (19) to the one of
the reconstruction PSF ($W_{Q}$). In order to characterize it more
intuitively, Fig. 10 plots the diffraction-limited imaging PSF $h(z)$ and the
reconstruction PSF $h^{(r)}(0,z)$ for $Q=7$ normalized to unity at their
maxima, respectively. To quantify the degree of the super-resolution, we
introduce the half-widths $W$ and $W_{Q}$ of these two PSFs measured at their
half maxima. Then the super-resolution factor $J$ can be obtained by

$\displaystyle J=\frac{W}{W_{Q}}.$ (53)

As depicted in Fig. 10, it is easy to find that $W=1.90$, $W_{Q}=0.25$, and
$J=7.6$ which reproduces the result in Ref. [36].

Figure 10: Point-spread function as a function of $z$. The blue-dashed line
denotes the traditional optical imaging system and the red one indicates the
super-resolving-Fourier-microscopy imaging system.

### Derivation of the SNR in the presence of photon loss

The detailed derivation of SNR with loss is shown as follows. Assume that the
focus mode experiences photon loss. We mimic this loss process by passing
through a fictional two-port beam splitter with vacuum state injected in the
other input port. After the focus state passes through the beam splitter, the
output mode could be described by

$\displaystyle\hat{a}^{\rm{w}}_{b,L}$
$\displaystyle=p\hat{a}^{\rm{w}}_{b}+q\hat{a}^{\rm{V}},$ (54)

where the complex number $p$ and $q$ denote transmission and reflection
coefficients ($|q|^{2}$ indicates the photon loss rate), $\hat{a}^{\rm{V}}$
represents the annihilation operator of vacuum state.

Inserting Eq. (28) into (54), one can obtain

$\displaystyle\hat{a}^{\rm{w}\dagger}_{b,L}$
$\displaystyle=p^{\ast}\sum_{a^{\prime}=1}^{M}|t_{a^{\prime}b}^{\ast}|(\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}+\alpha^{\ast})+p^{\ast}\sum_{b^{\prime}=1}^{M}r_{b^{\prime}b}^{\ast}\hat{a}_{b^{\prime}}^{\dagger}+q^{\ast}\hat{a}^{\rm{V}\dagger},$
$\displaystyle\hat{a}^{\rm{w}}_{b,L}$
$\displaystyle=p\sum_{a^{\prime}=1}^{M}|t_{a^{\prime}b}|(\hat{a}^{\rm{SV}}_{a^{\prime}}+\alpha)+p\sum_{b^{\prime}=1}^{M}r_{b^{\prime}b}\hat{a}_{b^{\prime}}+q\hat{a}^{\rm{V}},$
(55)

The photon number operator can be cast into

$\displaystyle\hat{n}^{\rm{w}}_{b,L}=$
$\displaystyle\hat{a}^{\rm{w}\dagger}_{b,L}\hat{a}^{\rm{w}}_{b,L}$
$\displaystyle=$ $\displaystyle
p^{\ast}q\hat{n}^{\rm{w}}_{b}+q^{\ast}\hat{a}^{\rm{V}\dagger}[p\sum_{a^{\prime}=1}^{M}|t_{a^{\prime}b}|(\hat{a}^{\rm{SV}}_{a^{\prime}}+\alpha)+p\sum_{b^{\prime}=1}^{M}r_{b^{\prime}b}\hat{a}_{b^{\prime}}]+[p^{\ast}\sum_{a^{\prime}=1}^{M}|t_{a^{\prime}b}^{\ast}|(\hat{a}^{\rm{SV}\dagger}_{a^{\prime}}+\alpha^{\ast})+p^{\ast}\sum_{b^{\prime}=1}^{M}r_{b^{\prime}b}^{\ast}\hat{a}_{b^{\prime}}^{\dagger}]q\hat{a}^{\rm{V}}.$
(56)

The corresponding mean photon number is then given by

$\displaystyle\langle\hat{n}^{\rm{w}}_{b,L}\rangle=$
$\displaystyle\langle\hat{a}^{\rm{w}\dagger}_{b,L}\hat{a}^{\rm{w}}_{b,L}\rangle=|p|^{2}\langle\hat{n}^{\rm{w}}_{b}\rangle,$
(57)

The variance of photon number is found to be

$\displaystyle\langle(\Delta\hat{n}^{\rm{w}}_{b,L})^{2}\rangle=|p|^{4}\langle(\Delta\hat{n}^{\rm{w}}_{b})^{2}\rangle+|p|^{2}|q|^{2}\langle\hat{n}^{\rm{w}}_{b}\rangle.$
(58)

The average Fano factor can be then obtained

$\displaystyle\overline{F_{L}}=$
$\displaystyle\dfrac{\overline{\langle(\Delta\hat{n}^{\rm{w}}_{b,L})^{2}\rangle}}{\overline{\langle\hat{n}^{\rm{w}}_{b,L}\rangle}}=|p|^{2}\overline{F}+|q|^{2}.$
(59)

The average SNR arrives at

$\displaystyle\frac{\overline{R_{L}}}{\overline{\langle\hat{n}_{b,L}^{\rm{w}}\rangle}}=1/\overline{F_{L}}.$
(60)

## References

  * [1] Van Putten, E.G., Akbulut, D., Bertolotti, J., Vos, W.L., Lagendijk, A. & Mosk, A.P. Scattering lens resolves sub-100 nm structures with visible light. Phys. Rev. Lett. 106, 193905 (2011).
  * [2] Choi, Y., Yang, T.D., Fang-Yen, C., Kang, P., Lee, K.J., Dasari, R.R., Feld, M.S. & Choi, W. Overcoming the diffraction limit using multiple light scattering in a highly disordered medium. Phys. Rev. Lett. 107, 023902 (2011).
  * [3] Choi, Y., Yoon, C., Kim, M., Choi, W. & Choi, W. Optical imaging with the use of a scattering lens. IEEE Journal of Selected Topics in Quantum Electronics 20, 61 (2013).
  * [4] Park, C. et al. Full-field subwavelength imaging using a scattering superlens. Phys. Rev. Lett. 113, 113901 (2014).
  * [5] Yilmaz, H., Van Putten, E.G., Bertolotti, J., Lagendijk, A., Vos, W.L. & Mosk, A.P. Speckle correlation resolution enhancement of wide-field fluorescence imaging. Optica 2, 424 (2015).
  * [6] Hong, P. & Zhang, G. Heisenberg-resolution imaging through a phase-controlled screen. Opt. Express 25, 22789 (2017).
  * [7] Hong, P. Two-photon imaging assisted by a dynamic random medium. Appl. Phys. Lett. 113, 101109 (2018).
  * [8] Leonetti, M., Grimaldi, A., Ghirga, S., Ruocco, G. & Antonacci G Scattering assisted imaging. Sci. Rep. 9, 4591 (2019).
  * [9] Daniel, A., Oron, D. & Silberberg, Y. Light focusing through scattering media via linear fluorescence variance maximization, and its application for fluorescence imaging. Opt. Express 27, 21778 (2019).
  * [10] Kanaev, A.V., Watnik, A.T., Gardner, D.F., Metzler, C., Judd, K.P., Lebow, P., Novak, K.M. & Lindle, J.R. Imaging through extreme scattering in extended dynamic media. Opt. Lett. 43, 3088 (2018).
  * [11] Wiersma, D.S. Disordered photonics. Nat. Photon. 7, 188 (2013).
  * [12] Vellekoop, I.M., Lagendijk, A. & Mosk, A P Exploiting disorder for perfect focusing. Nat. Photon. 4, 320 (2010).
  * [13] Li, D. & Yao, Y. Modulating quantum fluctuations of scattered lights in disordered media via wavefront shaping. Journal of Optical Society of America B 36, 3290 (2019).
  * [14] Caves, C.M. Quantum-mechanical noise in an interferometer. Phys. Rev. D 23, 1693 (1981).
  * [15] Xiao, M., Wu, L.A. & Kimble, H.J. Detection of amplitude modulation with squeezed light for sensitivity beyond the shot-noise limit. Opt. Lett. 13, 476 (1988).
  * [16] Treps, N., Andersen, U., Buchler, B., Lam, P.K., Maitre, A., Bachor, H.A. & Fabre, C. Surpassing the standard quantum limit for optical imaging using nonclassical multimode light. Phys. Rev. Lett. 88, 203601 (2002).
  * [17] Walls, D.F. Squeezed states of light. Nature 306, 141 (1983).
  * [18] Lvovsky, A.I. Squeezed light. Photonics: Scientific Foundations, Technology and Applications 1, 121 (2015).
  * [19] Andersen, U.L., Gehring, T., Marquardt, C.& Leuchs, G. 30 years of squeezed light generation. Physica Scripta 91, 053001 (2016).
  * [20] Walls, D.F. & Milburn, G.J. Quantum Optics. (Springer Science & Business Media, 2007).
  * [21] Barnett, S.M. & Radmore, P.M. Methods in theoretical quantum optics. (Oxford University Press, 2002).
  * [22] Beenakker, C.W.J. & Patra, M. Photon shot noise. Mod. Phys. Lett. B 13, 337 (1999).
  * [23] Caves, C.M. Quantum-mechanical noise in an interferometer. Phys. Rev. D 23, 1693 (1981).
  * [24] Yurke, B, McCall, S.L. & Klauder, J.R. SU(2) and SU(1,1) interferometers. Phys. Rev. A 33, 4033 (1986).
  * [25] Xiao, M., Wu, L.A. & Kimble, H.J. Precision measurement beyond the shot-noise limit. Phys. Rev. Lett. 59, 278 (1987).
  * [26] Horrom, T., Singh, R., Dowling, J.P. & Mikhailov, E.E. Quantum-enhanced magnetometer with low-frequency squeezing. Phys. Rev. A 86, 023803 (2012).
  * [27] Otterstrom, N., Pooser, R.C. & Lawrie, B.J. Nonlinear optical magnetometry with accessible in situ optical squeezing. Opt. Lett. 39, 6533 (2014).
  * [28] Aasi, J. et al. Enhanced sensitivity of the LIGO gravitational wave detector by using squeezed states of light. Nat. Photon. 7, 613 (2013).
  * [29] Barsotti, L., Harms, J. & Schnabel, R. Squeezed vacuum states of light for gravitational wave detectors. Reports on Progress in Physics 82, 016905 (2018).
  * [30] Mehmet, M. & Vahlbruch, H. High-efficiency squeezed light generation for gravitational wave detectors. Classical and Quantum Gravity 36, 015014 (2018).
  * [31] Kolobov, M.I. & Kumar, P. Sub-shot-noise microscopy: imaging of faint phase objects with squeezed light. Opt. Lett. 18, 849 (1993).
  * [32] Kolobov, M.I. & Sokolov, I.V. Spatial behavior of squeezed states of light and quantum noise in optical images. Zh. Eksp. Teor. Fiz 96, 1945 (1989).
  * [33] Brida, G., Genovese, M. & Berchera, I.R. Experimental realization of sub-shot-noise quantum imaging. Nat. Photon. 4, 227 (2010).
  * [34] Sokolov, I.V. & Kolobov, M.I. Squeezed-light source for superresolving microscopy. Opt. Lett. 29, 703 (2004).
  * [35] Kolobov, M.I. & Fabre, C. Quantum limits on optical resolution. Phys. Rev. Lett. 85, 3789 (2000).
  * [36] Beskrovnyy, V.N. & Kolobov, M.I. Quantum limits of super-resolution in reconstruction of optical objects. Phys. Rev. A 71, 043802 (2005).
  * [37] Chen, L., Bi, S.W. & Lu, B.Z. Experimental study on the imaging of the squeezed state light at 1064 nm. Laser Physics 21, 1202 (2011).
  * [38] Vellekoop, I.M. & Mosk, A.P. Focusing coherent light through opaque strongly scattering media. Opt. Lett. 32, 2309 (2007).
  * [39] Tay, J.W., Lai, P., Suzuki, Y. & Wang, L.V. Ultrasonically encoded wavefront shaping for focusing into random media. Sci. Rep. 4, 3918 (2014).
  * [40] Ojambati, O.S., Hosmer-Quint, J.T., Gorter, K.J., Mosk, A.P. & Vos, W.L. Controlling the intensity of light in large areas at the interfaces of a scattering medium. Phys. Rev. A 94, 043834 (2016).
  * [41] Fang, L., Zhang, X., Zuo, H. & Pang, L. Focusing light through random scattering media by four-element division algorithm. Opt. Commun. 407, 301 (2018).
  * [42] Vellekoop, I.M. & Mosk, A.P. Phase control algorithms for focusing light through turbid media. Opt. Commun. 281, 2071 (2008)
  * [43] Mosk, A.P., Lagendijk, A., Lerosey, G. & Fink, M. Controlling waves in space and time for imaging and focusing in complex media. Nat. Photon. 6, 283 (2012).
  * [44] Popoff, S., Lerosey, G., Fink, M., Boccara, A.C. & Gigan, S. Image transmission through an opaque material. Nat. Commun. 1, 81 (2010).
  * [45] Li, D., Yao, Y. & Li, M. Statistical distribution of quantum correlation induced by multiple scattering in the disordered medium. Opt. Commun. 446, 106 (2019).
  * [46] Ma, C., Wang, Y.S. & An, J.H. Floquet engineering of localized propagation of light in a waveguide array. Phys. Rev. A 97, 023808 (2018).
  * [47] Xu, Y., Zhang, H., Lin, Y. & Zhu, H. Light transport behaviours in quasi-1D disordered waveguides composed of random photonic lattices. J. Mod. Opt. 64, 1215 (2017).
  * [48] Xu, Y., Zhang, H., Lin, Y. & Zhu, H. Light transmission properties in inhomogeneously-disordered random media. Annalen der Physik 529, 1600225 (2017).
  * [49] Lodahl, P. Quantum correlations induced by multiple scattering of quadrature squeezed light. Opt. Express 14, 6919 (2006).
  * [50] Beenakker, C.W.J. Random-matrix theory of quantum transport. Rev. Mod. Phys. 69, 731 (1997).
  * [51] Goodman, J.W. Statistical Optics. (Wiley, 2015).
  * [52] Van Rossum, M.C.W. & Nieuwenhuizen, T.M. Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion. Rev. Mod. Phys. 71, 313 (1999).
  * [53] Lodahl, P. Quantum correlations induced by multiple scattering of quadrature squeezed light. Opt. Express 14, 6919 (2006).
  * [54] Popoff, S.M., Lerosey, G., Fink, M., Boccara, A.C. & Gigan, S. Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media. Phys. Rev. Lett. 104, 100601 (2010).
  * [55] Yoon, J., Lee, K., Park, J. & Park, Y.K. Measuring optical transmission matrices by wavefront shaping. Opt. Express 23, 10158 (2015).
  * [56] Zhang, Z., Zhang, B., Feng, Q., He, H. & Ding, Y. Rapid measurement of transmission matrix with the sequential semi-definite programming method. Chin. Phys. B 27, 84201 (2018).
  * [57] Jang, M., Horie, Y., Shibukawa, A., Brake, J., Liu, Y., Kamali, S.M., Arbabi, A., Ruan, H., Faraon, A. & Yang, C. Wavefront shaping with disorder-engineered metasurfaces. Nat. Photon. 12, 84 (2018).
  * [58] Moore, C.P., Blaikie, R.J. & Arnold, M.D. An improved transfer-matrix model for optical superlenses. Opt. Express 17, 14260 (2009).
  * [59] Mehmet, M., Ast, S., Eberle, T., Steinlechner, S., Vahlbruch, H. & Schnabel, R. Squeezed light at $1550$ nm with a quantum noise reduction of $12.3$ dB. Opt. Express 19, 25763 (2011).
  * [60] Paschotta, R., Collett, M., K$\ddot{\rm{u}}$rz P., Fiedler, K., Bachor, H.A. & Mlynek, J. Bright squeezed light from a singly resonant frequency doubler. Phys. Rev. Lett. 72, 3807 (1994).
  * [61] Xiao, H., Rokhlin, V. & Yarvin, N. Prolate spheroidal wavefunctions, quadrature and interpolation. Inverse Problems 17, 805 (2001).
  * [62] Bertero, M. & Pike, E.R. Resolution in diffraction-limitd imaging, a singular value analysis. J. Mod. Opt. 29, 727 (1982).
  * [63] Schuurmans, F.J.P., Vanmaekelbergh, D., van de Lagemaat, J. & Lagendijk, A. Strongly photonic macroporous gallium phosphide networks. Science 284, 141 (1999).
  * [64] Kolobov, M.I. Chapter of quantum limits of optical super-resolution in Quantum Imaging. (Springer Science & Business Media, 2007).

## Acknowledgements

We would like to thank Prof. S. Sun for valuable discussions. This work was
supported by Science Challenge Program (TZ2018003-3) and National Natural
Science Foundation of China (NSFC) (61875178, 11605166).

## Author contributions statement

All authors contributed to the numerical investigations, the discussions of
the results and the preparation of the manuscript.

## Additional information

Competing interests: The authors declare no competing interests.


# Towards Multilevel Modelling of Train Passing Events on the Staffordshire
Bridge

Lawrence A. Bull  University of Cambridge, Department of Engineering,
Cambridge, CB3 0FA, UK  Chiho Jeon  University of Cambridge, Department of
Engineering, Cambridge, CB3 0FA, UK  Mark Girolami  University of Cambridge,
Department of Engineering, Cambridge, CB3 0FA, UK  The Alan Turing Institute,
The British Library, 96 Euston Road, London, NW1 2DB, UK  Andrew Duncan
Imperial College London, Department of Mathematics, London, SW7 2AZ, UK  The
Alan Turing Institute, The British Library, 96 Euston Road, London, NW1 2DB,
UK  Jennifer Schooling  University of Cambridge, Department of Engineering,
Cambridge, CB3 0FA, UK  Miguel Bravo Haro  City, University of London,
Department of Engineering, London, EC1V 0HB, UK

(May 2023)

###### Abstract

We suggest a multilevel model, to represent aggregate train-passing events
from the Staffordshire bridge monitoring system. We formulate a combined model
from simple units, representing strain envelopes (of each train passing) for
two types of commuter train. The measurements are treated as a longitudinal
dataset and represented with a (low-rank approximation) hierarchical Gaussian
process. For each unit in the combined model, we encode domain expertise as
boundary condition constraints and work towards a general representation of
the strain response. Looking forward, this should allow for the simulation of
train types that were previously unobserved in the training data. For example,
trains with more passengers or freights with a heavier payload. The strain
event simulations are valuable since they can inform further experiments
(including FEM calibration, fatigue analysis, or design) to test the bridge in
hypothesised scenarios.

Keywords: Multilevel Models, Uncertainty Quantification, Bridge Monitoring,
Simulation, Digital Twins, Meta-Analysis.

## 1 Introduction

An increasing number of bridge-monitoring projects utilise streaming telemetry
data [1, 2, 3, 4]. It is vital that we develop appropriate statistical models
to represent and extract valuable insights from these large datasets, since
the bridges constitute critical infrastructure within modern transportation
networks. The process of monitoring engineered systems via streaming data is
typically referred to as Structural Health Monitoring (SHM) and while
successful applications have been emerging in recent years, a number of
challenges remain for practical implementation [5]. During model design, these
concerns usually centre around low variance data: that is, measurements are
not available for the entire range of expected operational, environmental, and
damage conditions. Consider a bridge following construction, this will have a
relatively small dataset that should only be associated with normal operation.
On the other hand, a structure with historical data might still not experience
low-probability events – such as extreme weather or landslides.

An obvious solution considers sharing data (or information) between
structures; this has been the focus of a large body of recent work [6, 7, 8].
Since no two structures are identical, simply combining the data (complete
pooling) is rarely sufficient111Although, in certain cases, complete pooling
works, especially for more simple operations such as novelty detection with
consistent data [6, 9]. A number of statistical and ML approaches become
suitable in this setting: domain adaptation [8], transfer learning [10],
multitask learning [11], and federated learning [12] are a few, recently
explored examples.

In this article, we present a model design for data from the landmark IB5
(Staffordshire) rail bridge digital twin [2, 13]. A convenient artefact of
rail bridge monitoring is the repeated and comparable train passing events,
which readily allow information to be shared in ways that are less feasible
for foot or automobile bridges. Train events occur frequently, with an average
of 60-70 trains crossing IB5 per day. Many of these correspond to
passenger/commuter trains, the focus of this initial work. We propose a
Gaussian process representation, to work towards encoding prior knowledge of
train events within a general, flexible representation.

The layout of this paper is as follows. Section 2 outlines the existing IB5
sensing system; Section 3 summarises a decision tree to classify train events;
Section 4 formulates the multilevel Gaussian process representation; Section 5
presents initial monitoring insights from the model; Section 6 offers
concluding remarks and summarises future work.

## 2 The Norton (IB5) Bridge

It is widely accepted that monitoring and analysing telemetry data from
bridges provides valuable insights. Though in practice, there are few examples
of permanent and reliable sensing systems, which allow for the remote
acquisition of streaming data. The Norton intersection bridge 5 (IB5) is one
example, with a sensing system that continuously monitors signals, such that
train events are recorded, stored, and made accessible via an API [14].

The IB5 structure itself is a half-through, E-type bridge owned by Network
Rail, with a single skew axis span of 26.8 m, shown in Figure 1. It carries
two rail lines on the West Coast Main Line near Crewe, UK, constructed as part
of a rail redevelopment project – the Stafford Area Improvement Programme.

Figure 1: The IB5 Norton bridge, a composite structure with steel I-girders
and a cast-in-situ reinforced concrete deck, spanning 26.8m

The bridge has been instrumented with numerous sensors. Fibre Bragg Gratings
(FBGs) were installed during construction. These were attached to the main and
cross beams, and embedded in the reinforced concrete deck. More recently,
accelerometers were deployed throughout the main load-bearing steel girder, as
well as optical, humidity, and temperature sensors (and video cameras) [15].

The acquisition system monitors all sensors in real-time, capturing a
synchronised batch of raw measurements with each train passing. Over 34,000
passing events have been captured at the time of writing.

### 2.1 Strain and BWIM Data

This work considers strain data and features from a Bridge Weigh in Motion
(BWIM) system, implemented in previous phases of the IB5 project [16]. FBG raw
measurements are changes in wavelength $\Delta\lambda$, which we consider from
the west main girder of the bridge only. This is converted to the change in
micro-strain ($\Delta\mu\varepsilon$) as follows [13],

$\displaystyle\Delta\mu\varepsilon=\frac{1}{k_{\varepsilon}}\cdot\frac{\Delta\lambda}{\lambda_{0}}\times
10^{-6}$ (1)

where $k_{\epsilon}$ is the gauge factor provided by the FBG manufacturer
(0.78), and $\lambda_{0}$ scales raw measurements ($\Delta\lambda$) to present
a relative change. Since we hope to compare the bridge response for many
events, the reference $\lambda_{0}$ is an average from the start of each time
series. In practical terms, this scales each response to zero at the start of
the strain envelope. Also, note that FBG sensors are highly sensitive to
temperature variations; however, here we assume that changes in temperature
are negligible for each event, and leave out the associated compensation
(which can be found in [13]).

Figure 2 presents an example of the strain data over time. The underlying
function has a convenient interpretation: the sum of strain influence lines
per axle, which can be used to estimate each axle weight contribution [17].
Currently, strain influence lines are only utilised within the BWIM system, to
classify passing events according to train type (Section 3). However, looking
forward, parametrising the function as a sum over (influence line) bases
presents an interesting modelling perspective – discussed as future work in
Sections 4 and 6.

Figure 2: Strain time-series data (blue) and the filtered signal (black).

Conceptually, the BWIM system turns the bridge into a weighing scale,
estimating each axle weight from in-situ monitoring data. The result is a
number of insightful features, three are considered in this work: (i) number
of axles, (ii) axle weights, (iii) axle spacing, and transformations thereof.

## 3 Train Event Classification

The focus of this article is to design a general representation of train
events. With this in mind, it makes sense to initially consider typical strain
responses, and then alter the representation, if required (ideally without
increasing model complexity). An obvious starting point is 16-axle commuter
train types; including 350, 220, and 221 trains [13]. Given the BWIM features
extracted from each train passing, alongside knowledge of train design, a
simple decision tree can be used to isolate the commuter train data. The
categorisation is presented in Figure 3.

16 axles?yesaxl-sep $<$ 5.3yes350 trainsno22x trainsnoother trains

Figure 3: Train classification using BWIM data: [axl-sep] is the average axle
separation, and 22x refers to either 220 or 221 trains.

In the current work, we consider both 350 and 22x trains, highlighted with
green and  purple text at the leaves of the tree in Figure 3. An example of
ten 350 passing events is shown in Figure 4. Note the input normalisation,
using passing distance rather than time, which horizontally scales each strain
signal to be comparable, irrespective of the train-speed.

Figure 4: An example of ten passing events for 350 trains. Different colour
markers for each train.

Observing Figure 4, it becomes clear that the strain data are characteristic
and repeatable (in this case, for 350 trains). In turn, their underlying
function should be ideal for both monitoring and finding a general
representation. Here we explore a basic Gaussian process representation, which
captures the characteristic response, alongside an appropriate uncertainty
quantification for the longitudinal data. Before inference, we z-score
normalise the inputs and divide the response by its scale, as shown in Figure
5. This allows for stable inferences via Hamiltonian Monte Carlo and
convenient prior definitions in the following section.

## 4 Model Design

We now describe the multilevel Gaussian process to represent the aggregated
data. The multilevel structure will naturally extend to trains of the same
type, while also presenting options to incorporate different designs.
Importantly, with some (mean function) modifications, the model should
represent the uncertainty of event-specific functions, alongside the
uncertainty due to variations between different events of the same train type,
i.e. type-sepcific functions.

We use the notation $\mathbf{y}_{k}$ to denote the strain series (response
vector) for event $k\in\\{1,2,\ldots,K\\}$ with the corresponding inputs
$\mathbf{x}_{k}$ (passing distance). Each event is then approximated as some
function $f_{k}$ with additive noise $\boldsymbol{\epsilon}_{k}$,

$\displaystyle\mathbf{y}_{k}=f_{k}(\mathbf{x}_{k})+\boldsymbol{\epsilon}_{k}$
(2)

Following a standard Bayesian approach [18] priors are placed over the latent
functions and variables. Where the prior encodes our knowledge of the
generating process before the model is informed by data; then, by conditioning
(or updating) the priors given measurements, we find the posterior
distribution. The posterior can then be used to make predictions, which
utilise a combination of (i) domain knowledge encoded in the prior and (ii)
patterns in the observed data.

First, a shared Gaussian process (GP) prior distribution is placed over each
latent function,

$\displaystyle f_{k}\sim\textrm{GP}\big{(}m(\cdot),\;k(\cdot,\cdot)\big{)}$
(3)

This prior is fully specified by its mean $m(\cdot)$ and covariance functions
$k(\cdot,\cdot)$. The mean captures the expected path of the response, while
the covariance describes the similarity between any two outputs, with respect
to their input (i.e. the function’s smoothness). Domain expertise and
knowledge of the underlying physics can (and should) be encoded within both of
these functions; although, for now, we keep their definitions simple,

$\displaystyle m(x_{i})$ $\displaystyle=0$ (4) $\displaystyle k(x_{i},x_{j})$
$\displaystyle=\alpha^{2}\left(1+\frac{\sqrt{3}\left|\mathbf{x}_{i}-\mathbf{x}_{j}\right|}{l}\right)\exp\left(-\frac{\sqrt{3}\left|\mathbf{x}_{i}-\mathbf{x}_{j}\right|}{l}\right)$
(5)

these correspond to a zero-mean and Matérn 3/2 kernel function. These prior
functions have hyperparameters: the process variance $\alpha$ and length
scale, $l$ which have their own (hyper) priors,

$\displaystyle\alpha\sim\textrm{Truncated-Normal}_{+}(1,1),\qquad
l\sim\textrm{Half-Normal}_{+}(1)$ (6)

The prior distribution for the additive noise is also defined, and collected
with other latent variables in $\boldsymbol{\theta}$,

$\displaystyle\epsilon_{ki}\sim\textrm{Normal}(0,\sigma_{k}),\qquad$
$\displaystyle\sigma_{k}\sim\textrm{Half-Normal}_{+}(.2)$ (7)
$\displaystyle\boldsymbol{\theta}\triangleq\big{\\{}\alpha,l,\\{\sigma_{k}$
$\displaystyle\\}_{k=1}^{K}\big{\\}}$ (8)

Note, the half-normal distributions are zero-centred and the truncated normal
has a minimum value of zero, since $\\{\alpha,l,\\{\sigma_{k}\\}\\}$ must be
non-negative.

Since the hyperparameters also have a Bayesian treatment, they are also
updated given the training data. Practically, this means that the shape of the
prior adapts during inference (via the mean and covariance functions). The
above priors are set to be weakly informative, in view of the normalised
space, postulating a signal-to-noise ratio of 5 (in terms of standard
deviation).

Importantly, $k$ subscripts are missing for $\\{\alpha,l\\}$, indicating that
hyperparameters are shared between all $k$ GPs (of both 350 and 22x events). A
convenient interpretation is that the shared prior captures the
characteristics of typical commuter trains, which should be useful in
monitoring procedures (especially novelty detection). Looking forward,
$m(\cdot)$ should be specified beyond a zero-mean assumption (with parameters
contributing to $\boldsymbol{\theta}$) to present a general mean function,
which might represent a generic 350 train or 22x train specifically. The semi-
parametric formulation Kennedy and O’Hagan [19] will be considered, alongside
a basis function interpretation [17].

### 4.1 A low-rank implementation to encode boundary conditions

The number of observations per event varies quite widely, depending on the
speed of the train (750-1500). Since Gaussian processes scale badly with large
data, we implement a low-rank approximation. Following Riutort-Mayol _et al_.
[20] another basis function approximation is used for the kernel function,
which reduces the order of complexity from ${O}(N^{3})$ to ${O}(N\cdot M+M)$,
where $M$ is the number of basis functions.

As the Matérn 3/2 kernel is a stationary function, it can be approximated
using an eigendecomposition. Solin and Särkkä [21] achieve this by
interpreting $k(\cdot,\cdot)$ as the kernel of a pseudo-differential operator.
The pseudo-differential operator is then approximated via Hilbert space
methods on a compact subset $\Omega\subset\textrm{R}^{D}$, subjected to
boundary conditions. A convenient side-effect is that the boundary conditions
can be utilised to constrain the approximation given domain knowledge, see
Jones _et al_. [22] for an engineering example.

Considering the z-score normalised inputs in Figure 5, we know that most
inputs should vary between approximately $[-2,2]$ while returning to zero-
strain at the boundaries of the envelope. This scaling naturally defines an
interval in which the approximated GP is valid,

$\displaystyle\Omega\in[-L,L]=[-3,3]$ (9)

Note the buffer for $L$ (such that $L>2$) is to support stable inferences and
ensure that test data do not fall outside the approximating window of
$\Omega$.

Within $\Omega$, the kernel is then be approximated as follows [20],

$\displaystyle
k\left(x_{i},x_{j}\right)\approx\sum_{m=1}^{M}S\left(\sqrt{\lambda_{m}}\right)\phi_{m}(x_{i})\phi_{m}\left(x_{j}\right)$
(10)

where $x_{i},x_{j}\in\Omega$; $S$ is the spectral density of the kernel
function; $\lambda_{m}$ are the eigenvalues and $\phi_{m}(x)$ are the
eigenvectors of the Laplacian operator in the given domain $\Omega$. In these
experiments, we approximate with 40 basis functions ($M=40$).

For a univariate Matérn 3/2 function, $S(\cdot)$ is as follows [21],

$\displaystyle S(\boldsymbol{\omega})=\alpha\frac{4\cdot
3^{3/2}}{l^{3}}\left(\frac{3}{l^{2}}+\boldsymbol{\omega}^{\top}\boldsymbol{\omega}\right)^{-2}$
(11)

Then, following [20], the following eigenvalue problem is considered, with
Dirichlet boundary conditions respecting $\Omega$,

$\displaystyle-\nabla^{2}\phi_{m}(x)$ $\displaystyle=\lambda_{m}\phi_{m}(x),$
$\displaystyle x\in\Omega$ (12) $\displaystyle\phi_{m}(x)$ $\displaystyle=0,$
$\displaystyle x\notin\Omega$ (13)

The solution of which presents the following analytical eigenvalues and
eigenvectors,

$\displaystyle\lambda_{m}$ $\displaystyle=\left(\frac{m\pi}{2L}\right)^{2}$
(14) $\displaystyle\phi_{m}(x)$
$\displaystyle=\sqrt{\frac{1}{L}}\sin\left(\sqrt{\lambda_{m}}(x+L)\right)$
(15)

Practically, Dirichlet conditions are justified, since the strain envelopes
are bounded to zero before and after each signal (in view of the normalisation
and strain-conversion). Further details of the equivalent linear
implementation can be found in [20].

## 5 Results: Towards Monitoring Procedures

The dataset comprises 20 passing events: ten 350 trains and ten 22x trains,
sampled at random from the historical data and plotted in Figure 5. We
intentionally include an outlying event, shown by the pink markers in the left
panel. These outlying data appear to relate to an erroneous speed measurement,
and therefore, an incorrect scaling of the inputs. This outlier (within the
350 trains) is used to demonstrate simple monitoring capabilities.

Figure 5: Ten 350-type and ten 22x-type train passing events: all modelled by
the multilevel GP regression. Events are plotted with different colours. Pink
markers in the left-hand plot correspond to an outlying event.

The models were written in probabilistic programming language Stan [23] and
inferred via MCMC with the no U-turn implementation of Hamiltonian Monte Carlo
[24]. Throughout, 1000 iterations are used for burn-in and inference. The
posterior predictive distribution of each GP is shown in Figure 6. The plots
indicate a good fit to the data, however, the representation at the start/end
of each strain envelope is quite poor, since the GP must represent the entire
function without a parametrised mean – as discussed, this will be the focus of
future work.

Figure 6: Posterior predictive distributions for the multilevel GP: teal
frames are for the 350 train types and purple frames are for 22x trains.

Rather than focus on prediction metrics here, we present how the correlation
structure might inform monitoring procedures. One simple approach plots the
correlation between the (expectation) of each $f_{k}$ posterior in the
combined model, presented in Figure 7. Intuitively, there are two blocks
within this confusion matrix: one associated with each train type.
Additionally, the abnormal 350-event (index 9 within the upper left block)
would clearly be identified as outlying: data which are not flagged by the
decision tree heuristic (Figure 3).

Figure 7: Pairwise Pearson’s correlation between the (expected) posterior
distribution of each event GP in the combined model.

## 6 Concluding Remarks

We have presented initial work towards modelling multilevel data structures,
actively streaming from the (IB5) Staffordshire rail bridge. Firstly, a pre-
existing bridge-weigh-in-motion system was used to classify train events into
groups of a specific vehicle type – here, we focus on commuter trains. For
repeated crossings, the data were treated as longitudinal (or panel) data. A
hierarchical (low-rank approximation) Gaussian process was then used to model
aggregate train events, working towards a general representation of the strain
response.

Looking forward, it would be useful to parametrise the mean function, to learn
a typical representation for each train type, alongside event-specific models.
A basis function formulation using strain influence lines will be
investigated. An interpretable mean function should enable further insights
regarding the variation of strain response, as well as the option for more
expressive simulations.

## Acknowledgements

The authors would like to thank Professor C.R. Middleton, Dr. F. Huseynov and
P.R.A. Fiddler from the University of Cambridge, who led the latest phase of
the instrumentation of the Norton Bridge, thanks to the support of the Centre
for Digital Built Britain’s (CDBB). CDBB is a core partner of the Construction
Innovation Hub, funded by UK Research and Innovation (UKRI) through the
Industrial Strategy Challenge Fund (ISCF).

LAB and MG acknowledge the support of the UK Engineering and Physical Sciences
Research Council (EPSRC) through the ROSEHIPS project (Grant EP/W005816/1). AD
is supported by Wave 1 of The UKRI Strategic Priorities Fund under the EPSRC
Grant EP/T001569/1 and EPSRC Grant EP/W006022/1, particularly the Ecosystems
of Digital Twins theme within those grants & The Alan Turing Institute.

Finally, this research was possible thanks to the support of the 2022 Trimble
Fund of The University of Cambridge.

## References

  * Butler _et al_. [2016] L. J. Butler, N. Gibbons, C. Middleton, and M. Z. Elshafie.  Integrated fibre-optic sensor networks as tools for monitoring strain development in bridges during construction.  In _International Association of Bridge and Structural Engineering_ , 2016.
  * Butler _et al_. [2018] L. J. Butler, W. Lin, J. Xu, N. Gibbons, M. Z. Elshafie, and C. R. Middleton.  Monitoring, modelling, and assessment of a self-sensing railway bridge during construction.  _Journal of Bridge Engineering_ , 23(10):04018076, 2018.
  * Cross [2012] E. Cross.  _On structural health monitoring in changing environmental and operational conditions_.  PhD thesis, University of Sheffield, 2012.
  * Dervilis _et al_. [2015] N. Dervilis, K. Worden, and E. Cross.  On robust regression analysis as a means of exploring environmental and operational conditions for shm data.  _Journal of Sound and Vibration_ , 347:279–296, 2015.
  * Sohn _et al_. [2003] H. Sohn, C. R. Farrar, F. M. Hemez, D. D. Shunk, D. W. Stinemates, B. R. Nadler, and J. J. Czarnecki.  A review of structural health monitoring literature: 1996–2001.  _Los Alamos National Laboratory, USA_ , 1:16, 2003.
  * Bull _et al_. [2021] L. Bull, P. Gardner, J. Gosliga, T. Rogers, N. Dervilis, E. Cross, E. Papatheou, A. Maguire, C. Campos, and K. Worden.  Foundations of population-based SHM, part I: Homogeneous populations and forms.  _Mechanical Systems and Signal Processing_ , 148:107141, 2021.
  * Gosliga _et al_. [2021] J. Gosliga, P. Gardner, L. Bull, N. Dervilis, and K. Worden.  Foundations of population-based SHM, part II: Heterogeneous populations–graphs, networks, and communities.  _Mechanical Systems and Signal Processing_ , 148:107144, 2021.
  * Gardner _et al_. [2021] P. Gardner, L. Bull, J. Gosliga, N. Dervilis, and K. Worden.  Foundations of population-based SHM, part III: Heterogeneous populations–mapping and transfer.  _Mechanical Systems and Signal Processing_ , 149:107142, 2021.
  * Bull _et al_. [2022] L. A. Bull, P. Gardner, T. J. Rogers, N. Dervilis, E. J. Cross, E. Papatheou, A. E. Maguire, C. Campos, and K. Worden.  Bayesian modelling of multivalued power curves from an operational wind farm.  _Mechanical Systems and Signal Processing_ , 169:108530, 2022.
  * Tsialiamanis _et al_. [2020] G. P. Tsialiamanis, D. Wagg, P. A. Gardner, N. Dervilis, and K. Worden.  On partitioning of an shm problem and parallels with transfer learning.  In _Topics in Modal Analysis & Testing, Volume 8: Proceedings of the 38th IMAC, A Conference and Exposition on Structural Dynamics 2020_, page 41. Springer Nature, 2020.
  * Bull _et al_. [2023] L. Bull, D. Di Francesco, M. Dhada, O. Steinert, T. Lindgren, A. K. Parlikad, A. Duncan, and M. Girolami.  Hierarchical bayesian modelling for knowledge transfer across engineering fleets via multitask learning.  _Computer-Aided Civil and Infrastructure Engineering_ , 38(7):821–848, 2023.
  * Anaissi _et al_. [2023] A. Anaissi, B. Suleiman, and W. Alyassine.  Personalised federated learning framework for damage detection in structural health monitoring.  _Journal of Civil Structural Health Monitoring_ , 13(2-3):295–308, 2023.
  * Lin _et al_. [2019] W. Lin, L. J. Butler, M. Z. Elshafie, and C. R. Middleton.  Performance assessment of a newly constructed skewed half-through railway bridge using integrated sensing.  _Journal of Bridge Engineering_ , 24(1):04018107, 2019.
  * Broo _et al_. [2022] D. G. Broo, M. Bravo-Haro, and J. Schooling.  Design and implementation of a smart infrastructure digital twin.  _Automation in Construction_ , 136:104171, 2022.
  * Fidler _et al_. [2022] P. Fidler, F. Huseynov, M. Bravo Haro, V. Vilde, J. Schooling, and C. Middleton.  Augmenting an existing railway bridge monitoring system with additional sensors to create a bridge weigh-in-motion system and digital twin.  In _SHMII-11: 11th International Conference on Structural Health Monitoring of Intelligent Infrastructure_ , 2022.
  * Huseynov _et al_. [2023] F. Huseynov, P. Fidler, M. Bravo Haro, V. Vilde, J. Schooling, and C. Middleton.  Setting up a real-time train load monitoring system in the uk using bridge weigh-in motion technology-a case study.  In _SHMII-11: 11th International Conference on Structural Health Monitoring of Intelligent Infrastructure_ , 2023.
  * OBrien _et al_. [2006] E. J. OBrien, M. Quilligan, and R. Karoumi.  Calculating an influence line from direct measurements.  In _Proceedings of the Institution of Civil Engineers-Bridge Engineering_ , volume 159, pages 31–34. Thomas Telford Ltd, 2006.
  * Rasmussen _et al_. [2006] C. E. Rasmussen, C. K. Williams, _et al_.  _Gaussian Processes for Machine Learning_ , volume 1.  Springer, 2006.
  * Kennedy and O’Hagan [2001] M. C. Kennedy and A. O’Hagan.  Bayesian calibration of computer models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 63(3):425–464, 2001.
  * Riutort-Mayol _et al_. [2023] G. Riutort-Mayol, P.-C. Bürkner, M. R. Andersen, A. Solin, and A. Vehtari.  Practical hilbert space approximate bayesian gaussian processes for probabilistic programming.  _Statistics and Computing_ , 33(1):17, 2023.
  * Solin and Särkkä [2020] A. Solin and S. Särkkä.  Hilbert space methods for reduced-rank gaussian process regression.  _Statistics and Computing_ , 30(2):419–446, 2020\.
  * Jones _et al_. [2023] M. R. Jones, T. J. Rogers, and E. J. Cross.  Constraining gaussian processes for physics-informed acoustic emission mapping.  _Mechanical Systems and Signal Processing_ , 188:109984, 2023.
  * Carpenter _et al_. [2017] B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell.  Stan: A probabilistic programming language.  _Journal of statistical software_ , 76(1), 2017.
  * Hoffman _et al_. [2014] M. D. Hoffman, A. Gelman, _et al_.  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  _J. Mach. Learn. Res._ , 15(1):1593–1623, 2014\.