On solving Schrödinger’s equation with classical action

Winfried Lohmiller and Jean-Jacques Slotine

Nonlinear Systems Laboratory

Massachusetts Institute of Technology

Cambridge, Massachusetts, 02139, USA

{wslohmil<EMAIL_ADDRESS>
###### Abstract

We show that the Schrödinger equation of quantum physics can be solved using
the classical Hamilton-Jacobi action dynamics, extending a key result of
Feynman applicable only to quadratic Lagrangians. This is made possible by two
developments. The first is incorporating geometric constraints directly in the
classical least action problem, in effect replacing in part the probabilistic
setting by the non-uniqueness of solutions of the constrained problem. For
instance, in the double slit experiment or for a particle in a box, spatial
inequality constraints create Dirac constraint forces, which lead to multiple
path solutions. The second development is a spatial rescaling of clocks,
specifically designed to achieve a general equivalence between Schrödinger and
Hamilton-Jacobi representations. These developments leave the results of
associated Feynman path integrals unchanged, but the computation can be
greatly simplified as only classical action is used and time-slicing is
avoided altogether. They also suggest a smooth transition between physics
across scales.

##  1 Introduction

Attempts to bridge the conceptual gap between classical and quantum physics
have a long and very distinguished history. Central among those is the path-
integral formulation of quantum mechanics, starting with Wiener’s work on
stochastic processes, Dirac’s discussion of the relation of classical least
action to quantum mechanics, Feynman’s fundamental paper [7, 8] on path
integral computation, and more recent important extensions such as Duru and
Kleinert’s time reparametrization [4, 10, 19, 20].

The current paper stems from the same general motivation. It starts by
deriving simple results on classical action optimization of Lagrangian
dynamics subject to spatial inequality constraints. Such constraints typically
imply multiple or set-valued [17] least action solutions of the optimization
problem. It then shows that quantum kernels, computed with the usual Feynman
integral along all possible zig-zag paths [8] and its many extensions [10],
can be computed instead using classical action, without requiring the
Lagrangian to be quadratic in velocity. In the double slit or $n$-slit
experiments, for instance, the spatial inequality constraints simply represent
the geometry of the slits and the multiple action solutions correspond to
diffraction behind the slits. The paper then derives more generally an
equivalence between Schrödinger and Hamilton-Jacobi representations, based on
a state transformation and a rescaling of time. Besides its computational
implications, this equivalence also suggests a smooth transition between
physics models across scales.

Recall that the classical motion is a local extremum of the action

$\phi({\bf q},t)\ =\int_{0}^{t}Ldt\ +\ \phi({\bf q},0)\ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ L=\frac{1}{2}\ \dot{\bf q}^{T}{\bf M}({\bf q})\dot{\bf q}-V({\bf
q},t)-E$ (1)

with ${\bf q}\in\mathbb{R}^{n}$ (see e.g. [8, 15]), inertia tensor ${\bf M}$,
potential energy $V$ with arbitrary constant offset $E$.

Note that the Lagrangian $L$ can in principle be extended with a linear term
in $\dot{\bf q}\ $, e.g. to describe electromagnetic effects, without changing
the main results of this paper. The action $\phi({\bf q},t)$ can be computed
from the Hamilton-Jacobi equation

$\ \ \ \ {\frac{\partial\phi}{\partial t}}+h=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ h={\frac{1}{2}}\ \frac{\partial\phi}{\partial{\bf q}}^{T}\ {\bf M}^{-1}({\bf
q})\frac{\partial\phi}{\partial{\bf q}}+V({\bf q},t)+E$ (2)

For Lipschitz continuous ${\bf M}({\bf q})$ and $V({\bf q},t),$ and with ${\bf
q}\in\mathbb{R}^{n}$, the Hamilton-Jacobi equation has a unique [17] solution
[6] for a given initial $\phi({\bf q},0)$. However, this deterministic unique
solution can change to multiple solutions if spatial constraints are
introduced. Also, note that the constant offset $\ E\ $ modifies $\ \phi\ $
but does not affect the related momentum and position dynamics.

Section 2 first extends the standard Euler-Lagrange or Hamiltonian dynamics
with ${\bf q}\in\mathbb{R}^{n}$ to the case of constrained positions ${\bf
q}\in\mathbb{G}^{n}\subset\mathbb{R}^{n}$, with $\mathbb{G}^{n}$ defined by
$k=1,...,K\ $ inequality constraints

$g_{k}({\bf q},t)\leq 0$ (3)

At the border $\partial\mathbb{G}^{n}$ of $\mathbb{G}^{n}$, a Dirac constraint
force ensures that the constraint is not violated. This non-Lipschitz
activation of the constraint force can lead to multiple [17], but
deterministic, action and path solutions.

Section 3 then shows that the wave function $\Psi({\bf q},t)$ of the time-
dependent Schrödinger equation [3, 11]

$0=\left[-i\hbar{\frac{\partial}{\partial
t}}-\frac{\hbar^{2}}{2}\Delta_{M}+V({\bf q},t)\right]\Psi({\bf q},t)$ (4)

with the Laplacian (Laplace-Beltrami operator) [15] in a metric ${\bf M}({\bf
q})$

$\displaystyle\Delta_{M}$ $\displaystyle=$
$\displaystyle{\frac{1}{\sqrt{\det{\bf M}}}}\frac{\partial}{\partial{\bf
q}}^{T}\left({\sqrt{\det{\bf M}}}\ {\bf M}^{-1}\frac{\partial}{\partial{\bf
q}}\right)$ (5)

can be solved using only classical action, where inequality constraints (3)
may generate multiple solutions.

To handle more general applications, the spatial scaling from the metric ${\bf
M}({\bf q})\ $ is extended with a temporal scaling $T({\bf q},t)$, inspired by
[4, 10], and used to cancel a specific term in the associated Schrödinger
equation. The Hamilton-Jacobi equation, the Schrödinger equation, and the
mapping between them which we introduce, are all coordinate-invariant tensor
results. We also show how contraction analysis [12] may be used to quantify
the divergence of classical paths during diffraction or measurement collapse.

Let us recall on an example how Feynman [8] solves the Schrödinger equation
(4) using the classical motion (2), assuming a quadratic Lagrangian and
without any constraint.

  * ###### Example 1.1

:  Consider the $n$-dimensional harmonic oscillator

$h=\frac{1}{2M}\frac{\partial\phi}{\partial{\bf
q}}^{T}\frac{\partial\phi}{\partial{\bf q}}+\frac{M\omega^{2}}{2}{\bf
q}^{T}{\bf q}$

with angular frequency $\omega$ and constant mass $M$. The classical action
(2) connecting an initial point ${\bf q}_{o},0$ to a final point ${\bf q},t$
is

$\phi=\frac{M\omega}{2\sin{\omega t}}\left(({\bf q}^{T}{\bf q}+{\bf
q}_{o}^{T}{\bf q}_{o})\cos{\omega t}-2{\bf q}^{T}{\bf q}_{o}\right)$

starting as a Dirac impulse at $t=0$. The exponential of this action solves
the time-dependent Schrödinger equation (4) as [8]

$\displaystyle\Psi$ $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}\phi-\int\frac{\Delta_{M}\phi}{2}dt}=e^{\frac{i}{\hbar}\phi-\int\frac{nM\omega}{2}\cot{\omega
t}dt}=\sqrt{\frac{M\omega}{2\pi i\hbar\sin{\omega t}}}^{n}\
e^{\frac{i}{\hbar}\phi}$ $\displaystyle=$
$\displaystyle\sum_{j_{1},...j_{n}=0}^{\infty}e^{-\frac{i}{\hbar}E_{j}t}\
\psi_{j_{1},...j_{n}}({\bf q})\ \psi_{j_{1},...j_{n}}({\bf q}_{o})$

where the second equation algebraically expands the wave function $\Psi$ in
powers $e^{-\frac{i}{\hbar}E_{j}t}$ of the eigenvalues
$E_{j}=\hbar\omega(j+\frac{1}{2}),j=j_{1}+...+j_{n}$ using Euler’s formula [4,
8]. The corresponding eigenfunctions are

$\psi_{j_{1},...j_{n}}=\sqrt[4]{\frac{M\omega}{\pi\hbar}}^{n}\frac{1}{\sqrt{2^{(j_{1}+...+j_{n})}j_{1}!...j_{n}!}}H_{j_{1}}...H_{j_{n}}e^{-\frac{M\omega}{2\hbar}{\bf
q}^{T}{\bf q}}$ (6)

with the Hermite polynomials $\ H_{j}(q_{j}\sqrt{\frac{M\omega}{\hbar}})\ $
[4, 8]. $\Box$

The goal of this paper is to extend this conversion, from classical action to
the Schrödinger wave function, to constrained positions ${\bf
q}\in\mathbb{G}^{n}$ and non-quadratic Lagrangians. This conversion will be
illustrated for the double and single slit experiment, particle in a box and
Coulomb potential in section 4.

##  2 Constrained Hamilton-Jacobi Dynamics

In this section, we extend the standard Euler-Lagrange or Hamilton-Jacobi
dynamic formulations to the case of spatial inequality constraints. First, let
us denote whether spatial inequality constraints (3) are active or not.

###### Definition 1

The constrained space $\mathbb{G}^{n}$ is defined by the $k=1,...,K$
inequality constraints

$g_{k}({\bf q},t)\leq 0$

The set of active constraints $\mathcal{A}({\bf q},t)\subset\\{1,...,K\\}$ is
the set of indices $j$ on the boundary $\partial\mathbb{G}^{n}$ of
$\mathbb{G}^{n}$, i.e., such that

$g_{k}({\bf q},t)=0$

The action $\phi$ (1) has a local extremum if the variation of the action (1)

$\displaystyle\delta\phi$ $\displaystyle=$
$\displaystyle\int_{0}^{t}\frac{\partial L}{\partial\dot{\bf q}}\delta\dot{\bf
q}+\frac{\partial L}{\partial{\bf q}}\delta{\bf q}dt=\left[\frac{\partial
L}{\partial\dot{\bf q}}\delta{\bf
q}\right]_{0}^{t}-\int_{0}^{t}\left[\frac{d}{dt}\frac{\partial
L}{\partial\dot{\bf q}}-\frac{\partial L}{\partial{\bf q}}\right]\delta{\bf
q}dt=\int_{0}^{t}\sum_{k\in\mathcal{A}}\lambda_{k}\frac{\partial
g_{k}}{\partial{\bf q}}\delta{\bf q}dt$

is only non-zero orthogonal to an active constraint, where the Lagrange
parameter $\lambda_{j}$ defines the magnitude of the cost gradient at the
active constraint.

The first term on the right-hand side is zero since $\delta{\bf q}$ is zero at
the start and end points. Between the end points $\delta{\bf q}$ can take on
any arbitrary value. Thus a least action solution satisfies

$\frac{d}{dt}\frac{\partial L}{\partial\dot{\bf q}}-\frac{\partial
L}{\partial{\bf q}}\ =\ \frac{d}{dt}\frac{\partial\phi}{\partial{\bf
q}}+\frac{\partial h}{\partial{\bf q}}\ =\
\sum_{k\in\mathcal{A}}\lambda_{k}\frac{\partial g_{k}}{\partial{\bf q}}$ (7)

This extends the usual Euler-Lagrange or Hamiltonian dynamics (see e.g. [15])
with Lagrangian collision forces activated by inequality constraints. While
constrained dynamics with equality constraints [16] would simply correspond to
Lagrange’s method of the first kind, inequality constraints (3) are
fundamentally different since they lead to Dirac constraint forces at the
collision with the border $\partial\mathbb{G}^{n}$. At the border ${\bf
q}\in\partial\mathbb{G}^{n}$ multiple solutions
$\frac{\partial\phi}{\partial{\bf q}}({\bf q},t,j),j\in\mathbb{Z}$ can occur,
which are continuous in ${\phi}({\bf q},t,j)$ due to (1). By contrast, for
Lipschitz continuous ${\bf M}({\bf q}),V({\bf q},t)$ and for ${\bf
q}\in\mathbb{G}^{n}\setminus\partial\mathbb{G}^{n}$ equation (7) has a unique
integrable solution $\frac{\partial\phi}{\partial{\bf q}}({\bf q},t,j)$ of the
Hamilton-Jacobi equation (2) [6] for a given $j\in\mathbb{Z}$ at the boundary
$\frac{\partial\phi}{\partial{\bf q}}({\bf q},t,j)$.

The constraint forces can be computed by assuring that each active constraint
$g_{\ell}({\bf q},t)=0,\dot{g}_{\ell}=\frac{\partial{g}_{\ell}}{\partial{\bf
q}}^{T}\dot{\bf q}+\frac{\partial{g}_{\ell}}{\partial 0}>0\ $ is not violated
at $\ t+dt$, i.e. for an instantaneous collision

$\ddot{g}_{\ell}=\frac{\partial{g}_{\ell}}{\partial{\bf q}}{\bf
M}^{-1}\left(\sum_{k\in\mathcal{A}}\lambda_{k}\frac{\partial
g_{k}}{\partial{\bf q}}\right)^{T}\leq\delta(g_{\ell})\dot{g}_{\ell}$ (8)

where $\delta(g_{\ell})$ is the Dirac delta function. The equation above
represents a partial elastic collision and can e.g. be solved with linear
programming.

Before we come to the main theorem, let us introduce d’Alembert’s time scaling
[4, 10]

$\tau=\int_{0}^{t}T({\bf q},t,j)dt\iff t=\int_{0}^{\tau}\frac{1}{T({\bf
q},t,j)}d\tau$ (9)

where the time integral is taken along the curve ${\bf q}(t)$ and the time is
invertible for bounded $T>0$. This enables in many cases the analytic
computation of the action of nonlinear Lagrangians as we will see in Lemma 1
and section 4. Note that $T({\bf q},t,j)$ is held constant at a gradient
conversion from action to momentum, since it only scales the dynamics along
and not across trajectories.

Let us now summarize the results above:

###### Theorem 1

The action (1) constrained according to Definition 1 in the scaled time (9) is
locally extremized by the momentum dynamics

$\displaystyle T\frac{d}{d\tau}\frac{\partial\phi}{\partial{\bf q}}\ +\
\frac{\partial h}{\partial{\bf q}}\ $ $\displaystyle=$
$\displaystyle\sum_{all\ k\in\mathcal{A}}\frac{\partial g_{k}}{\partial{\bf
q}}\lambda_{k}$ (10) $\displaystyle{\bf M}({\bf q})\ T\ \frac{d{\bf
q}}{d\tau}$ $\displaystyle=$ $\displaystyle\frac{\partial\phi}{\partial{\bf
q}}$

where the partially elastic collision force $\lambda_{k}\ $ fulfills (8). A
fully elastic collision force would in addition leave $h$ constant at the
collision instant.

  * •

At the border ${\bf q}\in\partial\mathbb{G}^{n}$, multiple momentum solutions
$\frac{\partial\phi}{\partial{\bf q}}({\bf q},t,j)$ with
$j\in\mathbb{J}\subset\mathbb{Z}$ can occur, leading to a set-valued and
continuous action function $\phi({\bf q},t,j)$.

  * •

Outside the border ${\bf q}\in\mathbb{G}^{n}\setminus\partial\mathbb{G}^{n}$
and for Lipschitz continuous ${\bf M}({\bf q}),\frac{V}{T}({\bf q},t)$ the
number of momentum $\frac{\partial\phi}{\partial{\bf q}}({\bf q},t,j)$ and
action ${\phi}({\bf q},t,j)$ solutions of the time scaled Hamilton-Jacobi
equation

$T\ {\frac{\partial\phi}{\partial\tau}}+h=0$ (11)

does not increase.

Note for any given $j\in\mathbb{J}$, the action ${\phi}({\bf q},t,j)$ is an
ordinary single-valued function.

The original formulation by Hamilton, Lagrange and Jacobi was derived in
$\mathbb{R}^{n}$, but not in $\mathbb{G}^{n}$, and thus cannot be used to
predict multiple solutions at a constraint. The computation of elastic
collision forces is well established [16], but to the authors’ knowledge has
not been used in the context of constrained momentum (10) and Hamiltonian
dynamics (11).

The time scaling (9) may be exploited to find an analytic solution for the
action $\phi$. For instance, by analogy with the harmonic oscillator of
Example 6 and [4], one can turn a nonlinear Hamiltonian (2) in $n$-dimensional
position coordinates ${\bf q}({\bf q}^{\prime})$

$\displaystyle\frac{h}{T}\ -\ \frac{V({\bf q},t)}{T}-\ \frac{E}{T}$
$\displaystyle=$ $\displaystyle\frac{1}{2T}\frac{\partial\phi}{\partial{\bf
q}}^{T}{\bf M}^{-1}({\bf q})\frac{\partial\phi}{\partial{\bf q}}$
$\displaystyle T({\bf q},t)$ $\displaystyle=$
$\displaystyle\frac{V}{C}=\frac{1}{{\bf q}^{\prime T}{\bf q}^{\prime}}$ (12)

with constant $M,C$, into a quadratic Hamiltonian

$\frac{h}{T}=\frac{1}{2M}\frac{\partial\phi}{\partial{\bf
q}^{\prime}}^{T}\frac{\partial\phi}{\partial{\bf q}^{\prime}}+C-E\ {\bf
q}^{\prime T}{\bf q}^{\prime}$

in $n+1$-dimensional coordinates ${\bf q}^{\prime}$ defined by

$\displaystyle{\bf M}^{-1}({\bf q})$ $\displaystyle=$
$\displaystyle\frac{T}{M}\frac{\partial{\bf q}}{\partial{\bf
q}^{\prime}}\frac{\partial{\bf q}}{\partial{\bf q}^{\prime}}^{T}$ (13)

The existence of the solution ${\bf q}^{\prime}$ of the p.d.e. above is given
by the Nash-Kuiper theorem [18]. Based on Example 6, this yields the following
construction.

###### Lemma 1

The time scaled Hamilton-Jacobi equation (11) of the Hamiltonian (12) in the
position coordinates ${\bf q}^{\prime}$, defined by (13) with
$E=-\frac{M}{2}\omega^{2}$, can be solved explicitly for $\
-\infty<\tau<\infty\ $ by the analytic action

$\phi({\bf q}^{\prime},\tau,j)=\frac{M\omega}{2\sin{\omega\tau}}\left(({\bf
q}^{{}^{\prime}T}{\bf q}^{\prime}+{\bf q}_{o}^{{}^{\prime}T}{\bf
q}^{\prime}_{o})\cos{\omega\tau}-2{\bf q}^{{}^{\prime}T}{\bf
q}^{\prime}_{o}\right)-C\tau$ (14)

where $\ j\in\mathbb{J}\ $. This gives the Laplacian $\ \Delta_{M}\phi=(n+1)\
\omega\cot{\omega\tau}\ $.

##  3 Constrained Hamiltonian Solution of the Schrödinger equation

We now build on the above result to show the equivalence of the Hamilton-
Jacobi and Schrödinger formulations, using the Schrödinger equation (4) at a
scaled time $\tau$ (9) [4]

$0=\left[-i\hbar
T{\frac{\partial}{\partial\tau}}-\frac{\hbar^{2}}{2}\Delta_{M}+V({\bf
q},t)\right]\Psi({\bf q},t)$ (15)

The Hamilton equation, as used at the times of Schrödinger and Feynman, only
had a single deterministic solution $\phi$ without diffraction at a constraint
$\partial\mathbb{G}^{n}$. Hence the Feynman path integral [8] had to consider
all possible stochastic zig-zag trajectories with a time slicing approach,
rather than just those extremizing (1). This stochastic complexity can be
avoided if one uses the multiple local extremal action solutions of Theorem 1.
The Schrödinger equation (15) can now be solved with the set-valued
exponential function of the set-valued action

$\displaystyle\psi({\bf q},t,j)\ $ $\displaystyle=$ $\displaystyle\
e^{\frac{i}{\hbar}\phi({\bf
q},t,j)-\int\frac{\Delta_{M}\phi}{2T}(\tau,j)d\tau+\frac{iE\tau}{\hbar}}\ \ \
\ \ \ \ \ \rm{with}\ \ j\in\mathbb{J}$ (16)

where the temporal scaling $T({\bf q},t,j)$ (9) was selected to assure that
$\frac{\Delta_{M}\phi}{2T}$ is independent of ${\bf q}$. Replacing the wave
function $\Psi$ in (15) with (16) leads to the equivalence of the Hamilton-
Jacobi equation (11)

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\left[-i\hbar
T{\frac{\partial}{\partial{\tau}}}-\frac{\hbar^{2}}{2}{\frac{1}{\sqrt{\det{\bf
M}}}}\frac{\partial}{\partial{\bf q}}^{T}\left({\sqrt{\det{\bf M}}}\ {\bf
M}^{-1}\frac{\partial}{\partial{\bf q}}\right)+V\right]\psi$ $\displaystyle=$
$\displaystyle\left[T{\frac{\partial\phi}{\partial{\tau}}}+{\frac{1}{2}}\
\frac{\partial\phi}{\partial{\bf q}}^{T}\ {\bf
M}^{-1}\frac{\partial\phi}{\partial{\bf q}}+V+E\right]\psi$

Note that the first equation is an operator equation, which becomes a product
in the second equation thanks to the exponential form of (16).

###### Theorem 2

Using the mapping (16) and selecting a temporal scaling (9) such that $\
\frac{\Delta_{M}\phi}{T}(\tau,j)$ is independent of ${\bf q}$, the set-valued
wave function $\psi({\bf q},t,j),j\in\mathbb{J}$ of the Schrödinger equation
(15) is equivalent to the set-valued action $\phi({\bf q},t,j),j\in\mathbb{J}$
(11) of momentum dynamics of Theorem 1.

The overall single-valued wave function solution is the equally weighted,
normalized superposition

$\Psi({\bf q},\tau)=\langle\Psi({\bf q},\tau)|\Psi({\bf
q}_{o},0)\rangle=\frac{\Psi^{\prime}}{\int_{\mathbb{G}^{n}}\sqrt{\Psi^{\prime}\Psi^{\prime\ast}}\
d{\bf q}}\ \ \ \ \ \ \ \ \ \ \Psi^{\prime}=\sum_{j\in\mathbb{J}}\psi({\bf
q},t,j)$ (17)

Equations (16, 17 ) asssume an initial Dirac distribution $\psi({\bf
q}_{o},0)=\delta({\bf q}_{o})$. Since (15) is linear, it extends immediately
to any arbitrary initial normalized distribution $\Psi({\bf q}_{o},0)$

$\Psi({\bf
q},t)=\int_{-\infty}^{\infty}...\int_{-\infty}^{\infty}\langle\Psi({\bf
q},t)|\Psi({\bf q}_{o},0)\rangle\Psi({\bf q}_{o},0)dq_{1},...,dq_{n}$ (18)

Equations (16) and (18) extend the non zig-zag Feynman path integral [8],
taken along the classical extremal action solution for a quadratic Lagrangian,
to

  * •

a constrained space $\mathbb{G}^{n}$ rather then $\mathbb{R}^{n}$, as the
constrained extremal action ${\phi}$ of Theorem 1 already includes diffraction
in ${\bf q}$.

  * •

a Laplacian $\frac{\Delta_{M}\phi}{T}(\tau,j)$ which only has to be
independent of ${\bf q}$, rather than being required to be quadratic.

The general Feynman path integral with zig-zag paths and time-slicing [4, 8,
10] is of course applicable to constrained space and non-quadratic
Lagrangians. However, the same results are obtained here using only the set-
valued action and classical paths (10), which are much fewer and are
continuous everywhere.

Contraction Interpretation

Let us briefly discuss the above results from a contraction analysis
perspective [12], which will allow to quantify analytically the divergence of
classical paths in the context e.g. of diffraction or measurement collapse
phenomena.

The Laplace-Beltrami operator $\Delta_{M}$ of (5) can be derived from the more
general second covariant derivative with respect to the metric ${\bf M}({\bf
q})$, using the related Christoffel term [15],

$\displaystyle\Delta_{M}\ $ $\displaystyle=$ $\displaystyle tr({\bf
M}^{-1}\nabla_{M}^{2})$ (19) $\displaystyle\nabla_{M}^{2}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}}{\partial q_{i}\partial
q_{j}}-\sum_{k=1}^{n}\gamma_{ij}^{k}\ \frac{\partial}{\partial q_{k}}\ \ \ \
\rm{where}\ \ \ \ \gamma_{ij}^{k}({\bf
x},t)=\frac{1}{2}\sum_{l=1}^{n}\left(\frac{\partial M_{il}}{\partial
q^{j}}+\frac{\partial M_{jl}}{\partial q^{i}}-\frac{\partial M_{ij}}{\partial
q^{l}}\right)M_{kl}^{-1}$

Earlier work [12, 13] shows that the contraction rate between two neighbouring
trajectories of the position dynamics $\ {\bf M}({\bf q})\dot{\bf
q}=\frac{\partial\phi}{\partial{\bf q}}\ $ is then given by

$\displaystyle\left({\bf\delta q}^{T}{\bf M}\ {\bf\delta
q}\right)\lambda_{\min}\ \leq\ \frac{1}{2}\ \frac{d}{dt}\left({\bf\delta
q}^{T}{\bf M}\ {\bf\delta q}\right)\ =\ {\bf\delta q}^{T}\nabla_{M}^{2}\phi\
{\bf\delta q}\ \leq\ \left({\bf\delta q}^{T}{\bf M}\ {\bf\delta
q}\right)\lambda_{\max}$

where $\lambda$ is the generalized eigenvalue of $\nabla_{M}^{2}\phi$ with
respect to ${\bf M}$, ${\rm det}(\nabla_{M}^{2}\phi-\lambda{\bf M})=0$.
Therefore, the geodesic distance $s=\min\int_{{\bf q}_{1}}^{{\bf
q}_{2}}\sqrt{\delta{\bf q}^{T}{\bf M}\delta{\bf q}}\ $ between two arbitrary
trajectories ${\bf q}_{1}(t)$ and ${\bf q}_{2}(t)$ verifies

$s(0)e^{\ \int\limits_{o}^{t}\lambda_{\min}({\bf q},t)dt}\leq s(t)\leq s(0)\
e^{\ \int\limits_{o}^{t}\lambda_{\max}({\bf q},t)dt}$ (20)

Since $\nabla_{M}^{2}\phi$ is a tensor [15], the contraction rates $\lambda$
are invariant under any spatial coordinate transformation ${\bf
q}^{\prime}({\bf q})$. This also holds for $\Delta_{M}\phi({\bf
q})=\Delta_{M}\phi({\bf q}^{\prime})$.

In addition the Schrödinger equation is indifferent, and hence the
differential variation between two neighbouring wave functions in (15) is
preserved [14]. Using the mapping (16),

$0\ =\ {\frac{\partial}{\partial\tau}}\int_{\mathbb{G}^{n}}|\delta\Psi|^{2}\
dq_{1}...dq_{n}\ =\
{\frac{\partial}{\partial\tau}}\int_{\mathbb{G}^{n}}e^{-\int\Delta_{M}\phi({\bf
q},\tau,j)d\tau}\ \delta\phi^{2}\ dq_{1}...dq_{n}$ (21)

we see that the weighted action variation above is also preserved. Computing
the corresponding geodesic distance between two arbitrary solutions then turns
this local stability result in a global result [14]. Thus, the Schrödinger
(15) and the Hamilton-Jacobi (11) equations are both indifferent.

##  4 Examples of Constrained Hamiltonian Solutions of the Schrödinger
equation

The above results are now illustrated for the double slit, single slit,
particle in a box, and Coulomb potential examples, where we compute the wave
function from the action. In all applicable examples the Lagrangian collision
force (8) simply corresponds to the momentum change at the collision. Finally,
note that in principle all examples could also be inverted to compute the
action from the wave function.

  * ###### Example 4.1

:  Consider the double slit experiment in Figure 1, with Hamiltonian
$h=\frac{1}{2M}\frac{\partial\phi}{\partial{\bf
q}}^{T}\frac{\partial\phi}{\partial{\bf q}}$ and constant mass $M$. We define
the real action on ${\bf
q}=(q_{1},q_{2})^{T}\in\mathbb{G}^{2}=\mathbb{R}^{2}\setminus\mathbb{E}^{2}$
rather than on $\mathbb{R}^{2}$, thus excluding the red double slit wall
obstacle $\mathbb{E}^{2}$ in Figure 1. As in [8], we assume the two slots are
infinitesimally thin, and consider the two cases of a given initial position
${\bf q}_{o}$ and a given initial momentum ${\bf p}_{o}$ – these cases are
mutually exclusive from the Heisenberg uncertainty principle.

    1. 1.

Consider first the case on the left side of Figure 1 where ${\bf q}_{o}$ is
given at $t=0$ and where the particle passes through the slit at $t_{j}$.
Based on Theorem 1, the local extremized continous set-valued action $j=1,2$
from ${\bf q}_{o}$ to ${\bf q}$ is given by

${\phi}({\bf q},t,j)=\left\\{\begin{array}[]{ll}\frac{M}{2t}({\bf q}-{\bf
q}_{o})^{T}({\bf q}-{\bf q}_{o})&\text{ for \ }q_{1}<0\\\
\frac{M}{2(t-t_{j})}({\bf q}-{\bf q}_{j})^{T}({\bf q}-{\bf
q}_{j})+\frac{M}{2t_{j}}({\bf q}_{j}-{\bf q}_{o})^{T}({\bf q}_{j}-{\bf
q}_{o})&\text{ for \ }q_{1}\geq 0\end{array}\right.$

with one corner in the slit, marked as star in Figure 1. Note each ${\phi}$ is
continuous, but not Lipschitz continuous, at $t=0$ and $t=t_{j}$. The
trajectory can bend due to the constraint Dirac collision force in any
direction at the corner, which has normals in all directions.

This can also be seen from the second covariant derivative (19) and the
implied contraction rate (20) whose position independence allows to select the
time-scaling $T=1$

$\nabla_{M}^{2}\phi=\left\\{\begin{array}[]{ll}\frac{M}{t}{\bf I}&\text{ for \
}q_{1}<0\\\ \frac{M}{t-t_{j}}{\bf I}&\text{ for \ }q_{1}\geq
0\end{array}\right.\ \ \ \ \ \ \ \ \ \ \
\lambda_{\max}=\lambda_{\min}=\left\\{\begin{array}[]{ll}\frac{1}{t}&\text{
for \ }q_{1}<0\\\ \frac{1}{t-t_{j}}&\text{ for \ }q_{1}\geq
0\end{array}\right.$

Hence the geodesic distance $s=M\sqrt{({\bf q}_{2}-{\bf q}_{1})^{T}({\bf
q}_{2}-{\bf q}_{1})}$ between two trajectories ${\bf q}_{1}(t)$, ${\bf
q}_{2}(t)$ evolves as (20)

$s(t)=s(0)\ e^{\int\limits_{o}^{t}\lambda_{\max}({\bf
q},t)dt}=s(0)\left\\{\begin{array}[]{ll}c_{o}\ t&\text{ for \ }q_{1}<0\\\
c_{j}\ (t-t_{j})&\text{ for \ }q_{1}\geq 0\end{array}\right.$

with constant diffraction factors $c_{j}$ . At the initial point ${\bf q}_{o}$
and at the slit ${\bf q}_{j}\ $, the position measurement collapses the action
and leads to a radially divergent action or wave.

The set-valued quantum wave function can be computed from Theorem 2 and the
set-valued action

${\Psi}^{\prime}({\bf q},t,j)=e^{\frac{i}{\hbar}\phi({\bf
q},t,j)-\int\frac{\Delta_{M}\phi}{2}dt}=\left\\{\begin{array}[]{ll}\frac{M}{2\pi
i\hbar t}e^{\frac{i}{\hbar}\phi}&\text{for \ }q_{1}<0\\\ \frac{M}{2\pi
i\hbar(t-t_{j})}e^{\frac{i}{\hbar}\phi}&\text{for \ }q_{1}\geq
0\end{array}\right.$

which are Dirac impulses at $t=0$ and $t=t_{j}$. At $t_{j}$ the energy can
change due to the partial elastic collision force (8) and the wave function
collapses due to the slit. Finally, both solutions $j=1,2$ are equally
weighted and superposed according to (17).

    2. 2.

Consider now the case on the right side of Figure 1 where ${\bf
p}_{o}=(p_{o},0)$ is given at $t=0$. Based on Theorem 1 the set-valued action
is given for $j=1,2$ with constant energy $E=\frac{p_{o}^{2}}{2M}$ and
$Et_{1}={\bf p}_{o}^{T}{\bf q}_{1}={\bf p}_{o}^{T}{\bf q}_{2}$ by

${\phi}({\bf q},t,j)=\left\\{\begin{array}[]{ll}{\bf p}_{o}^{T}{\bf q}({\bf
q}^{\prime})-Et&\text{for \ }q_{1}<0\\\ p_{o}r_{j}-E(t-t_{1})&\text{for \
}q_{1}\geq 0\end{array}\right.$ (22)

where we use polar coordinates with metric

${\bf q}^{\prime}=\left(\begin{array}[]{c}r_{j}\\\
\theta_{j}\end{array}\right)=\left(\begin{array}[]{c}\sqrt{({\bf q}-{\bf
q}_{j})^{T}({\bf q}-{\bf q}_{j})},\\\
\arctan\frac{q_{1}-q_{1j}}{q_{2}-q_{2j}}\end{array}\right),\ \ {\bf
M}^{\prime}=\frac{\partial{\bf q}}{\partial{\bf
q}^{\prime}}^{T}\frac{\partial{\bf q}}{\partial{\bf
q}^{\prime}}=\left(\begin{array}[]{cc}1&0\\\ 0&r_{j}^{2}\end{array}\right)$

Note ${\phi}({\bf q}({\bf q}^{\prime}),t,j)$ is continuous even at the slit
but in the slit the momentum $\frac{\partial\phi}{\partial{\bf q}}({\bf
q},t,j)$ is discontinuous, i.e. it can change directions from ${\bf p}_{o}$ in
any direction $0\leq\theta_{j}\leq\pi$ with a fully elastic collision force in
Theorem 1.

This can also be seen from the second covariant derivative (19) of the action
and the implied contraction rate (20)

$\displaystyle\nabla_{M^{\prime}}^{2}\phi$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}{\bf 0}&\text{ for \ }q_{1}<0\\\
\left(\begin{array}[]{cc}0&0\\\ 0&\dot{r}_{j}r_{j}\end{array}\right)&\text{
for \ }q_{1}\geq 0\end{array}\right.$ (27) $\displaystyle\lambda_{\min r_{j}}$
$\displaystyle=$ $\displaystyle\lambda_{\max r_{j}}=0$
$\displaystyle\lambda_{\min\theta_{j}}$ $\displaystyle=$
$\displaystyle\lambda_{\max\theta_{j}}=\left\\{\begin{array}[]{ll}0&\text{ for
\ }q_{1}<0\\\ \frac{\dot{r}_{j}}{r_{j}}&\text{ for \ }q_{1}\geq
0\end{array}\right.$ (30)

where the $r_{j}$ and $\theta_{j}$ dynamics is decoupled, leading to two
independent contraction rates. The radial dynamics is indifferent. Howeve the
angular distance
$s=\min\int_{\theta_{j1}}^{\theta_{j2}}|\delta\theta_{j}r_{j}|$ between two
trajectories ${\bf q}_{1}(t)$, ${\bf q}_{2}(t)$ evolves with (20) as

$s(t)=s(0)\left\\{\begin{array}[]{ll}1&\text{ for \ }q_{1}<0\\\ c_{j}\
r_{j}&\text{ for \ }q_{1}\geq 0\end{array}\right.$

with constant diffraction factors $c_{j}$. Thus, behind the slit the lateral
distance scales with $r_{j}\ $, corresponding to diffraction with constant
momentum.

The quantum wave functions can now be computed with Theorem 2 from the set-
valued action (22), Laplacian and time-scaling (9)

$\displaystyle\Delta_{M^{\prime}}^{2}\phi$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&\text{ for \ }q_{1}<0\\\
\frac{\dot{r}_{j}}{r_{j}}&\text{ for \ }q_{1}\geq 0\end{array}\right.$ (33)
$\displaystyle T_{j}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}1&\text{ for \ }q_{1}<0\\\
\frac{\dot{r}_{j}}{r_{j}}&\text{ for \ }q_{1}\geq 0\end{array}\right.\ \ \ \ \
\tau=\left\\{\begin{array}[]{ll}\tau&\text{ for \ }q_{1}<0\\\
\ln{\frac{r_{j}}{r_{jo}}}&\text{ for \ }q_{1}\geq 0\end{array}\right.$ (38)

with constant $r_{o}$ as

${\Psi}^{\prime}({\bf q},t,j)=e^{\frac{i}{\hbar}\phi({\bf
q},t,j)-\int\frac{\Delta_{M^{\prime}}\phi}{2T_{j}}d\tau}=\left\\{\begin{array}[]{ll}e^{\frac{i}{\hbar}({\bf
p}_{o}^{T}{\bf q}-Et)}&\text{for \ }q_{1}<0\\\
\frac{r_{jo}}{r_{j}}e^{\frac{i}{\hbar}(p_{o}r_{j}-E(t-t_{1}))}&\text{for \
}q_{1}\geq 0\end{array}\right.$

Finally, both solutions $j=1,2$ are equally weighted and superposed according
to (17) which leads for $q_{1}\geq 0$ to the well known two-slit Frauenhofer
wave function

$\Psi^{\prime}({\bf q},t)={\psi}({\bf q},t,1)+{\psi}({\bf q},t,2)\ =\
e^{-\frac{i}{\hbar}E(t-t_{1})}r_{o}\left(\frac{e^{\frac{i}{\hbar}p_{o}r_{1}}}{r_{1}}+\frac{e^{\frac{i}{\hbar}p_{o}r_{2}}}{r_{2}}\right)$
(39)

While both results match well-known quantum results [8], the key difference is
that they are now fully consistent with the classical Hamilton-Jacobi dynamics
of Theorem 1 for ${\bf q}\in\mathbb{G}^{2}$. $\Box$

Figure 1: Two local least action solutions $\phi({\bf q},t,1),\phi({\bf
q},t,2)$ of the double slit experiment

  * ###### Example 4.2

:  Consider the single slit experiment in Figure 2 with Hamiltonian
$h=\frac{M}{2}\frac{\partial\phi}{\partial{\bf
q}}^{T}\frac{\partial\phi}{\partial{\bf q}}$ and constant mass $M$ for a slit
of finite width. We define the action on ${\bf
q}\in\mathbb{G}^{2}=\mathbb{R}^{2}\setminus\mathbb{E}^{2}$ rather than on
$\mathbb{R}^{2}$, thus excluding the red double slit wall obstacle
$\mathbb{E}^{2}$ in Figure 2.

According to Example 4.4 or [11] a particle corresponds to a wave function
decaying as $e^{-\alpha q^{2}}$ or $e^{-\alpha|q|}$ where the origin of $q$ is
at the center of the particle. This means that if e.g. a Gaussian electron
passes a corner at a geometric distance of $q$, it still collides with the
corner with state magnitude $e^{-\alpha q^{2}}$ or $e^{-\alpha|q|}$. Hence
particles which pass through a slit of finite size always collide with the
corner, but with a small state magnitude. According to Theorem 1, the least
actions $\phi({\bf q},t,j),j\in\mathbb{N}^{+}$ of initial momentum ${\bf
p}_{o}$ are given by (22), with corners at the slit marked as star in Figure
2. At the slit, the trajectory can bend in any direction due to the Dirac
elastic collision force constraint, which has normals in all directions.

Finally all solutions (22), $j\in\mathbb{N}^{+}$ are equally weighted and
superposed according to (17) which leads to the well known single-slit
Frauenhofer wave function (39)

$\Psi^{\prime}({\bf q},t)=\sum_{j\in\mathbb{N}^{+}}{\psi}({\bf
q},t,j)=e^{-\frac{i}{\hbar}E(t-t_{1})}r_{o}\sum_{j\in\mathbb{N}^{+}}\frac{e^{\frac{i}{\hbar}p_{o}r_{j}}}{r_{j}}$

Again, while this result matches the computation in e.g. [8], the key
difference is that it is fully consistent with classical physics using the
least action (1) of Theorem 1 for ${\bf q}\in\mathbb{G}^{2}$. $\Box$

Figure 2: $j\in\mathbb{N}^{+}$ local least action solutions $\phi({\bf
q},t,j)$ of the single slit experiment

  * ###### Example 4.3

:  Consider the particle in a box of Figure 3 with Hamiltonian
$h=\frac{1}{2M}\frac{\partial\phi}{\partial{q}}^{2}$ and constant mass $M$. We
define the action on ${\bf
q}\in\mathbb{G}^{1}=\mathbb{R}^{1}\setminus\mathbb{E}^{1}$ or $l\geq q\geq 0$
rather than on $\mathbb{R}^{1}$, thus excluding the upper and lower wall in
Figure 3.

The set-valued action solution (11) of Theorem 1 is illustrated in Figure 3
for $j\in\mathbb{Z}$ reflections, where $j$ is positive if the particle moton
at $t$ is up and $0$ or negative if the particle motion at $t$ is down,

${\phi}(q,t,j)=\left\\{\begin{array}[]{ll}p(jl+q)-Et&\text{for \
}j\in\mathbb{N^{+}}\\\ p(jl+l-q)-Et&\text{for \
}j\in\mathbb{Z}\setminus\mathbb{N^{+}}\\\ \end{array}\right.$

with constant momentum $p$ and energy $E=\frac{1}{2M}p^{2}$. The set-valued
wave function (15) of Theorem 2 is

${\psi}(q,t,j)=e^{\frac{i}{\hbar}{\phi}(q,t,j)-\int\frac{\Delta_{M}\phi}{2}dt}\
=\ \left\\{\begin{array}[]{ll}e^{\frac{i}{\hbar}(p(jl+q)-Et)}&\text{for \
}j\in\mathbb{N^{+}}\\\ e^{\frac{i}{\hbar}(p(jl+l-q)-Et)}&\text{for \
}j\in\mathbb{Z}\setminus\mathbb{N^{+}}\\\ \end{array}\right.$

Superposing according to (17) of Theorem 2 all reflections with the Euler
formula leads to

${\Psi}^{\prime}=\sum_{j\in\mathbb{Z}}\psi(q,t,j)=\sum_{j\in\mathbb{N^{+}}}\frac{1}{2}e^{\frac{i}{\hbar}(p\
j\ l-E\ t)}\sin{\frac{pq}{\hbar}}$

The symmetry
${\Psi}^{\prime}(\frac{l}{2}+q,t)=\pm{\Psi}^{\prime}(\frac{l}{2}-q,t)$ implies
the set-valued $\frac{p(j)}{\hbar}=\frac{\pi
j}{l},E_{j}=\frac{1}{2M}p^{2}(j)$. Normalising according to (17) finally leads
to

${\Psi}=\frac{\Psi^{\prime}}{\int_{\mathbb{G}^{n}}\sqrt{\Psi^{\prime\ast}\Psi^{\prime}}\
d{\bf
q}}=\sum_{j\in\mathbb{N^{+}}}e^{-\frac{i}{\hbar}E_{j}t}\sqrt{\frac{2}{l}}\sin\frac{\pi
jq}{l}$

This result matches the particle eigenfunction e.g. in [9, 11], with the key
difference that it is fully consistent with the constrained Hamilton-Jacobi
equation of Theorem 1 for ${q}\in\mathbb{G}^{1}$. $\Box$

Figure 3: $j\in\mathbb{Z}$ local least actions $\phi({\bf q},t,j)$ of the
particle in a box experiment

From this point of view, the key difference between classical least action
dynamics and quantum mechanics is the fact that classical action is defined in
unconstrained Riemann spaces ${\bf q}\in\mathbb{R}^{n}$ whereas the set-valued
least action of Theorem 1 is defined in a constrained Riemann space ${\bf
q}\in\mathbb{G}^{n}\subset\mathbb{R}^{n}$. Due to the non-Lipschitz Dirac
collision force of Theorem 1, these constraints imply multiple classical
action solutions which imply multiple wave functions. The superposition of
these wave functions leads to the final quantum state.

The following example illustrates that the approach is also applicable to a
single action and wave function solution, similar to Feynman’s harmonic
oscillator of Example 6. Note that our main result in Theorem 1 allows both
contexts to be combined.

  * ###### Example 4.4

:  Consider a particle in a $3$-dimensional Coulomb potential with time-scaled
Hamiltonian

$\frac{h}{T}\ -\ \frac{V}{T}-\
\frac{E}{T}=\frac{1}{2mT}\frac{\partial\phi}{\partial{\bf
q}}^{T}\frac{\partial\phi}{\partial{\bf
q}}=\frac{1}{2M}\frac{\partial\phi}{\partial{\bf
q}^{\prime}}^{T}\frac{\partial\phi}{\partial{\bf q}^{\prime}}\ \ \ \ \ \ \ \ \
\ \ T=\frac{V}{C}=\frac{1}{\sqrt{{\bf q}^{T}{\bf q}}}=\frac{1}{{\bf q}^{\prime
T}{\bf q}^{\prime}}$ (40)

with constant $C,M=4m,E=-\frac{M}{2}\omega^{2}$, $3$-dimensional Cartesian
position ${\bf q}({\bf q}^{\prime})$ centered at the origin of the Coulomb
potential and the $4$-dimensional ’square root coordinates’ ${\bf q}^{\prime}$
[4, 10] which fulfill Lemma 1

$\displaystyle q_{1}$ $\displaystyle=$ $\displaystyle
2q_{1}^{\prime}q_{3}^{\prime}+2q_{2}^{\prime}q_{4}^{\prime}$ $\displaystyle
q_{2}$ $\displaystyle=$
$\displaystyle-2q_{1}^{\prime}q_{2}^{\prime}+2q_{3}^{\prime}q_{4}^{\prime}$
$\displaystyle q_{3}$ $\displaystyle=$ $\displaystyle-q_{1}^{\prime
2}+q_{2}^{\prime 2}+q_{3}^{\prime 2}-q_{4}^{\prime 2}$

The action of the time-scaled Hamilton-Jacobi dynamics in Theorem 1 of the
Hamiltonian (40) is then given by the analytic action (14) in Lemma 1. Indeed
the Hamiltonian in ${\bf q}^{\prime}$ or ${\bf q}$ have equivalent
trajectories. Hence in the following we assume from the start the dynamics in
${\bf q}^{\prime}$ with time $\tau$.

The resulting normalized wave function (16, 17) of Theorem 2 is then

$\Psi({\bf q}^{\prime},\tau)=e^{\frac{i}{\hbar}\phi({\bf
q},\tau)-\int\frac{\Delta_{M}\phi(\tau)}{2}d\tau+\frac{iE}{\hbar}\tau}=\sum_{j_{1},...j_{4}=0}^{\infty}e^{-\frac{i}{\hbar}(E_{j}^{\prime}-E-C)\tau}\psi_{j_{1},...j_{4}}({\bf
q}^{\prime})\psi_{j_{1},...j_{4}}({\bf q}^{\prime}_{o})d\tau$

with the eigenvalues $E_{j}^{\prime}=\hbar\omega(j_{1}+...+j_{4}+2)$ and
eigenfunctions (6) from the harmonic oscillator Example 6. Note that
$\phi({\bf q}^{\prime},{\tau})=\phi(-{\bf q}^{\prime},{\tau})$ and hence
$\Psi({\bf q}^{\prime},E)=\Psi(-{\bf q}^{\prime},E)$ are point symmetric
functions. This point symmetry excludes odd polynomials in
$\psi_{j_{1},...j_{4}}({\bf q}^{\prime})$ such that we can write
$j_{1}+...+j_{4}=2j-2$ with the quantum number $j\in\mathbb{N}^{+}$. Hence we
can rewrite the above as

$\Psi({\bf q}^{\prime},\tau)=\sum_{\begin{subarray}{c}j\in\mathbb{N}^{+}\\\
j_{1}+...+j_{4}=2j-2\end{subarray}}e^{\frac{i}{\hbar}(C+E-2\hbar
j\sqrt{-\frac{2E}{M}})\tau}\psi_{j_{1},...j_{4}}({\bf
q}^{\prime})\psi_{j_{1},...j_{4}}({\bf q}^{\prime}_{o})d\tau$

whose Fourier transformation w.r.t. energy $E$ is with
$E_{j}=-\frac{m}{2}\frac{C^{2}}{\hbar^{2}j^{2}}$

$\displaystyle\mathscr{F}(\Psi)({\bf q}^{\prime},E)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}e^{-\frac{iE}{\hbar}\tau}\Psi({\bf
q}^{\prime},\tau)d\tau$ $\displaystyle=$ $\displaystyle
2\pi\sum_{\begin{subarray}{c}j\in\mathbb{N}^{+}\\\
j_{1}+...+j_{4}=2j-2\end{subarray}}\delta\frac{E-E_{j}}{\hbar}\
\psi_{j_{1},...j_{4}}({\bf q}^{\prime})\ \psi_{j_{1},...j_{4}}({\bf
q}^{\prime}_{o})$

where we used $\int_{-\infty}^{\infty}e^{i\beta\tau}d\tau=2\pi\delta(\beta)$,
$\forall\beta$. An inverse Fourier transformation then turns the dependence
from $E$ back into the actual time $t$

$\Psi({\bf q}^{\prime},t)=\mathscr{F}^{-1}(\Psi)({\bf
q}^{\prime},E)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{\frac{iE}{\hbar}t}\mathscr{F}(\Psi)dE=\sum_{\begin{subarray}{c}j\in\mathbb{N}^{+}\\\
j_{1}+...+j_{4}=2j-2\end{subarray}}e^{\frac{i}{\hbar}E_{j}t}\
\psi_{j_{1},...j_{4}}({\bf q}^{\prime})\ \psi_{j_{1},...j_{4}}({\bf
q}^{\prime}_{o})$

This result matches the 3-dimensional Coulomb result in [4], with the key
difference that it is based only on the classical Hamiltonian trajectories of
Theorem 1. Neither zig-zag trajectories nor time-slicing of the Feynman path
integral [4] are necessary, which considerably simplifies the computation.
Note that [4] algebraically convert the Cartesian Hermite polynomials above to
the more commonly used spherical Laguerre polynomials. $\Box$

All our examples avoid zig-zag trajectories or time-slicing of the Feynman
path integral [8, 4], which considerably simplifies the computations. This is
not a statement on the nature of physical reality, but simply on the fact that
the same actual results or measurements can be computed just from classical
actions.

##  5 Concluding Remarks

Theorem 1 extends the classical Lagrangian dynamics within an unconstrained
space $\mathbb{R}^{n}$ to a constrained space
$\mathbb{G}^{n}\subset\mathbb{R}^{n}$. At the border $\partial\mathbb{G}^{n}$
of $\mathbb{G}^{n}$, a Dirac Lagrangian constraint force ensures that the
constraint is not violated. This non-Lipschitz activation of the constraint
force can lead to multiple, but deterministic, constrained least action
solutions.

This result allows a classical least action interpretation of standard
experiments of quantum physics. In the double-slit experiment, for instance, a
main difference between quantum mechanics and the classical least action is
the fact that the former is defined in unconstrained Riemann spaces
$\mathbb{R}^{n}$, leading to a single deterministic solution of the least
action problem. By contrast, the non-Lipschitz Dirac collision force,
introduced by inequality constraints $\mathbb{G}^{n}\subset\mathbb{R}^{n}$,
leads to a set of multiple solutions of a deterministic least action
principle, which is intepreted as diffraction in quantum physics. To freely
paraphrase Einstein, God may not play dice but perhaps the world is just
ambiguous.

The approach of this paper also features a natural equivalence of the
microscopic linear Schrödinger equation (4) with the macroscopic Hamilton
equation (11), and thus a smooth transition across scales. In principle, this
may also extend to general relativity, which is also based on Hamiltonian
dynamics, with ${\bf M}$ a non-definite metric and $t$ replaced by proper
time. The equivalence was already known in [8] for the specific case of a
quadratic Lagrangian in $\mathbb{R}^{n}$, and equations (16), (18) generalize
it to non-quadratic Lagrangians and to constrained actions. For high-
dimensional problems as e.g. in computational quantum chemistry, the piecewise
differentiability of the classical action should also make numerical or
machine learning techniques more readily applicable.

## References

  * [1] Aspect A., Bell’s inequality test more ideal than ever, Nature 398(189), 1999
  * [2] Bell J.S., On the Einstein-Poldolsky-Rosen paradox, Physics 1(195), 1964.
  * [3] Cohen-Tannoudji C., Diu B., and Laloe F., Quantum Mechanics, 2nd Ed., Wiley, 2019.
  * [4] Duru I.H., Kleinert H., Quantum Mechanics of H-atoms from path integrals, Fortschritte der Physik, 1982.
  * [5] Einstein A., Podolsky B., and Rosen N., Can the quantum-mechanical description of physical reality be considered complete? Phys. Rev. 47 (777), 1935.
  * [6] Evans L.C., Partial Differential Equations, American Mathematical Society, 1991.
  * [7] Feynman R.P., Space-Time Approach to Non-Relativistic Quantum Mechanics, Review of Modern Physics, 20(2), 1948.
  * [8] Feynman R.P. and Hibbs A.R., Quantum Mechanics and Path Integrals, McGraw- Hill, New York, 1965
  * [9] Janke W., Kleinert H., Summing Paths for a Particle in the Box, Lettera al Luovo Cimento, 1979.
  * [10] Kleinert, H., Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, World Scientific, 2009.
  * [11] Liboff R.L., Introductory Quantum Mechanics, 4th Ed., Addison Wesley, 2002.
  * [12] Lohmiller W. and Slotine J.J.E., On Contraction Analysis for Nonlinear Systems, Automatica, 34(6), 1998.
  * [13] Lohmiller W. and Slotine, J.J.E., Contraction Theory with Inequality Constraints, arXiv:2306.06628, 2023.
  * [14] Lohmiller W. and Slotine, J.J.E., Contraction Analysis of Nonlinear Distributed Systems, International Journal of Control, 2005.
  * [15] Lovelock D. and Rund H., Tensors, Differential Forms, and Variational Principles, Dover, 1989.
  * [16] Magnus K. and Mueller H.H., Grundlagen der technischen Mechanik, Teubner, 1990.
  * [17] https://en.wikipedia.org/wiki/set-valued_function
  * [18] https://en.wikipedia.org/wiki/Nash_embedding_theorems
  * [19] Pelster, A. and Wunderlin, A, On the generalization of the Duru-Kleinert-propagator transformations, Zeitschrift für Physik B Condensed Matter, 89, 1992.
  * [20] Sakoda, S., On the effective potential of Duru-Kleinert path integrals, J. Math. Phys. 58, 062111, 2017.
  * [21] Tolman R.C., Relativity, Thermodynamics, and Cosmology, Dover Publications, 1934.


# Tracking moving objects through scattering media via speckle correlations

Y. Jauregui-Sánchez, Harry Penketh, and Jacopo Bertolotti

(Physics and Astronomy Department, University of Exeter, Stocker Road, Exeter
EX4 4QL, UK
Corresponding author<EMAIL_ADDRESS>
Scattering can rapidly degrade our ability to form an optical image, to the
point where only speckle-like patterns can be measured. Truly non-invasive
imaging through a strongly scattering obstacle is difficult, and usually
reliant on a computationally intensive numerical reconstruction. In this work
we show that, by combining the cross-correlations of the measured speckle
pattern at different times, it is possible to track a moving object with
minimal computational effort and over a large field of view.

Most optical imaging systems, including our eyes, use a lens system to
reproduce on the detector an appropriately modified (e.g. magnified or
filtered) version of the intensity profile on the object plane. This approach
requires light to propagate in straight lines, and as soon as the medium
between the object and the detector is inhomogeneous, the resulting image is
distorted or blurred [1]. If the scattering is weak enough, there is a
significant fraction of the light that is left unscattered and that can be
used to form a sharp image [2, 3, 4], but for strongly scattering media it
becomes impossible to filter out the unwanted signal [5]. Knowledge of the
exact distortion, or strong priors on what the image should look like, allow
for the detected image to be corrected and the desired information restored
[6, 7, 8]. On the other hand, reconstructing an unknown image through a
strongly scattering medium is a difficult and still largely open problem [9,
10].

A possible approach to this problem is based on the idea that there are
universal correlations in the scattered light, that do not depend on the fine
details of the medium. In particular, the optical memory effect [11, 12] tells
us that the speckle pattern we get when coherent light passes through a slab
of scattering material [13, 14] is strongly correlated with the speckle
pattern we get if we illuminate the scattering slab from a different angle
[15, 16, 17, 18]. This is a powerful tool for imaging, as it gives us
information on the light on the hidden side of the scattering medium without
ever having to measure anything there. In particular it allows us to measure
the autocorrelation of the hidden object [19, 20, 21, 22, 23, 24, 25]. The
autocorrelation can be inverted to obtain the shape of the object using a
Gerchberg-Saxton algorithm [26, 27, 28], but this is notoriously a
computationally intensive process. If the objects in the scene are moving,
there is no time to perform the inversion before the scene has changed.

In this work we show that, while the imaging problem is hard, tracking the
motion of objects hidden behind a scattering medium is much less so. In
particular, the proposed technique has a negligible computational cost, and
thus can be performed in real-time, and tracking can be performed well beyond
the memory effect range that otherwise limits the field of view [29].

Figure 1: Schematic representation of the experimental set-up. (a) An object
hidden at a distance $L_{o}$ behind a scattering medium emits light (e.g.
because it is illuminated). The scattered light forms a speckle pattern, which
is detected by a camera at a distance $L_{s}$ from the medium (b). While the
speckle itself looks random, its autocorrelation shows information about the
hidden object (c). In the present experiment we used $L_{o}=430$ mm and
$L_{s}=50$ mm.

If we consider an object either emitting or reflecting light hidden behind a
strongly scattering layer, each point $\mathbf{x}_{\text{o}}$ will radiate an
amount $O(\mathbf{x}_{\text{o}})$ of light towards the scattering screen (in
the following we will assume isotropic emission for simplicity). For a narrow-
band enough detection, the light coming from each point will create a speckle
pattern $S(\mathbf{x}_{\text{o}},\mathbf{x}_{\text{d}})$ at position
$\mathbf{x}_{\text{d}}$ on the detector (see Fig. 1). If the spatial coherence
is low enough, these speckle patterns will not interfere, but simply sum with
each other [30], and the total measured intensity on the detector will be

$I(\mathbf{x}_{\text{d}})=\int
O(\mathbf{x}_{\text{o}})\,S(\mathbf{x}_{\text{o}},\mathbf{x}_{\text{d}})\text{d}^{2}\mathbf{x}_{\text{o}}.$
(1)

In this form there is no obvious way to retrieve $O$, but if we autocorrelate
the measured intensity we obtain

$\displaystyle[I\star I](\boldsymbol{\Delta x}_{\text{d}})=\int
I(\mathbf{x}_{\text{d}})\,I(\mathbf{x}_{\text{d}}+\boldsymbol{\Delta
x}_{\text{d}})\text{d}^{2}\mathbf{x}_{\text{d}}=$ (2) $\displaystyle=$
$\displaystyle\int\left[\left(\int
O(\mathbf{x}_{\text{o}})\,S(\mathbf{x}_{\text{o}},\mathbf{x}_{\text{d}})\text{d}^{2}\mathbf{x}_{\text{o}}\right)\right.$
$\displaystyle\left.\left(\int
O(\mathbf{y}_{\text{o}})\,S(\mathbf{y}_{\text{o}},\mathbf{x}_{\text{d}}+\boldsymbol{\Delta
x}_{\text{d}})\text{d}^{2}\mathbf{y}_{\text{o}}\right)\right]\,\text{d}^{2}\mathbf{x}_{\text{d}}=$
$\displaystyle=$ $\displaystyle\int
O(\mathbf{x}_{\text{o}})\,O(\mathbf{y}_{\text{o}})\Big{(}\left[S\star
S\right](\mathbf{x}_{\text{o}},\mathbf{y}_{\text{o}},\boldsymbol{\Delta
x}_{\text{d}})\Big{)}\text{d}^{2}\mathbf{x}_{\text{o}}\text{d}^{2}\mathbf{y}_{\text{o}},$

where $\star$ represents the cross-correlation product. Making an ensemble
average $\left\langle\cdot\right\rangle$ over disorder, we can rewrite this in
a compact form in terms of the speckle correlations
$\mathcal{C}=\frac{\left\langle\delta S\star\delta
S\right\rangle}{\overline{S}^{2}}$ [15]:

$\displaystyle\,\left\langle I\star I\right\rangle(\boldsymbol{\Delta
x}_{\text{d}})$ $\displaystyle\propto\left[O\star O\right]\otimes\mathcal{C}$
(3)

where $\delta S=S-\overline{S}$, $\overline{S}$ is the average intensity of
the speckle on the object plane, and $\otimes$ is the convolution product. See
the Supplementary Information: S1 for the step by step calculation, and a
discussion of the meaning of each term.

Equation (3) holds for any correlation between the speckle generated by
incoherent sources at a distance $\boldsymbol{\Delta
x}_{\text{o}}=\mathbf{x}_{\text{o}}-\mathbf{y}_{\text{o}}$, and in particular
it works for the optical memory effect [12]:

$\mathcal{C}\simeq e^{-k^{2}\left(\boldsymbol{\Delta
x}_{\text{o}}-\boldsymbol{\Delta
x}_{\text{d}}\right)^{2}\sigma^{2}}\cdot\left(\frac{k\lvert\boldsymbol{\Delta\theta}\rvert
L}{\sinh\left(k\lvert\boldsymbol{\Delta\theta}\rvert L\right)}\right)^{2},$
(4)

where we assumed that the illumination is a Gaussian beam with variance
$\sigma^{2}$, $\boldsymbol{\Delta\theta}\simeq\boldsymbol{\Delta
x}_{\text{o}}/L_{o}$ (with $L_{o}$ the distance between the object plane and
the scattering medium, see Fig. 1a), and $L$ the scattering medium thickness.
As per equation (3), $\mathcal{C}$ behaves like the point spread function of
our measurement of $O\star O$. Its first term tells us that the
autocorrelation is sharply peaked around $\boldsymbol{\Delta
x}_{\text{o}}=\boldsymbol{\Delta x}_{\text{d}}$, and the second that the
correlation decreases rapidly for large $\boldsymbol{\Delta x}_{\text{o}}$,
resulting in a finite range where this correlation can be exploited (the
isoplanatic patch) [9]. Within this range, we can therefore measure $O\star O$
with a resolution that depends on the wavelength and the width of the
illumination beam, but not on the scattering properties of the medium. And
while it is not possible to analytically invert an autocorrelation, it is
possible to do so numerically under very mild assumptions [31, 32, 33].

If the scene one wants to image is not static, it is possible to repeat the
process above for each desired frame, but since the autocorrelation inversion
is a computationally intensive process, prone to get stuck in local minima
[28], this can quickly become unwieldy. That said, in many cases one might be
satisfied in knowing how the scene moved at any given moment, instead of
imaging every single frame. In this case it is possible to significantly
simplify the problem.

If the whole scene moved rigidly, it is possible to reconstruct accurately by
how much by performing the correlation of the measured light at time $t_{1}$
with the measured light at time $t_{0}$. This results in the very same
autocorrelation $O\star O$, just shifted by exactly the amount the scene moved
[34, 35]. But going beyond rigid movements requires some care. Let’s assume
that the object we want to measure is made of a number of smaller components,
each moving independently from the others: $O=\sum_{j}o_{j}$. If we try to
correlate the intensity at time $t_{1}$ with the intensity at time $t_{0}$ we
obtain

$\displaystyle I(t_{1})\star
I(t_{0})\propto\left[\left(\sum_{j}o_{j}(t_{1})\right)\star\left(\sum_{j}o_{j}(t_{0})\right)\right]\otimes\mathcal{C}=$
(5) $\displaystyle=\left[\sum_{j}o_{j}(t_{0})\star o_{j}(t_{1})+\sum_{j\neq
i}o_{j}(t_{0})\star o_{i}(t_{1})\right]\otimes\mathcal{C}.$

The problem with equation (5) is that, while the first term tells us how much
a given sub-object moved, the second term contains all cross correlations
between all the sub-objects comprising the scene, including information about
all the sub-objects that never moved. If we assume for simplicity that only a
single sub-object $o_{k}(t)$ is moving, we can split the above equation into a
static and a moving part

$\displaystyle I(t_{1})\star I(t_{0})\propto\left[\sum_{i,j\neq k}o_{i}\star
o_{j}+\sum_{j}o_{k}(t)\star o_{j}\right]\otimes\mathcal{C},$ (6)

i.e. a part of the cross-correlation will stay static, and another will move
following the moving object. If the number of overall sub-objects $o_{j}$ is
small, the moving part of the autocorrelation is easy to see and interpret
(see Fig. 2 and Supplementary Video 1), but if there are many sub-objects the
static ones will dominate the cross-correlation, making the moving part
difficult to spot. Furthermore, the moving part now contains the cross-
correlations with the many non-moving objects, which makes it even harder to
parse the image (see Fig. 3b, and Supplementary Video 2). A solution to this
problem is to calculate $I(t_{1})\star I(t_{0})-I(t_{1})\star I(t_{1})$
instead. If only the $k$th sub-object moves, this reads as

$I(t_{1})\star I(t_{0})-I(t_{1})\star
I(t_{1})\propto\left[\sum_{j}\left[o_{k}(t_{0})-o_{k}(t_{1})\right]\star
o_{j}\right]\otimes\mathcal{C}.$ (7)

If there are more sub-objects moving, each will yield a term similar to this
one. This subtraction has several advantages: first of all, it removes most of
the information on the objects that are not moving from the picture, leaving
only information on their relative distance from the moving sub-objects.
Furthermore, even if two sub-objects are too far from each other to be
contained within the memory range, we still have information about how the
sub-object is moving with respect to the sub-objects near to it. This allows
us to track the moving sub-object(s) beyond the memory range (see Fig. 3 and
Supplementary Video 2). It is worth noting that one can perform the cross-
correlation between the most recent frame and any other previous frames, not
necessarily the previous one. This degree of freedom can be used to maximize
the visibility of the moving object, depending on the speed of the moving
object(s) and the acquisition frame rate. In fact, if the object did not move
enough between the frames we correlate, $o_{k}(t_{0})-o_{k}(t_{1})$ in
equation 7 will be close to zero, producing a weak signal (see Supplementary
Video 3).

Figure 2: Experimental results of tracking moving objects via speckle
correlations in the case of few objects in the scene. (a) Displayed intensity
patterns on the DMD at different times. The moving object is highlighted in
red in the drawing. (b) Autocorrelation of the speckle pattern measured at
different times, and (c) cross-correlation between the measured speckle at
time $t_{0}$ and that at time $t_{n}$ (see Supplementary Video 1).

As an experimental validation of the technique described above, we generated a
moving scene on a digital micromirror device (DMD, Texas Instruments DLP9000),
and illuminated it with the light from a Red He-Ne laser beam (HNLS008R,
Thorlabs), passed through a rotating diffuser (DG20-1500 Ground Glass
Diffuser, Thorlabs) to reduce its spatial coherence (see Fig. 1a). Between the
DMD and the digital camera (Allied Vision Manta G-125B, Edmund Optics) we
place an optical diffuser (DG10-220-MD Ground Glass Diffuser, Thorlabs) to
scatter all the light and produce a speckle pattern (see Fig. 1b). The optical
diffuser was changed for two layers of adhesive sellotape in the second part
of the experiment to limit the memory effect range. It is important to notice
that the spatially incoherent illumination is not collimated, so the light
reflected by the DMD opens in a wide cone, and no image of the scene is formed
on the scattering medium.

If the scene is composed of few objects like in Fig. 2a, the time evolution of
the simple autocorrelation of the measured speckle can tell us a lot about how
it changed. For instance in Fig. 2b (and Supplementary Video 1) we can see
that something did not move, as the top and bottom part of the autocorrelation
do not change, while something else both got further away and rotated. The
problem with looking at the time evolution of the autocorrelation is that,
even with such a simple case, it is difficult to extract much more
information, unless one has very strong priors on the underlying scene. A
significant improvement can be obtained by looking at the cross-correlation of
the most recent measured speckle with a previously measured one (in this case,
the one measured at $t=t_{0}$). As shown in Fig. 2c, the movement of a single
object is visible as a part of the autocorrelation detaching. As per equation
6, the moving part is composed of the cross-correlation of the moving object
with all the other objects (including itself), which in this simple case
immediately tells us that the underlying scene is composed of 3 objects, of
which only one is moving, and it is both translating and rotating. Notice that
the amount we see it moving in the cross-correlation is linearly proportional
to real lateral movement, with a scaling factor $L_{s}/L_{o}$. Movements in
the direction of the scattering medium are not directly visible in the cross-
correlation, as this method has an infinite depth of field, but a closer
object will appear larger, giving a clear visual cue. It is important to
notice that, due to the geometry of the imaging system, the axes of the
speckle image, and thus of all autocorrelations and cross-correlations, are
reversed with respect of the original hidden scene (similarly to a pinhole
camera). This is easy to correct by flipping the images if deemed
inconvenient. It is also important to notice that, as the autocorrelation of a
real function is always centrosymmetric, it is only ever possible to detect
relative motion with respect with the position of the objects at the reference
time $t_{0}$, and not absolute movement.

Figure 3: (a) Three representative still frames of the moving scene on the
DMD, made of a large number of static dots, and a single moving star. (b) The
cross-correlation between the measured speckle patterns at $t_{0}$ and the
speckle patterns at successive times. The cross-correlation decreases with
distance from the centre due to the finite memory effect range. (c) Plot of
$I(t_{n})\star I(t_{n-2})-I(t_{n})\star I(t_{n})$, where $I(t_{n})$ is the
measured speckle pattern at frame $n$. The red arrow follows the cross-
correlations between one of the static background dots and the moving star. In
the first panel the background dot is at the edge of the memory effect range
from the moving star, and it is thus barely visible, becoming brighter while
it gets closer to the centre of the picture (i.e. closer to the moving star).

If the scene becomes more complex, with many objects present (e.g. the object
shown if Fig. 3a), the simple cross-correlation becomes difficult to interpret
(see Fig. 3b, and Supplementary Video 2). In the case of a single (or a few)
objects moving over an otherwise static background, the static part will
dominate the cross-correlations, thus making it impossible to see the
movement. If we use equation 7 and look at the difference between the cross-
correlation and the autocorrelation, we effectively subtract all the static
background, leaving just the informative part of the picture, i.e. the cross-
correlation of the difference between the moving object at time $t_{0}$ and
$t_{1}$ with the static background. As $o_{k}(t_{0})-o_{k}(t_{1})$ is negative
in the area now occupied by the moving object, and positive in the area it
previously occupied, the resulting image (see Fig. 3c) has zero mean, and the
negative part show the direction of movement of the object. It is important to
notice that the moving object (or objects, if there is more than one) is
always in the centre of the picture, with the (static) background moving
around it, essentially showing the movement from the frame of reference of the
moving object. Once again, this can be easily corrected numerically if deemed
inconvenient. For this experiment we changed the optical diffuser for two
layers of sellotape to limit the memory effect range, as visible in Fig. 3b,
and show that it is possible to track a moving object well beyond the memory
range. As shown in Fig. 3c (and Supplementary Video 2), only the cross-
correlations with background objects closer to the moving one that the memory
effect range are visible (see equation 4). Therefore, as long as the moving
object is close enough to at least one static background object, it is
possible to track its movements, effectively mosaicing the memory range (see
Supplementary Video 3 and 4).

In this work we show that, while imaging through a strongly scattering medium
is still (despite the large amount of work done in the field) a largely
unsolved problem, tracking the movement of an object admits a much simpler
solution. In particular, an appropriate linear combination of the cross-
correlations of the speckle-like images at different times, produces a human-
interpretable sequence of frames, which allows the tracking in real time and
with minimal computation. The motion to be tracked is by no means limited to
simple lateral translations, but includes rotations, changes of size,
deformations etc. The main limiting factor of this approach (like all non-
invasive imaging approaches) is the small amount of light that is going to
pass through the scattering medium in any realistic situation. Sensitive and
fast detectors will be required for real-world applications [36].

Methods Summary A He-Ne laser beam was expanded and passed through a spinning
diffuser to reduce its spatial coherence. The resulting divergent beam was
used to illuminate a moving scene generated on a DMD (or a reflecting sample,
see Supplementary Information: S3). The reflected light passed through a
220-grit ground glass diffuser (DG10-220-MD, Thorlabs), and the resulting
speckle-like intensity pattern was collected lenslessly by a CMOS camera with
an exposure time of 3 seconds. To obtain the results in Fig. 3, the scattering
layer was changed into two layers of adhesive sellotape, in order to reduce
the memory effect range, while keeping a comparable transmission. See
Supplementary Information: S2 for further details.

Supplementary Information is available in the online version of the paper.

Data accessibility Raw data for all the figures, and the MATLAB code used to
process them, is available at https://doi.org/10.5281/zenodo.6124320.

Acknowledgements The research was financially supported by the Engineering and
Physical Sciences Research Council (EPSRC) and Dstl via the University Defence
Research Collaboration in Signal Processing (EP/S000631/1).

Author Information JB developed the original idea, and performed the
calculations and numerical simulations. YJS designed and set up the
experimental apparatus. YJS and HP performed the experiments and the data
analysis. All authors wrote the manuscript. The authors declare no competing
financial interests.

## References

  * [1] S. Rotter and S. Gigan, “Light fields in complex media: Mesoscopic scattering meets wave control,” Rev. Mod. Phys., vol. 89, p. 015005, 2017.
  * [2] D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science, vol. 254, p. 1178, 1991.
  * [3] J. Pawley, ed., Handbook Of Biological Confocal Microscopy.  Springer, 2006.
  * [4] W. R. Zipfel, R. M. Williams, and W. W. Webb, “Nonlinear magic: multiphoton microscopy in the biosciences,” Nature biotechnology, vol. 21, no. 11, pp. 1369–1377, 2003.
  * [5] A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Smart optical coherence tomography for ultra-deep imaging through highly scattering media,” Science Advances, vol. 2, p. e1600370, 2016.
  * [6] R. Davies and M. Kasper, “Adaptive optics for astronomy,” Annual Review of Astronomy and Astrophysics, vol. 50, p. 305, 2012.
  * [7] I. M. Vellekoop and A. P. Mosk, “Universal optimal transmission of light through disordered materials,” Phys. Rev. Lett., vol. 101, no. 12, p. 120601, 2008.
  * [8] S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nat. Commun., vol. 1, p. 81, 2010\.
  * [9] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nat. Photonics, vol. 6, p. 283, 2012.
  * [10] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, “Deep optical imaging within complex scattering media,” Nature Reviews Physics, vol. 2, p. 141, 2020.
  * [11] S. Feng, C. Kane, P. A. Lee, and A. D. Stone, “Correlations and fluctuations of coherent wave transmission through disordered media,” Phys. Rev. Lett., vol. 61, pp. 834–837, 1988.
  * [12] I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Phys. Rev. Lett., vol. 61, pp. 2328–2331, 1988.
  * [13] J. W. Goodman, “Some fundamental properties of speckle,” J. Opt. Soc. Am., vol. 66, pp. 1145–1150, 1976.
  * [14] J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications.  SPIE, 2020.
  * [15] E. Akkermans and G. Montambaux, Mesoscopic Physics of Electrons and Photons.  Cambridge University Press, 2007.
  * [16] B. Judkewitz, R. Horstmeyer, I. Vellekoop, P. I.N., and C. Yang, “Translation correlations in anisotropically scattering media,” Nature Physics, vol. 11, p. 684, 2015.
  * [17] G. Osnabrugge, R. Horstmeyer, I. Papadopoulos, B. Judkewitz, and I. Vellekoop, “Generalized optical memory effect,” Optica, vol. 4, p. 886, 2017.
  * [18] S. Schott, J. Bertolotti, J.-F. Léger, L. Bourdieu, and S. Gigan, “Characterization of the angular memory effect of scattered light in biological tissues,” Optics Express, vol. 23, no. 10, p. 13505, 2015.
  * [19] L. I. Goldfischer, “Autocorrelation function and power spectral density of laser-produced speckle patterns,” J. Opt. Soc. Am., vol. 55, pp. 247–253, 1965.
  * [20] A. Labeyrie, “Attainment of diffraction limited resolution in large tele- scopes by fourier analysing speckle patterns in star images,” Astron. & Astrophys., vol. 6, pp. 85–87, 1970.
  * [21] J. Bertolotti, E. G. Van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non-invasive imaging through opaque scattering layers,” Nature, vol. 491, no. 7423, pp. 232–234, 2012.
  * [22] O. Katz, E. Small, and Y. Silberberg, “Looking around corners and through thin turbid layers in real time with scattered incoherent light,” Nat. Photonics, vol. 6, no. 8, pp. 549–553, 2012.
  * [23] M. Hofer, C. Soeller, S. Brasselet, and J. Bertolotti, “Wide field fluorescence epi-microscopy behind a scattering medium enabled by speckle correlations,” Opt. Express, vol. 26, no. 8, pp. 9866–9881, 2018.
  * [24] G. Stern and O. Katz, “Noninvasive focusing through scattering layers using speckle correlations,” Opt. Lett., vol. 44, pp. 143–146, Jan 2019.
  * [25] D. Wang, S. K. Sahoo, X. Zhu, G. Adamo, and C. Dang, “Non-invasive super-resolution imaging through dynamic scattering media,” Nat. Commun., vol. 12, no. 1, pp. 1–9, 2021.
  * [26] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik, vol. 35, pp. 1–6, 1972.
  * [27] J. R. Fienup, “Reconstruction of an object from the modulus of its fourier transform,” Opt. Lett., vol. 3, no. 1, pp. 27–29, 1978.
  * [28] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt., vol. 21, no. 15, pp. 2758–2769, 1982.
  * [29] A. Boniface, J. Dong, and S. Gigan, “Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix,” Nat Commun, vol. 11, p. 6154, 2020.
  * [30] J. W. Goodman, Statistical Optics.  John Wiley & Sons, 1985.
  * [31] J. C. Dainty, Laser Speckle and Related Phenomena.  Springer, 1984.
  * [32] O. Katz, P. Heidmann, M. Fink, and S. Gigan, “Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations,” Nat. Photonics, vol. 8, pp. 784–790, 2014.
  * [33] H. Yilmaz, E. G. van Putten, J. Bertolotti, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Speckle correlation resolution enhancement of wide-field fluorescence imaging,” Optica, vol. 2, pp. 424–429, 2015.
  * [34] M. Cua, E. H. Zhou, and C. Yang, “Imaging moving targets through scattering media,” Opt. Express, vol. 25, no. 4, pp. 3935–3945, 2017.
  * [35] C. Guo, J. Liu, T. Wu, L. Zhu, and X. Shao, “Tracking moving targets behind a scattering medium via speckle correlation,” Appl. Opt., vol. 57, no. 4, pp. 905–913, 2018.
  * [36] S. Gigan, O. Katz, H. B. de Aguiar, E. R. Andresen, A. Aubry, J. Bertolotti, E. Bossy, D. Bouchet, J. Brake, S. Brasselet, Y. Bromberg, H. Cao, T. Chaigne, Z. Cheng, W. Choi, T. Čižmár, M. Cui, V. R. Curtis, H. Defienne, M. Hofer, R. Horisaki, R. Horstmeyer, N. Ji, A. K. LaViolette, J. Mertz, C. Moser, A. P. Mosk, N. C. Pégard, R. Piestun, S. Popoff, D. B. Phillips, D. Psaltis, B. Rahmani, H. Rigneault, S. Rotter, L. Tian, I. M. Vellekoop, L. Waller, L. Wang, T. Weber, S. Xiao, C. Xu, A. Yamilov, C. Yang, and H. Yılmaz, “Roadmap on wavefront shaping and deep imaging in complex media,” arXiv, p. 2111.14908, 2021.