11institutetext: Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
11email<EMAIL_ADDRESS>22institutetext: DOTA, ONERA, Université
Paris Saclay, F-91123 Palaiseau, France 33institutetext: IFREMER, Laboratoire
Detection, Capteurs et Mesures (LDCM), Centre Bretagne, ZI de la Pointe du
Diable, CS 10070, 29280, Plouzane, France 44institutetext: INAF - Osservatorio
Astronomico di Arcetri

# Modeling noise propagation in Fourier-filtering wavefront sensing,
fundamental limits and quantitative comparison

V. Chambouleyron 11 O. Fauvarque 33 C. Plantet 44 J-F. Sauvage, 2211 N.
Levraud, 1122 M. Cissé, 1122 B. Neichel 11 T. Fusco, 2211

###### Abstract

Context. Adaptive optics (AO) is a technique allowing to drastically improve
ground-based telescopes angular resolution. The wavefront sensor (WFS) is one
of the key components of such systems, driving the fundamental performance
limitations.

Aims.  In this paper, we focus on a specific class of WFS: the Fourier-
filtering wavefront sensors (FFWFS). This class is known for its extremely
high sensitivity. However, a clear and comprehensive noise propagation model
for any kind of FFWFS is lacking.

Methods. Considering read-out noise and photon noise, we derive a simple and
comprehensive model allowing to understand how these noises propagates in the
phase reconstruction in the linear framework.

Results. This new noise propagation model works for any kind of FFWFS, and
allows to revisit the fundamental sensitivity limit of these sensors.
Furthermore, a new comparison between widely used FFWFS is held. We focus on
the two main used FFWFS classes: the Zernike WFS (ZWFS) and the pyramid WFS
(PWFS), bringing new understanding of their behavior.

###### Key Words.:

wavefront sensing – Fourier-filtering wavefront sensor – noise propagation

##  1 Introduction

Adaptive optics (AO) is a technique allowing to compensate for the impact of
atmospheric turbulence on telescopes that has become essential for a large
number of astrophysical applications. Motivated in particular by the hunt for
exoplanets (Kasper, 2012), AO systems pushing the limits of performance are
now being developed, called extreme AO system (XAO). These systems rely on
high-order deformable mirrors with fast real time computation. The fundamental
limits of such instruments are based on the quality of the measurements
provided by the optical device at the heart of this technique: the wavefront
sensor (WFS). One key aspect driving the XAO instruments is the WFS
sensitivity, that can limit the number of controlled modes and the speed of
the loop. The Fourier Filtering WFS (FFWFS) represents a wide class of sensors
of particular interest thanks to their superior sensitivity. From a general
point of view, a FFWFS consists of a phase mask located in an intermediate
focal plane which performs an optical Fourier filtering. This filtering
operation allows the conversion of phase information at the entrance pupil
into amplitude at the pupil plane, where a quadratic sensor is used to record
the signal (Vérinaud, 2004; Guyon, 2005; Fauvarque et al., 2016). The FFWFS
can exhibit an extra optical stage called modulation, which consists in making
the electromagnetic field spin around the tip of the filtering mask during one
integration time of the detector. This well-known operation allows to increase
dynamic range at the expense of the sensitivity (Ragazzoni, 1996).

The goal of this paper is to introduce a simple and comprehensive model for
noise propagation through any FFWFS, allowing to quantitatively compare them
within a fixed and clear framework. Inspired by an approach already employed
for the Shack-Hartmann (Rigaut & Gendron, 1992), we first derive a noise
propagation model in the linear regime that can apply to any FFWFS (section
2). We then show all the potential of this model by bringing a new
comprehension of the behavior of two widely used classes of FFWFS in section
3: the Zernike WFS (ZWFS) and the pyramid WFS (PWFS).

##  2 Noise propagation and sensitivity in the linear regime

###  2.1 Reduced intensities: a standard way to process FFWFS signal

In this paper, the FFWFS signal is processed in a standard way which leads to
a quantity named reduced intensities. If we call $I(\phi)$ the intensities
recorded on the detector for a given phase $\phi$, we build the reduced
intensities $\Delta I(\phi)$ the following way:

$\Delta I(\phi)=\frac{I(\phi)}{N_{ph}}-I_{0}$ (1)

where $N_{ph}$ is the number of photo-electrons available for the measurement
and $I_{0}$ the reference intensities. We specify that the normalization by
the total flux is absolutely required otherwise a change in flux would be
interpreted as phase aberration. In this paper, we assume working around a
flat wavefront and we therefore define the reference intensities $I_{0}$ as
the normalized intensities recorded on the WFS detector for a flat wavefront
in the pupil:

$I_{0}=\frac{I(\phi=0)}{N_{ph}}$ (2)

The reduced intensities provide the simplest way to process FFWFS signal by
taking in account all the pixels on the detector, which is a procedure
independent from the shape of the mask and structure of the modulation
(Fauvarque, 2017).

###  2.2 Noise propagation in the linear regime

To derive the noise propagation model, we make some assumptions: (i) we make
the small-phase approximation, which allow us to assume that the FFWFS are
working in their linear regime (ii) We also assume to be working with a
monochromatic light of wavelength $\lambda$. In that framework, the FFWFS
response to a modal basis $[\phi_{1},...,\phi_{i},...,\phi_{N}]$ (for which
the rms value of each mode is set to 1 in any chosen unit) is fully
characterized by the interaction matrix $\mathcal{D}$ built through a push-
pull calibration process:

$\mathcal{D}=[\delta I(\phi_{1}),...,\delta I(\phi_{i}),...,\delta
I(\phi_{N})]$ (3)

where $\delta I(\phi_{i})$ is constructed the following way:

$\delta I(\phi_{i})=\frac{\Delta I(\epsilon\phi_{i})-\Delta
I(-\epsilon\phi_{i})}{2\epsilon}$ (4)

and where $\epsilon$ is a scalar small enough to remain in the linear regime
of the sensor. By computing the reconstructor $\mathcal{D}^{{\dagger}}$
defined as the pseudo-inverse of $\mathcal{D}$ (i.e
$\mathcal{D}^{{\dagger}}=(\mathcal{D}^{t}\mathcal{D})^{-1}\mathcal{D}^{t}$)
and assuming small-phase approximation, we have the following a relationship
that allows us to reconstruct the incoming phase, providing the fact the
sensing conditions are the same as calibration (same wavelength $\lambda$,
same pupil, etc…):

$\phi=\mathcal{D}^{{\dagger}}\Delta I(\phi)\\\ $ (5)

When the signal delivered by the camera is affected by noise, one can write:

$I_{b}(\phi)=I(\phi)+b(\phi)$ (6)

where $b$ is the noise, which may depend on the phase to be measured. By
applying the reduced intensities computation on these noisy intensities
(assuming a noise-free $I_{0}$), one get:

$\Delta I_{b}(\phi)=\Delta I(\phi)+\frac{b(\phi)}{N_{ph}}=\Delta
I(\phi)+n(\phi)$ (7)

where $n(\phi)$ can be seen as a ”noise-to-signal ratio”. The smaller its
value, the less reduced intensities will be affected by noise. This quantity
propagates in the phase reconstruction the following way:

$\begin{split}\phi+\xi&=\mathcal{D}^{{\dagger}}(\Delta I(\phi)+n(\phi))\\\
\xi&=\mathcal{D}^{{\dagger}}n(\phi)\end{split}$ (8)

where $\xi$ is therefore the phase estimation error due to noise. One can
study the statistical behavior of this phase error with respect to the
statistical behavior of noise. For that, we write $<.>$ the operator averaging
on noise realizations. We can compute the error covariance matrix:

$<\xi\xi^{t}>=\mathcal{D}^{{\dagger}}<n(\phi)n^{t}(\phi)>\mathcal{D}^{{\dagger}t}$
(9)

where $<n(\phi)n^{t}(\phi)>$ is the noise covariance matrix (Rigaut & Gendron,
1992). Starting from this equation, it is possible to establish a noise
propagation model for all FFWFS sensors.

####  2.2.1 Sensitivity to read-out noise

We first focus on the behavior with respect to the read-out noise. This noise
depends only on the detector, it is therefore independent from the incoming
phase and the FFWFS characteristics (namely the mask and the modulation).
Assuming that each pixel noise is decorrelated with the others and that noise
is the same for each pixel with a standard deviation $\sigma_{ron}$ (in
$e^{-}$/px/frame), one can write:

$<n(\phi)n^{t}(\phi)>=<nn^{t}>=\frac{\sigma_{ron}^{2}}{N_{ph}^{2}}\mathbf{I}$
(10)

where $\mathbf{I}$ is the identity matrix. Equation 9 can be rewritten as:

$<\xi\xi^{t}>=\frac{\sigma_{ron}^{2}}{N_{ph}^{2}}\mathcal{D}^{{\dagger}}\mathcal{D}^{{\dagger}t}=\frac{\sigma_{ron}^{2}}{N_{ph}^{2}}(\mathcal{D}^{t}\mathcal{D})^{-1}$
(11)

By using the modal basis generated by the singular-value decomposition of
$\mathcal{D}$, $\mathcal{D}^{t}\mathcal{D}$ becomes diagonal. Let’s assume
this assertion to be still true for any modal basis. This assumption is quite
strong because it is stating (to some extend) that there is no cross talk
between measurements, which is not always the case for the considered FFWFS.
One should carefully analyze the elements of
($\mathcal{D}^{t}\mathcal{D})^{-1}$ to fully grasp the inversion process at
stake (Rigaut & Gendron, 1992). Nonetheless, assuming
$\mathcal{D}^{t}\mathcal{D}$ diagonal allows to drastically simplify the
study, because it is then possible to compute the variance of each mode due to
noise propagation with the diagonal terms of equation 11. For a given mode
$\phi_{i}$, it leads to:

$\sigma^{2}_{\phi_{i}}=\frac{\sigma_{ron}^{2}}{(\mathcal{D}^{t}\mathcal{D})_{i,i}\times
N_{ph}^{2}}$ (12)

which gives in standard deviation:

$\sigma_{\phi_{i}}=\frac{\sqrt{N_{sap}}\times\sigma_{ron}}{\sqrt{N_{sap}}\times\sqrt{(\mathcal{D}^{t}\mathcal{D})_{i,i}}\times
N_{ph}}$ (13)

In this equation, we have included $N_{sap}$ which is the number of sub-
apertures used for the measurement in the pupil. It allows to define the
sensitivity $s(\phi_{i})$ to read-out noise for a given mode $\phi_{i}$ as:

$s(\phi_{i})=\sqrt{N_{sap}}\times\sqrt{(\mathcal{D}^{t}\mathcal{D})_{i,i}}$
(14)

Adding the quantity $\sqrt{N_{sap}}$ allows to compare different masks for any
given sampling, providing that the considered mode $\phi_{i}$ is well sampled.
Equation 14 can be reformulated the following way:

$s(\phi_{i})=\sqrt{N_{sap}}\times\big{|}\big{|}\delta
I(\phi_{i})\big{|}\big{|}_{2}$ (15)

the sensitivity to RON for a given mode $\phi_{i}$ being therefore the 2-norm
of its corresponding column in the interaction matrix. One can understand from
equation 13 that a higher sensitivity will lead to less noise propagation.
Therefore, we aim at building WFS reaching the higher sensitivity possible.

####  2.2.2 Sensitivity to photon noise

Photon noise follows a Poisson distribution, therefore the variance of the
distribution is equal to the average of the intensities. Hence, we can assume
that $<n(\phi)n^{t}(\phi)>$ is diagonal and write:

$<n(\phi)n^{t}(\phi)>=\frac{\textbf{diag}(I(\phi))}{N_{ph}^{2}}$ (16)

where $\textbf{diag}(I(\phi))$ is a square matrix whose diagonal elements are
the values of $I(\phi)$. However by definition this noise depends on the
illumination pattern on the detector. It will therefore depend on the
modulation, the mask and also the incoming phase. Because we are working in
the linear regime, we can assume that:

$I(\phi)=I_{0}+\Delta I(\phi)\approx I_{0},\ \text{for }\ \phi<<1$ (17)

which means that the intensities on the detector can be approximated by the
ones corresponding to a flat wavefront. Because $I_{0}$ is already a
normalized map, $\textbf{diag}(I(\phi))=N_{ph}\times\textbf{diag}(I_{0})$. We
then can write:

$<n(\phi)n^{t}(\phi)>=<nn^{t}>=\frac{\textbf{diag}(I_{0})}{N_{ph}}$ (18)

and by injecting this formula in equation 9:

$<\xi\xi^{t}>_{\tau}=\frac{1}{N_{ph}}\mathcal{D}^{{\dagger}}\textbf{diag}(I_{0})\mathcal{D}^{{\dagger}t}$
(19)

Using the fact that
$\mathcal{D}^{{\dagger}}\mathcal{D}^{{\dagger}t}=(\mathcal{D}^{t}\mathcal{D})^{-1}$,
one can rearrange the previous equation:

$<\xi\xi^{t}>_{\tau}=\frac{1}{N_{ph}}\big{(}(\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D})^{t}\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D}\big{)}^{-1}$
(20)

As it was done for the RON noise, assuming the diagonality of the following
matrix:

$(\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D})^{t}\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D}$

allows to simply write the phase error for a given mode $\phi_{i}$:

$\sigma_{\phi_{i}}=\frac{1}{\sqrt{[(\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D})^{t}\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D}]_{i,i}}\times\sqrt{N}}$
(21)

From this last formula, we call sensitivity to photon noise $s_{\gamma}$ the
following quantity:

$s_{\gamma}(\phi_{i})=\sqrt{[(\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D})^{t}\textbf{diag}(1/\sqrt{I_{0}})\mathcal{D}]_{i,i}}$
(22)

In a more comprehensive way, the photon noise sensitivity to the mode
$\phi_{i}$ is calculated by the following operation:

$s_{\gamma}(\phi_{i})=\big{|}\big{|}\frac{\delta
I(\phi_{i})}{\sqrt{I_{0}}}\big{|}\big{|}_{2}$ (23)

One could wonder if there should be a special procedure for non-illuminated
pixels in $I_{0}$. In that case, these pixels exhibits a null value in the
interaction matrix and add no contribution to the sensitivity computation. The
formulas given 13 and 21 give a relationship between the phase error, the
number of photons available for the measurement, the RON, the number of sub-
apertures and the intrinsic properties of the studied FFWFS. The noise
propagation for a given mode $\phi_{i}$ is then quantified through the formula
:

$\boxed{\sigma^{2}_{\phi_{i}}=\frac{N_{sap}\times\sigma_{ron}^{2}}{s^{2}(\phi_{i})\times
N_{ph}^{2}}+\frac{1}{s^{2}_{\gamma}(\phi_{i})\times N_{ph}}}$ (24)

Where $\sigma^{2}_{\phi_{i}}$ is given in $u^{2}$, where $u$ is the unit used
to normalize the modes to rms 1 when building the interaction matrix (usually
radians or nanometers).We mention that $s(\phi_{i})$ and
$s_{\gamma}(\phi_{i})$ have no reason to be equal. When comparing different
FFWFS in term of sensitivity, we thus need to compare their RON sensitivity
and their photon noise sensitivity. This kind of formula was already derived
in (Guyon, 2005) for simplified and idealized cases of 4PWFS and ZWFS.
Providing their interaction matrix and the linear regime, the formula 24
provides a unified scheme that works for any kind of modulation and mask, and
doesn’t use an analytical model: the computed sensitivities are therefore more
accurate with this derivation.

###  2.3 Maximum sensitivity limit

The columns of the interaction matrix are encoding the derivative of the
signal with respect to the phase. Hence, the sensitivity expressions given
equations 15 and 23 are actually describing a metric called Fisher information
(Plantet, 2015). This metric is bounded by the Cramer-Rao bound, and it is
then possible to show that both sensitivities defined here (for RON and photon
noise) have a maximum value of 2 (Paterson, 2008; Bouchet et al., 2021).

$\begin{split}0\leq\ &s\leq 2\\\ 0\leq\ &s_{\gamma}\leq 2\end{split}$ (25)

These boundaries are therefore useful to provide an absolute reference to
compare FFWFS sensitivities. We add the fact that the Cramer-Rao bound
actually works for any kind of WFS, and these boundaries extend out of the
FFWFS scope.

##  3 Sensitivity comparison of different FFWFS

The goal of this section is to use the metrics defined previously to compare
well-known and widely used FFWFS. We will focus on two main class of FFWS: the
Zernike WFS (ZWFS) and the Pyramid WFS (PWFS). We are going to show that the
use of the RON sensitivity in one hand and the photon noise sensitivity in the
other hand will highlight some important properties of these sensors. For our
comparison, we compute sensitivities with respects to the Fourier modes which
are simply defined by the sum of a cosine and a sine carrying a given spatial
frequency f (Fauvarque et al., 2019). The following quantity then encodes the
sensitivity:

$s_{\textbf{f}}=\frac{1}{\sqrt{2}}\sqrt{s\Big{(}\cos_{\textbf{f}}\Big{)}^{2}+s\Big{(}\sin_{\textbf{f}}\Big{)}^{2}}\,.$
(26)

###  3.1 Zernike class

The ZWFS is known for its extremely high sensitivity while having a low
dynamic range, which makes it a unique sensor for second stage adaptive optics
(AO) systems or quasi-static aberrations calibration sensor (Jensen-Clem et
al., 2012; N’Diaye et al., 2016; Ruane et al., 2020). This sensor is composed
of a focal plane mask made of a phase shifting dot fully described by two
parameters: its diameter $p$ and phase shift $d$. In a previous paper
(Chambouleyron et al., 2021), we extensively described and analyzed the impact
the influence of the mask parameters on sensitivity, but only for RON. We
extend this study by applying our metric developed for photon noise
sensitivity to have better understanding of the behavior of this sensor.
Usually, the ZWFS is used in what we will call here the classical ZWFS
configuration: $p=1.06\ \lambda/D$ and $d=\pi/2$ (Vigan, A. et al., 2019). We
remind the reader that the linear signal for the ZWFS is fully located into
the footprint of the pupil image onto the detector plane, meaning that all
photons outside the footprint are left unused in the linear regime (Fauvarque,
2017).

####  3.1.1 Diameter

It was shown that increasing the ZWFS diameter can significantly increase the
RON sensitivity for high spatial frequencies, at the expense of low spatial
frequency located inside the dot perimeter (Chambouleyron et al., 2021). Our
goal is to assess if this statement also applies for the photon noise
sensitivity. In this study, we fix the dot phase shift to $\pi/2$. Figure 1
shows the sensitivity for RON and photon noise computed in the case of three
different diameters: $p=1.06\ \lambda/D$, $p=2\ \lambda/D$ and $p=5\
\lambda/D$. We see that the behavior for RON described in Chambouleyron et al.
(2021) is also seen for photon noise: high sensitivity gain for high spatial
frequencies and drop in sensitivity for frequencies located inside the dot. We
remind that changing mask parameters also affects the reference intensities
$I_{0}$, which are plotted on figure 1 for the considered configurations.
Finally, all conclusions that were brought in the Z2WFS study are still valid
for the photon noise sensitivity.

Figure 1: Senstivities to RON and photon noise for different configuration of
the ZWFS. Phase shift fixed at $d=\pi/2$ and diameter $p=1.06\ \lambda/D$,
$p=2\ \lambda/D$ and $p=5\ \lambda/D$. We can see that the two kind of
sensitivities follow the same trends.

####  3.1.2 Phase-shift

We now focus on the impact of the dot phase shift on sensitivity. It is
already know that for RON sensitivity, the maximum is reached for $d=\pi/2$ no
matter what we choose as diameter size $p$. Is this assertion also true for
photon noise sensitivity ? As a first step to answer this question, we fixed
$p=1.06\lambda/D$ and let the depth vary between 0 and $2\pi$. We focus on the
sensitivity for a given spatial frequency corresponding to 6 cycles in the
diameter. The choice of this spatial frequency is arbitrary, but motivated by
the fact that its footprint in the focal plane is located far from the phase
shifting dot boundaries and not inside the phase-shifting area. The ZWFS
behavior with respect to this frequency is then representative of all the
high-order modes. The results are given Figure 2. We can see that for the RON
sensitivity, the maximum value is indeed reached at $d=\pi/2$. We also notice
the coronagraphic configuration corresponding to $d=\pi$ for which the
intensities on the detector don’t have any linear dependence with the entrance
phase (Roddier, 1997). One can also notice the symmetric behavior of this
curve, which is perfectly logical because of the $2\pi$ wrapping of the phase.
For the photon noise sensitivity, the curve exhibits a slightly different
behavior: the maximum value is not reached for $\pi/2$, but for $d=0.8\pi$.
This configuration is very instructive, because it gives us a deep insight on
how optimal RON sensitivity configurations can differ from optimal photon
noise sensitivity configurations.

Figure 2: Sensitivities for a given spatial frequency corresponding to 6
cycles in the diameter for $p=1.06\ \lambda/D$ while varying phase shift $d$.
The maximum values are reached for different value of $d$.

To understand why $d=0.8\pi$ brings better performance with respect to photon
noise, we have to focus on the structure of the reference intensities. In the
case of a WFS sensitivity for photon noise, photons are the source of signal
but also of noise. We showed the importance of the reference intensities in
the previous section (Equation 22). In the configuration $p=1.06\lambda/D$ and
$d=0.8\pi$, the sensor exhibits a behavior close to the coronagraphic one.
Therefore, the reference intensities have a limited number of photons in the
pupil footprint, as it is shown in Figure 3 which shows the percentage of
photons available in the pupil footprint of the reference intensities while
increasing phase shift for diameter $p=1.06\lambda/D$. A way to understand the
photon noise sensitivity curve’s shape presented in figure 2: assuming uniform
reference intensities for the configuration $p=1.06\lambda/D$ while changing
the depth $d$, it can be interpreted as the ratio between the RON sensitivity
curve given on the same figure and the pupil footprint flux curve shown figure
3. The maximum found at $d=0.8\pi$ can therefore be seen as a trade-off
between high signal ($d=\pi/2$) and low noise produced by photons ($d=\pi$).
This configuration is a nice example of ”dark-WFS” concept, a class of sensor
for which reference intensities would gather a limited number of photons
(Ragazzoni, 2015).

Figure 3: Percentage of photons kept in the pupil footprint for a
configuration with diameter $p=1.06\lambda/D$ while changing phase shift $d$.
The reference intensities for $d=0.8\pi$ are represented on the right.

We then plot the sensitivity for all the spatial frequencies for different
configurations, as it was done in the previous subsection when varying the
diameter $p$. We continue exploring the instructive $d=0.8\pi$ configuration
as an example by comparing here the classical ZWFS configurations with two
other set of parameters for the dot shape: $p=1.06\lambda/D$ with $d=0.8\pi$,
and $p=2\lambda/D$ with $d=0.8\pi$. Results are given in Figure 4. For the
$p=1.06\lambda/D$ with $d=0.8\pi$ configuration, the RON sensitivity is lower
for all spatial frequencies compared to classical ZWFS while it is always
higher for photon noise sensitivity. However, when increasing dot diameter to
a value of $p=2\lambda/D$ while having $d=0.8\pi$ decrease the gain for photon
noise sensitivity for all frequencies. This can be easily explained by the
fact that in this configuration, the large size of the dot prevent from a
coronagraphic behavior and does not remove a significant amount of photons
from the pupil footprint in reference intensities.

Figure 4: Sensitivities to RON and photon noise for different configuration of
the ZWFS. The ”dark” ZWFS configuration correspond to $p=1.06\lambda/D$ and
$d=0.8\pi$.

One could push further the analysis by studying the photon noise sensitivity
behavior for a full set of different diameter $p$ and phase shift $d$ and for
all spatial frequencies. This study is out of the scope of this paper. We
demonstrated the benefit to use two different sensitivity metrics for RON and
photon noise and therefore introduced the required tools to lead a fine study
of the ZWFS class.

###  3.2 Pyramid class

In this section, we study another FFWFS class: the PWFS. This class holds on a
pyramidal shape of the phase mask, which splits the electromagnetic field in
different beams (Ragazzoni, 1996). For the PWFS, the extra optical module
called ”modulation” presented in the introduction is added. We study the
impact of the following parameters on sensitivities: (i) the number of faces:
it determines the number of pupil images formed on the detector (ii) the angle
of the faces: it determines the relative separation between pupil images.

####  3.2.1 Number of faces

We start this study by analyzing the impact of the number of faces on PWFS
sensitivities. To that end, we use five different configurations: a three-
sided pyramid (3PWFS), the most widely used: a four-sided pyramid (4PWFS), a
five-sided pyramid (5PWFS), a six-sided pyramid (6PWFS) and finally a pyramid
with a infinite number of faces, forming what we call an axicon. As an
example, we show Figure 5 the phase masks and the corresponding reference
intensities of the 4PWFS, 5PWFS and axicon in the non-modulated case.

Figure 5: PWFS with different number of faces. Top: arg(mask) for 4PWFS, 5PWFS
and axicon. Botton: Reference intensities associated for the non-modulated
case.

We assess the impact of the number of faces on sensitivities by considering
two modulation cases: a modulated case with a given modulation radius
$r_{mod}=3\ \lambda/D$ and the non-modulated case $r_{mod}=0\ \lambda/D$.
Because the considered masks are not axi-symmetric (except for axicon), we
could analyze their sensitivities in the frequencies space by extracting 2D
sensitivity maps. For an easier comparison, we will extract radial cuts from
these maps. For that, we choose a direction which does not lie on any edges
for the considered configurations.

RON sensitivity: We analyze first the RON sensitivity. Figure 6 shows the
results of our study for different PWFS cases. We point out that the
sensitivity scale stops at 1 for better readability, but we remind that the
maximum sensitivity is 2.

Figure 6: RON sensitivities for the PWFS when changing the number of faces.
Left: non-modulated cases. Right: modulation at $3\ \lambda/D$.

In the non-modulated case, we note that reducing the number of faces increase
the sensitivity. The gain in sensitivity is logically proportional to the
number of useful pixels: if the linear signal is reaching more pixels, the
PWFS will let propagate more RON. Therefore, the plateau for the 3PWFS is 4/3
higher than the 4PWFS for example. When modulating (right panel Figure 6), we
observe the same relative behavior: the less faces, the better. We can see the
well-known loss of sensitivity for spatial frequencies located in the
modulation radius. We also spot a global loss for all spatial frequencies. One
can also note the different behavior for the axicon: its sensitivity is
decreasing drastically compared to other configurations.

Photon noise sensitivity: We now analyze the impact of the number of faces
with respect to photon noise sensitivity. Results are given Figure 7. Let’s
start our analysis with the modulated case (right panel). As for the RON
sensitivity, we see that the sensitivity increases when reducing the number of
faces. When modulating, the distribution of photons is mainly uniform in the
pupil images. The ratio between sensitivities is then proportional to the
square root of the number of faces ratio, because signal-to-noise ratio
increase with the square root of the number of photons. For instance, the
3PWFS exhibits a increased sensitivity of $\sqrt{4}/\sqrt{3}$ with respect to
the 4PWFS. However, when looking at the non-modulated cases (Figure 7, left
panel), all the different PWFS configuration exhibit the same sensitivities.
This behavior can be explained by the fact that in the non-modulated case,
diffraction effects by the edges of the masks have a significant impact on the
intensities recorded on the detector.

This study shows that the 3PWFS exhibits a better sensitivity than the other
configurations for the RON sensitivity (as expected and already tackled in
previous studies (Schatz et al., 2021)). We show here for the first time that
this superiority is also true for photon noise sensitivity, but only when
modulated.

We also emphasize a key aspect when comparing PWFS to ZWFS: the non-modulated
PWFS have a better photon noise sensitivity than the classical ZWFS ($p=1.06\
\lambda/D$ and $d=\pi/2$). This observation seems important, as the classical
ZWFS is often presented as the most sensitive FFWFS. We see here that in its
classical configuration, that is far from being the case because its
sensitivities are not reaching the maximum value of 2.

Figure 7: Photon noise sensitivities for the PWFS when changing the number of
faces. Left: non-modulated cases. Right: modulation at $3\ \lambda/D$.

####  3.2.2 Angle

We now study the impact of the pyramid angle on sensitivities. For this study,
we choose a non-modulated 4PWFS. We already know that decreasing PWFS angle
can change its behavior as soon as the pupils are overlapping, creating a
configuration called the flattened-PWFS (FPWFS) which can be interpreted as a
mix between a PWFS and a shearing-interferometer. We will call a classical
PWFS the one showing no overlapping. It was already demonstrated (Fauvarque et
al., 2015) that the FPWFS class exhibit oscillating sensitivity curve, whose
period depends on the separation between pupil images (illustration given
figure 8).

Figure 8: The Flattened PWFS: periodicity of the oscillating sensitivity
depends on the pupil overlapping ratio.

RON sensitivity: Following the same procedure as before, we start this study
with the sensitivity to RON. This aspect of the FPWFS has been already
studied. We take this example of 5 configurations with different overlapping
ratio (percentage of pupil diameter), and plot their respective sensitivity
curves figure 9. We find back the well-know oscillating curves. One can notice
that there are some configurations which show a better sensitivity to RON with
respect to classical PWFS for all spatial frequencies. However, one could
wonder if such configurations are still more sensitive than the classical PWFS
when considering the photon noise sensitivity. Thanks to our new propagation
model, we can assess this question for the first time.

Figure 9: RON sensitivities for different FPWFS configurations (non-modulated
cases).

Photon noise sensitivity: We compute sensitivity to photon noise with the same
overlapping configurations studied just before. Sensitivity curves are given
figure 10. These sensitivity curves show a very different story: this time, we
can see that there is no configuration always above the classical PWFS. If
some configuration have a better sensitivity for some frequencies, it is
always at the expense of others. The curves are oscillating around the
classical PWFS one. These results bring a new understanding of the FPWFS
class, showing that the FPWFS configurations are not as interesting as one
could have expected for the photon noise sensitivity. The sensitivity curves
given figures 9 and 10 can be compared to the ones computed for the ZWFS class
given figure 2, where one can notice that in its optimal configurations, the
ZWFS class is out-reaching the PWFS for a large set of spatial frequencies.

Figure 10: Photon noise sensitivities for different FPWFS configurations (non-
modulated cases).

##  4 Conclusion and perspectives

In this paper, we have derived a simple and comprehensive noise propagation
model for all FFWFS. This model relies on the small phases approximation and
allows to define two kinds of sensitivities which describe FFWFS behavior with
respect to noise. There is the RON sensitivity which is computed through the
interaction matrix, and the photon noise sensitivity which is computed through
the interaction matrix and the reference intensities. These two bounded
quantities can be used to link the number of photons available for the phase
measurement and the estimation error variance. It provides a unified model
which allows to quantitatively compare all FFWFS and that goes beyond results
given in Guyon (2005). We demonstrated the usefulness of this model by
applying it to ZWFS and PWFS, highlighting for the first time some interesting
properties of these sensors and exhibiting unprecedented comparisons.

The fact that we now have well-defined criteria for the RON and photon noise
sensitivities can lead us to try optimizing mask shape to reach sensitivities
close to the maximum value of 2 (Chambouleyron et al., 2022).

These noise propagation model was derived assuming linear regime. We did not
intend here to analyze dynamic range of the studied FFWFS. Therefore this
study is just giving a part of the full description these sensors. On top of
that, it is well known that FFWFS are not working in their linear regime when
working in closed-loop, leading to non-linear effects than can be described at
first order by a loss of sensitivity for each mode (the so-called optical gain
effect Chambouleyron et al. (2020)). One could therefore define the effective
sensitivities, which would be the sensitivities defined in this paper,
weighted by the loss of sensitivities induced by non-linearities.

##  5 acknowledgements

This work benefited from the support of the WOLF project ANR-18-CE31-0018 of
the French National Research Agency (ANR). It has also been prepared as part
of the activities of OPTICON H2020 (2017-2020) Work Package 1 (Calibration and
test tools for AO assisted E-ELT instruments). OPTICON is supported by the
Horizon 2020 Framework Programme of the European Commission’s (Grant number
730890). Authors are acknowledging the support by the Action Spécifique Haute
Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES. Vincent
Chambouleyron PhD is co-funded by ”Région Sud” and ONERA, in collaboration
with First Light Imaging. Finally, part of this work is supported by the LabEx
FOCUS ANR-11-LABX-0013.

## References

  * Bouchet et al. (2021) Bouchet, D., Dong, J., Maestre, D., & Juffmann, T. 2021, Phys. Rev. Applied, 15, 024047
  * Chambouleyron et al. (2020) Chambouleyron, V., Fauvarque, O., Janin-Potiron, P., et al. 2020, 644, A6
  * Chambouleyron et al. (2022) Chambouleyron, V., Fauvarque, O., Plantet, C., et al. 2022, in Adaptive Optics Systems VIII, ed. L. Schreiber, D. Schmidt, & E. Vernet, Vol. 12185, International Society for Optics and Photonics (SPIE), 121852T
  * Chambouleyron et al. (2021) Chambouleyron, V., Fauvarque, O., Sauvage, J. F., et al. 2021, A&A, 650, L8
  * Fauvarque (2017) Fauvarque, O. 2017, PhD thesis, thèse de doctorat dirigée par Fusco, Thierry et Neichel, Benoît Physique et sciences de la matière. Astrophysique et cosmologie Aix-Marseille 2017
  * Fauvarque et al. (2019) Fauvarque, O., Janin-Potiron, P., Correia, C., et al. 2019, J. Opt. Soc. Am. A, 36, 1241
  * Fauvarque et al. (2015) Fauvarque, O., Neichel, B., Fusco, T., & Sauvage, J.-F. 2015, Opt. Lett., 40, 3528
  * Fauvarque et al. (2016) Fauvarque, O., Neichel, B., Fusco, T., Sauvage, J.-F., & Girault, O. 2016, Optica, 3, 1440
  * Guyon (2005) Guyon, O. 2005, The Astrophysical Journal, 629, 592–614
  * Jensen-Clem et al. (2012) Jensen-Clem, R., Wallace, J. K., & Serabyn, E. 2012, in 2012 IEEE Aerospace Conference, 1–7
  * Kasper (2012) Kasper, M. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8447, Adaptive Optics Systems III, ed. B. L. Ellerbroek, E. Marchetti, & J.-P. Véran, 84470B
  * N’Diaye et al. (2016) N’Diaye, M., Vigan, A., Dohlen, K., et al. 2016, in SPIE Astronomical Telescopes + Instrumentation, Vol. 9909, Edinburgh, United Kingdom, 99096S
  * Paterson (2008) Paterson, C. 2008, Journal of Physics: Conference Series, 139, 012021
  * Plantet (2015) Plantet, C. 2015, Optics Express, 23, 28619
  * Ragazzoni (1996) Ragazzoni, R. 1996, Journal of Modern Optics, 43, 289
  * Ragazzoni (2015) Ragazzoni, R. 2015, AO4ELT-4
  * Rigaut & Gendron (1992) Rigaut, F. & Gendron, E. 1992, A&A, 261, 677
  * Roddier (1997) Roddier. 1997, Publications of the Astronomical Society of the Pacific, 109, 815
  * Ruane et al. (2020) Ruane, G., Wallace, J. K., Steeves, J., et al. 2020, Journal of Astronomical Telescopes, Instruments, and Systems, 6
  * Schatz et al. (2021) Schatz, L., Males, J. R., Correia, C., et al. 2021, arXiv e-prints, arXiv:2109.06386
  * Vérinaud (2004) Vérinaud, C. 2004, Optics Communications, 233, 27
  * Vigan, A. et al. (2019) Vigan, A., N´Diaye, M., Dohlen, K., et al. 2019, A&A, 629, A11


We describe a new phenomenon in the study of the real-time path integral,
where complex classical paths hit singularities of the potential and need to
be analytically continued beyond the space for which they solve the boundary
value problem. We show that the behavior is universal and central to the
problem of quantum tunneling. These analytically continued complex classical
paths enrich the study of real-time Feynman path integrals.

# Crossing singularities in the saddle point approximation

Job Feldbrugge<EMAIL_ADDRESS>Higgs Centre for Theoretical Physics,
James Clerk Maxwell Building, Edinburgh EH9 3FD, UK  Dylan L. Jow
<EMAIL_ADDRESS>Canadian Institute for Theoretical Astrophysics,
University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
Department of Physics, University of Toronto, 60 St. George Street, Toronto,
ON M5S 1A7, Canada  Dunlap Institute for Astronomy & Astrophysics, University
of Toronto, AB 120-50 St. George Street, Toronto, ON M5S 3H4, Canada  Ue-Li
Pen ${}^{2,\,3,\,4,\,}$<EMAIL_ADDRESS>Institute of Astronomy and
Astrophysics, Academia Sinica, Astronomy-Mathematics Building, No. 1, Section
4, Roosevelt Road, Taipei 10617, Taiwan  Perimeter Institute for Theoretical
Physics, 31 Caroline St. North, Waterloo, ON, Canada N2L 2Y5  Canadian
Institute for Advanced Research, CIFAR program in Gravitation and Cosmology

Introduction: In the sum over histories formulation of quantum physics,
evolution arises through the constructive interference of histories around
classical paths. Classically allowed transitions are governed by real
classical paths. Classically forbidden transitions – such as quantum tunneling
– are often associated with complex solutions of the equations of motion
interpolating between boundary conditions. Indeed, complex classical paths –
sometimes known as instantons – are frequently studied in quantum mechanics,
quantum field theory, and theories of quantum gravity. However, it is often an
open question whether a complex classical solution is relevant to the real-
time Feynman path integral and contributes to the corresponding amplitude. As
we show in this letter, the situation is even more intricate, as the main
contribution to the path integral may come from analytically continued complex
classical paths that do not solve the classical boundary value problem. In
this letter, we demonstrate this phenomenon for several quantum systems and
discuss its various implications.

We consider tunneling of a wave packet through a barrier centered at $x=0$.
This can be considered an initial value problem of the time-dependent
Schrödinger equation, with an initial Gaussian, minimum-uncertainty wave
packet on the left with positive momentum. For high-momentum wave packets,
most of the probability passes over the barrier, similar to the classical
equation of motion. For small momenta, with energies less than the height of
the barrier, most of the wave packet reflects back, and a small amount
tunnels.

The time-dependent Schrödinger equation,

$\displaystyle i\hbar\frac{\partial\psi_{t}(x)}{\partial
t}=\hat{H}\psi_{t}(x),$ (1)

is solved by the propagator

$\displaystyle\psi_{T}(x_{1})$
$\displaystyle=\int\psi_{0}(x_{0})G[x_{0},x_{1};T]\mathrm{d}x_{0}\,,$ (2)
$\displaystyle G[x_{0},x_{1};T]$
$\displaystyle=\int_{x(0)=x_{0}}^{x(T)=x_{1}}e^{iS[x]/\hbar}\mathcal{D}x\,.$
(3)

The second equality (3) is the Feynman path integral formulation, which for
small $\hbar$ is solved by stationary points of the action $\delta S/\delta
x=0$ resulting in Euler-Lagrange (EL) equations of motion

$\displaystyle m\ddot{x}=-V^{\prime}(x)\,,$ (4)

with the mass $m$, the Heaviside theta function $\Theta_{H}$ and the dots
denoting time derivatives, satisfying boundary conditions $x(0)=x_{0},\
x(T)=x_{1}$ Feynman:1948 ; Feynman:1965 ; Feynman:1948b . In the Picard-
Lefshetz (PL) picture, one deforms the real space path integral (3) contour
into the sum of complex path segments (thimbles), each of which is non-
oscillatory Turok:2014 ; Feldbrugge:2023 . For small $\hbar$, each thimble is
dominated by its saddle point. This saddle may be complex, even though the
original integral was defined over real paths. It is tempting to consider
complex solutions (saddles) of (4) satisfying the real boundary conditions for
complex energies and evaluate the propagator (3) using the corresponding
action Schulman:2012 . Two technical complications arise: 1. while PL saddles
are solutions to the EL equations, the converse is generally not true, and 2.,
no constructive example has been found where the propagator from the
Schrödinger equation agrees with the complex saddle.

In this paper we revisit this problem, systematically examining the real and
complex paths for fixed $x_{0}$ while varying $x_{1}$. The Rosen-Morse barrier
(see below) has a particularly simple exact solution (saddle) set, and
corresponding exact actions at those saddles. This allows us to examine the
propagator under these continuous deformations.

Several conceptual issues arise: 1. complex saddles depend on variable
rescalings, and 2. the evaluation of the action at a complex saddle gives a
different value from the analytic continuation of the action along the
deformation path of (real) $x_{1}$. We find these are reconciled due to branch
cuts and singularity crossings for the saddles.

This Letter addresses this novel approach of treating branch cuts in the
action, explicated by the exact Rosen-Morse barrier solution, but
generalizable to generic potentials. The analysis is simplified using the
propagator in Energy $E$, conjugate to time $T$.

The Rosen-Morse barrier: The symmetric Rosen-Morse Rosen:1932  or modified
Pöschl-Teller system Poschl:1933 , describing the evolution of a non-
relativistic particle in the potential

$\displaystyle V(x)=V_{0}\,\text{sech}^{2}(x)\,,$ (5)

is one of the few fully solvable non-linear quantum systems. For simplicity,
we will restrict ourselves to the barrier problem with
$V_{0}\geq\hbar^{2}/(8m)$. The real-time energy propagator of a non-
relativistic particle with mass $m$ propagating form $x_{0}$ to $x_{1}$ at
energy $E$ is solved in closed form Kleinert:1992 ,

$\displaystyle K[x_{1},x_{0};E]$
$\displaystyle=\int_{0}^{\infty}\int_{x(0)=x_{0}}^{x(1)=x_{1}}e^{i(S[x]+ET)/\hbar}\mathcal{D}x\,\mathrm{d}T$
(6) $\displaystyle=-\frac{im\Gamma(-ik_{E}-N)\Gamma(-ik_{E}+N+1)}{\hbar}$
$\displaystyle\phantom{=}\times P_{N}^{ik_{E}}(\tanh
x_{>})P_{N}^{ik_{E}}(-\tanh x_{<})\,,$ (7)

where the path integral ranges over the continuous paths interpolating between
$x(0)=x_{0}$ and $x(1)=x_{1}$ parametrized by $\lambda\in[0,1]$, and the
action

$\displaystyle
S[x]=\int_{0}^{1}\left[\frac{mx^{\prime}(\lambda)^{2}}{2T}-TV(x(\lambda))\right]\mathrm{d}\lambda\,,$
(8)

with the propagation time $T$, and the prime denoting derivatives with respect
to the parameter $\lambda$. The closed-form solution includes the associated
Legendre function $P_{\lambda}^{\mu}(x)$ with the parameter
$N=-\tfrac{1}{2}+\tfrac{i}{2\hbar}\sqrt{8mV_{0}-\hbar^{2}}\,,$ the minimum
$x_{<}$ and the maximum $x_{>}$ of the boundary conditions $(x_{0},x_{1})$,
and dimensionless momentum $k_{E}=\sqrt{2mE}/\hbar$. The Legendre function
with degree $\text{Re}[\lambda]=-1/2$ is known as the conical or Mehler
function.

The associated classical system is given by the Euler-Lagrange and Hamilton-
Jacobi equations

$\displaystyle\frac{\delta S}{\delta x}=0\,,\quad E+\frac{\partial S}{\partial
T}=0\,,$ (9)

which for the Rosen-Morse system assumes the form

$\displaystyle\frac{mx^{\prime\prime}}{T^{2}}=\frac{2V_{0}\tanh
x}{\cosh^{2}x}\,,\quad
E=\frac{mx^{\prime}(0)^{2}}{2T^{2}}+\frac{V_{0}}{\cosh^{2}x_{0}}\,.$ (10)

The Euler-Lagrange equation is solved by

$\displaystyle\sinh
x_{C}(\lambda)=c\sqrt{\frac{E-V_{0}}{E}}\sinh\left[\sqrt{\frac{2E}{m}}(T\lambda-C)\right]\,,$
(11)

with the sign $c=\pm 1$ and the shift parameter $C$, with the associated
classical action

$\displaystyle S_{C}=ET-\sqrt{2mV_{0}}$
$\displaystyle\times\tanh^{-1}\left[\sqrt{\frac{V_{0}}{E}}\tanh\left[\sqrt{\frac{2E}{m}}(T\lambda-C)\right]\right]\Bigg{|}_{0}^{1}\,.$
(12)

We solve the initial value problem, where the Hamilton-Jacobi equation yields
the initial velocity

$\displaystyle
v_{0}=\pm\sqrt{\frac{2}{m}\left[E-V_{0}\,\text{sech}^{2}x(0)\right]}\,T\,,$
(13)

for complex $T$ to find the solutions to the boundary value problem associated
with the energy propagator. Candidate classical solutions lay on the curves
for which the final position $x(\lambda=1)$ is real (see the green curves in
the left panel of fig. 2).

(a) $x_{1}$

(b) $\sinh x_{1}$

Figure 1: The final position of the classical paths starting at $x_{0}=-5$ for
an energy $E=0.9$, mass $m=1$, and a barrier strength $V_{0}=1$ in the complex
$T$ plane. The green lines indicate times that result into a real final
position.

(c) Real classical paths

(d) Complex classical paths

Figure 2: The classical paths of the energy propagator at $E=0.9$ for a
particle starting at $x_{0}=-5$ with mass $m=1$ in a barrier with strength
$V_{0}=1$. Left: the real classical paths for the final positions
$x_{1}=-2,-1,-0.5,-x_{c}$. Right: the complex classical paths for the final
positions $x_{1}=-0.3,-0.15,0,0.15,0.2$ and the analytically continued paths
for $x_{1}=0.21,0.3,1,2,3,4,5$. The orange and brown paths are conjugate
pairs. The black dot represents the turning point, the background is the
absolute valued squared of the potential and the green points are the
singularities of the potential.

When the energy is below the top of the barrier ($E<V_{0}$) and the initial
and final position lay in the same classically allowed region (either
$x_{0},x_{1}\leq-x_{c}$ or $x_{c}\leq x_{0},x_{1}$ with the turning point
$x_{c}=\cosh^{-1}\sqrt{V_{0}/E}$) there exist two real classical solutions to
the boundary value problem. The red and blue points in fig. 2 represent the
direct and bouncing solution of the boundary value problem (see the left panel
of fig. 2). When keeping the initial position fixed in the left region and
letting the final position $x_{1}$ approach the turning point $x_{c}$, the two
real classical paths coalesce. For larger final positions, no classical
solution exists, as the real paths have formed a conjugate pair of complex
classical paths (see the brown and orange points in the left panel of fig. 2).
According to Picard-Lefschetz theory, the complex classical path for which the
exponent of the path integral $i(S+ET)$ has a negative real part (the orange
path) remains relevant after this transition. The other complex path (marked
by the brown points) is irrelevant to the path integral.

As we increase the final position further, the complex classical path moves
further into the complex plane along the vertical green line, until it
collides with a discontinuity in the final position in the complex $T$-plane.
For larger $x_{1}$, the solution to the boundary value problem ceases to
exist. This dramatic transition occurs as the complex classical path
intersects one of the singularities of the analytic continuation of the
potential at $i\pi(\tfrac{1}{2}+n)$ with $n\in\mathbb{Z}$ (see the right panel
of fig. 2). At this singularity crossing, the naive saddle point approximation
of the path integral Schulman:2012 ,

$\displaystyle
K\sim\sum_{(x_{C},T_{C})}\sqrt{-\frac{\partial^{2}S_{C}/\partial x_{0}\partial
x_{1}}{\partial^{2}S_{C}/\partial T^{2}}}e^{i(S_{C}+ET_{C})/\hbar}\,,$ (14)

ranging over the relevant classical paths fails to approximate the tunneling
behavior of the energy propagator. To continue the saddle point approximation
through this singularity crossing, we interpret the discontinuity (see the
left panel of fig. 2) as a branch cut and analytically continue the saddle
point beyond the space of classical paths satisfying the boundary conditions.
For the Rosen-Morse system, this analytical continuation is most easily
implemented by working with the hyperbolic sine of the final position $\sinh
x_{1}$ (see the right panel of fig. 2). The hyperbolic sine function unveils
the analytic continuation saddle point and highlights the occurrence of a
complex caustic, where the horizontal and vertical green lines intersect away
from the real line, corresponding to $x_{1}$ moving into the classically
allowed region on the right. As far as we are aware, this is the first example
of a system with a complex caustic, where complex saddle points coalesce.
After the complex caustic, the relevant saddle point moves to the left along
the horizontal green line.

Beyond the singularity crossing, the complex $T$ only serves as a label of the
analytically continued saddle point, as the corresponding initial value
problem does not terminate at the boundary condition $x(1)=x_{1}$. Instead,
the solution to the equation of motion approaches $x(1)=-x_{1}-i\pi$ (see the
right panel of fig. 2), as $\sinh(-x_{1}-i\pi)=\sinh(x_{1})$. Note that when
approximating the energy propagator for the configuration $x_{1}=-x_{0}$ for
small $x_{0}$, the solution to the initial value problem linearly interpolates
between $x_{0}$ and $x_{0}-i\pi$. To describe the tunneling behavior, the
classical action also needs to be analytically continued at the singularity
crossing. In the classical action (12), the singularity crossing corresponds
to a branch-cut of the arctanh function, yielding a correction
$i\sqrt{2mV_{0}}\pi$ to the classical action of the solution of the initial
value problem.

We compare the exact energy propagator for a particle below the top of the
barrier in fig. 3. In the classically allowed region, we observe an
interference pattern that can be understood in terms of the direct and
bouncing real classical paths. At the turning point, $x_{1}=x_{c}$, the saddle
point approximation diverges in a fold caustic. For larger $x_{1}$ the real
saddle point approximation vanishes. The complex saddles reasonably
approximate the exact propagator until the singularity crossing, where the
green line drops to zero. After the singularity crossing, the saddle point
approximation with the analytically continued classical path diverges in a
second complex caustic, when $x_{1}$ enters the classically allowed region
after which the approximation approaches the non-zero exact propagator (this
value is proportional to the tunneling amplitude). From this analysis, it is
clear that the singularity crossing and the corresponding analytic
continuation are central to quantum tunneling in the real-time Feynman path
integral. The same can be said for quantum reflections when the energy of the
particle exceeds the potential strength, $E>V_{0}$.

Indeed, as we show in Feldbrugge:2023Rosen , the WKB tunneling rate is only
recovered when including the correction $i\sqrt{2mV_{0}}\pi$ in the saddle
point approximation. Note that the saddle point approximation and the exact
result do not coincide as we are working at finite $\hbar$ and with an energy
close to the top of the barrier. The saddle point approximation converges to
the exact propagator away from the caustics for lower energies and reduced
Planck constants, however, these changes will also suppress the tunneling
amplitude and energy propagator for large $x_{1}$.

Figure 3: A comparison of the energy propagator $E=0.9$ as a funciton of
$x_{1}$ and the saddle point approximations for the initial position
$x_{0}=-5$, for a particle of mass $m=1$ in a barrier of strength $V_{0}=1$
and the reduced Planck constant $\hbar=0.5$. The energy propagator (blue), the
real saddle point approximation (dotted black), the saddle point approximation
with complex solutions to the boundary value problem (dashed green), and the
saddle point approximation including the analytic continuations of the complex
saddle points (red).

The Gaussian barrier: The path integral of the Rosen-Morse potential exhibits
a singularity crossing when a relevant complex classical path intersects a
singularity of the analytically continued action. Similar behavior can be
observed for rational potentials, such as the Lorentzian potential
$V(x)=V_{0}/(1+x^{2})$ with singularities at $x=\pm i$. However, singularity
crossings are not restricted to real-time path integrals of theories whose
analytically continued potential has singularities at finite distances in the
complex plane.

To illustrate this, consider a particle interacting with a Gaussian barrier
$V(x)=V_{0}\exp(-x^{2})$ for $V_{0}>0$, which has an essential singularity at
complex infinity. Let’s consider the real-time path integral

$\displaystyle
G[x_{1},x_{0};T]=\int_{x(0)=0}^{x(1)=x_{1}}e^{iS[x]/\hbar}\mathcal{D}x\,,$
(15)

with the associated Euler-Lagrange equation

$\displaystyle\frac{\delta S}{\delta x}=0:\quad
mx^{\prime\prime}=2V_{0}T^{2}xe^{-x^{2}}\,,$ (16)

with the boundary conditions $x(0)=x_{0}$ and $x(1)=x_{1}$ for fixed real
propagation time $T$. We look for solutions to the boundary value problem by
solving the initial value problem in the space of complex initial velocities
$x^{\prime}(0)=v_{0}$ (see fig. 4).

Figure 4: Classical paths for the Gaussian barrier in the complex initial
velocity plane for the initial position $x_{0}=-5$, particle mass $m=1$, the
barrier strength $V_{0}=1$, and the total propagation time $T=10$. The green
curves indicate initial velocities that lead to a real final position.

When the initial and final positions reside on the same side of the barrier,
there either exists a single or three real classical paths separated by a
caustic. When there are three real classical paths, the particle can either
travel directly between the boundary points, travel with a light bounce, or
spend a significant amount of time near the top of the barrier before rolling
back to the final position. Keeping the initial position fixed on the left of
the barrier, and moving the final position to the top of the barrier, two of
the three real classical paths coalesce in a caustic at the turning point. For
larger $x_{1}$, the real classical paths form complex conjugate pairs of which
one remains relevant to the Feynman path integral. As we see in fig. 4 the
complex path again hits discontinuities of the final position as a function of
the complex initial velocity. As the complex path approaches the discontinuity
(see fig. 4), the classical path travels to complex infinity and back in
finite time (see the paths in fig. 5). The real-time Feynman path integral of
the Gaussian barrier thus includes complex paths undergoing singularity
crossings where the path runs off to complex infinity.

Figure 5: Complex classical paths (orange and green are conjugate paths) for
the Gaussian barrier in the complex initial velocity plane for the initial
position $x_{0}=-5$, particle mass $m=1$, the barrier strength $V_{0}=1$, the
total propagation time $T=10$, and the final position $x_{1}=-1.6,$ $-1,$
$-0.5,$ $-0.25,$ $-0.1,$ $-0.03,$ $-0.02,$ $-0.01$

The discretized path integral: Singularity crossings of complex classical
paths are universal in the continuum path integral but absent in the simplest
discretized version. Consider the time-discretized path integral

$\displaystyle\int e^{\frac{i\epsilon
T}{\hbar}\sum_{j=0}^{n-1}\left[\frac{m}{2}\left(\frac{y_{j+1}-y_{j}}{\epsilon
T}\right)^{2}-V(y_{j})\right]}\mathrm{d}y_{1}\mathrm{d}y_{2}\dots\mathrm{d}y_{n-1}\,,$
(17)

with $n$ equially space timesteps $y_{j}=x(\epsilon j)$ for $j=0,1,\dots,n$ of
size $\epsilon=1/n$ with $y_{0}=x_{0}$ and $y_{n}=x_{1}$. This discretization
resembles a multi-plane lens system in wave optics
Feldbrugge:2023MNRAS.520.2995F . The variation of the action yields the
discrete equation of motion

$\displaystyle\frac{y_{j+1}-y_{j}}{\epsilon T}=\frac{y_{j}-y_{j-1}}{\epsilon
T}-\frac{\epsilon T}{m}V^{\prime}(y_{i})\,,$ (18)

for $j=1,2,\dots,n-1$. This set of equations has a unique recursive solution
for $y_{2},\dots,y_{n}$ given the initial position $y_{0}$ and the discrete
initial velocity $\bar{v}_{0}=(y_{1}-y_{0})/\epsilon$ resembling the initial
value problem,

$\displaystyle y_{j}$
$\displaystyle=x_{0}+j\epsilon\bar{v}_{0}-\frac{\epsilon^{2}T^{2}}{m}\sum_{k=1}^{j-1}(j-k)V^{\prime}(y_{k})\,,$
(19)

for $j=2,\dots,n$. When the potential is analytic and free of branch cuts,
such as $V(x)=\text{sech}^{2}(x)$, so is the final position $y_{n}$ as a
function of the initial position and discretized initial velocity.
Consequently, the solutions of the discretized Rosen-Morse system never
undergo singularity crossings, making it qualitatively different from the
continuum theory.

While the discretized theory does not have singularity crossings, the final
position $y_{n}$ does develop a large set of singularities in the complex
$\bar{v}_{0}$ plane and a large set of additional complex solutions to the
boundary value problem, rapidly increasing with the number of discretization
steps. When tracking a singularity crossing event in the continuum theory, we
find that the corresponding discretized path develops large jumps,
invalidating the discretization condition. Further research is required to
verify whether the absence of singularity crossings in the discretized theory
is replaced by additional relevant complex discretized classical paths.

Implications: The real-time Feynman path integral studied in terms of
classical paths not only includes real and complex but also analytically
continued paths. Beyond the singularity crossing, the complex classical path
no longer solves the boundary value problem, yet can dominate the path
integral. We believe that previous attempts to express classically forbidden
phenomena with complex paths have been frustrated by the lack of this
observation. When interpreting such an analytically continued path in terms of
an initial value problem, the termination point is altered and additional
terms need to be added to the classical action and functional determinant for
the saddle point approximation to hold.

This phenomenon has large implications in the study of instantons,
particularly in quantum gravity where complex solutions of the Einstein field
equations are commonly studied Coleman:1980 ; Hawking:1982 ; Hartle:1983 ;
Vilenkin:1982 ; Feldbrugge:2017 . Since general relativity predicts curvature
singularities, complex solutions will likely undergo singularity crossings
affecting the path integral for gravity. Witten recently proposed a complex
space-time selection criterion Witten:2021  based on the study of quantum
fields on complex metrics Louko:1997 ; Kontsevich:2021  which will need to be
reconsidered in light of these singularity crossings; how does one embed a
quantum field on an analytically continued spacetime, when the solution to the
Einstein field equations satisfying the boundary conditions does not exist?

## Acknowledgements

The work of JF is supported by the STFC Consolidated Grant ‘Particle Physics
at the Higgs Centre,’ and, respectively, by a Higgs Fellowship and the Higgs
Chair of Theoretical Physics at the University of Edinburgh.

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) license to any Author Accepted Manuscript version arising
from this submission.

## References

  * (1) R. P. Feynman, “Space-Time Approach to Non-Relativistic Quantum Mechanics,” Reviews of Modern Physics 20 no. 2, (Apr., 1948) 367–387.
  * (2) R. P. Feynman and A. R. Hibbs, Quantum mechanics and path integrals.  International series in pure and applied physics. McGraw-Hill, New York, NY, 1965.
  * (3) R. P. Feynman, “Space-Time Approach to Non-Relativistic Quantum Mechanics,” Reviews of Modern Physics 20 no. 2, (Apr., 1948) 367–387.
  * (4) N. Turok, “On quantum tunneling in real time,” New Journal of Physics 16 no. 6, (June, 2014) 063006, arXiv:1312.1772 [quant-ph].
  * (5) J. Feldbrugge and N. Turok, “Existence of real time quantum path integrals,” Annals of Physics 454 (July, 2023) 169315, arXiv:2207.12798 [hep-th].
  * (6) L. Schulman, Techniques and Applications of Path Integration.  Dover Books on Physics. Dover Publications, 2012.
  * (7) N. Rosen and P. M. Morse, “On the Vibrations of Polyatomic Molecules,” Physical Review 42 no. 2, (Oct., 1932) 210–217.
  * (8) G. Pöschl and E. Teller, “Bemerkungen zur Quantenmechanik des anharmonischen Oszillators,” Zeitschrift fur Physik 83 no. 3-4, (Mar., 1933) 143–151.
  * (9) H. Kleinert and I. Mustapic, “Summing the spectral representations of Pöschl-Teller and Rosen-Morse fixed-energy amplitudes,” Journal of Mathematical Physics 33 no. 2, (Feb., 1992) 643–662.
  * (10) J. Feldbrugge, D. Jow, and U.-L. Pen, “Complex classical paths in quantum reflections and tunneling,” (forthcoming) .
  * (11) J. Feldbrugge, “Multiplane lensing in wave optics,” Monthly Notices of the RAS 520 no. 2, (Apr., 2023) 2995–3006.
  * (12) S. Coleman and F. de Luccia, “Gravitational effects on and of vacuum decay,” Phys. Rev. D 21 no. 12, (June, 1980) 3305–3315.
  * (13) S. W. Hawking, “The boundary conditions of the universe,” in Astrophysical Cosmology Proceedings, pp. 563–572.  Jan., 1982.
  * (14) J. B. Hartle and S. W. Hawking, “Wave function of the Universe,” Phys. Rev. D 28 no. 12, (Dec., 1983) 2960–2975.
  * (15) A. Vilenkin, “Creation of universes from nothing,” Physics Letters B 117 no. 1-2, (Nov., 1982) 25–28.
  * (16) J. Feldbrugge, J.-L. Lehners, and N. Turok, “Lorentzian quantum cosmology,” Phys. Rev. D 95 no. 10, (May, 2017) 103508, arXiv:1703.02076 [hep-th].
  * (17) E. Witten, “A Note On Complex Spacetime Metrics,” arXiv e-prints (Nov., 2021) arXiv:2111.06514, arXiv:2111.06514 [hep-th].
  * (18) J. Louko and R. D. Sorkin, “Complex actions in two-dimensional topology change,” Classical and Quantum Gravity 14 no. 1, (Jan., 1997) 179–203, arXiv:gr-qc/9511023 [gr-qc].
  * (19) M. Kontsevich and G. Segal, “Wick rotation and the positivity of energy in quantum field theory,” arXiv e-prints (May, 2021) arXiv:2105.10161, arXiv:2105.10161 [hep-th].

##  I Supplemenatary material

While the real solutions of the time-discretized theory

$\displaystyle\int e^{\frac{i\epsilon
T}{\hbar}\sum_{j=0}^{n-1}\left[\frac{m}{2}\left(\frac{y_{j+1}-y_{j}}{\epsilon
T}\right)^{2}-V(y_{j})\right]}\mathrm{d}y_{1}\mathrm{d}y_{2}\dots\mathrm{d}y_{n-1}\,,$
(20)

approach the continuum theory

$\displaystyle
G[x_{1},x_{0};T]=\int_{x(0)=x_{0}}^{x(1)=x_{1}}e^{\frac{i}{\hbar}\int_{0}^{1}\left[\frac{m\dot{x}^{2}}{2T}-TV(x)\right]\mathrm{d}\lambda}\mathcal{D}x$
(21)

as the discretization is refined, the complex structure of the solutions is
dramatically different (see fig. 6).

In the continuum theory, the final position as a function of complex initial
velocity develops a set of branch points and branch cuts (see fig. 6(a)).
Generally, there exists an infinite set of complex solutions to the boundary
value problem roughly homogeneously distributed in the complex plane. We find
that only a select few of these complex paths are relevant to the Feynman path
integral Feldbrugge:2023Rosen .

The final position in the discretized classical paths diverges for several
complex discretized initial velocities $\bar{v}_{0}$. Note that the number of
singularities increases as we increase the number of discretization steps $n$.
Each singularity comes with a set of complex solutions to the boundary value
problem, corresponding with the intersections of the red and green curves in
fig. 6. It is an open problem which of these classical paths is relevant to
the discretized path integral. One would hope that the complex classical path
resembling the relevant continuous complex path dominates the discretized path
integral.

However, this is probably not the full story, as the singularity crossings
present in the continuum theory are absent in the discretized version. In
particular, as we track the discretized classical path closest to the relevant
continuum complex classical path while increasing the final position $x_{1}$,
we find that the discretized classical path develops large jumps when the
continuum classical path crosses the branch cut (see fig. 7). The path
contains a gradual evolution followed by a large jump to the final position.
These phenomena need to be studied in more detail in the future, as large
jumps normally signal places where the discretization assumption is no longer
valid. However, it could be that a good approximation of the continuum path
integral is recovered with an additional set of discretized complex classical
paths that are absent in the continuum theory.

(a) Continuum

(b) $n=2$

(c) $n=3$

(d) $n=4$

(e) $n=5$

(f) $n=6$

Figure 6: Comparison of the continuum and discretized classical paths. We plot
the final position $x_{1}$ in the complex (continuum and discretized) initial
velocity plane. We consider classical paths of a particle with mass $m=1$
starting at $x_{0}=-5$ propagating for a time $T=10$ while interacting with a
Rosen-Morse barrier with strength $V_{0}=1$. The red curves correspond to the
classical paths terminating at a real position. The orange curves correspond
to their analytic continuations beyond the singularity crossings. The green
lines represent the classical paths that terminate at a position whose real
part is $-1$. The intersections of the red ($Im[x_{1}]=0$) and green lines
($Re[x_{1}]=-1$) are classical paths propagating between $x_{0}=-5$ and
$x_{1}=-1$. Figure 7: Discretized complex classical path for $n=6$ as the
final position crosses past the singularity crossing of the continuum theory.