# Fast reconstruction of programmable integrated interferometers

Boris Bantysh  1,2,* Konstantin Katamadze  1 Andrey Chernyavskiy  1,2 and
Yurii Bogdanov1,2  1Valiev Institute of Physics and Technology, Russian
Academy of Sciences, 117218, Moscow, Russia,
2Russian Quantum Center, Skolkovo, Moscow 143025, Russia
<EMAIL_ADDRESS>
###### Abstract

Programmable linear optical interferometers are important for classical and
quantum information technologies, as well as for building hardware-accelerated
artificial neural networks. Recent results showed the possibility of
constructing optical interferometers that could implement arbitrary
transformations of input fields even in the case of high manufacturing errors.
The building of detailed models of such devices drastically increases the
efficiency of their practical use. The integral design of interferometers
complicates its reconstruction since the internal elements are hard to
address. This problem can be approached by using optimization algorithms [Opt.
Express 29, 38429 (2021)]. In this paper, we present a novel efficient
algorithm based on linear algebra only, which does not use computationally
expensive optimization procedures. We show that this approach makes it
possible to perform fast and accurate characterization of high-dimensional
programmable integrated interferometers. Moreover, the method provides access
to the physical characteristics of individual interferometer layers.

††journal: oe††articletype: Research Article

##  1 Introduction

The use of classical and quantum properties of light opens up new
possibilities for information processing methods [1]. Programmable optical
circuits play an important role in classical and quantum optical
communications [2, 3, 4, 5], in quantum computing [6, 7, 8], and in machine
learning problems [9, 10, 11]. The key element of such devices is a
programmable linear optical interferometer. When control parameters are fixed,
it performs a linear transformation of $N$ input fields:

$\matrixquantity(E_{1}^{\text{(out)}}&E_{2}^{\text{(out)}}&\cdots&E_{N}^{\text{(out)}})^{T}=V\cdot\matrixquantity(E_{1}^{\text{(in)}}&E_{2}^{\text{(in)}}&\cdots&E_{N}^{\text{(in)}})^{T},$
(1)

where $V$ is $N\times N$ complex transfer matrix of the interferometer. By
changing the values of the control parameters one can change $V$. This
dependence is determined by the interferometer internal structure, i.e. by the
position of beam splitters (both two- and multichannel) and phase shifters
[12, 13, 14, 15]. In the absence of losses, such structures are capable of
implementing an arbitrary unitary transformation $V$ of the fields. Moreover,
a proper choice of the architecture can achieve this even in the presence of
high manufacturing errors [14, 15, 16].

One can tune the interferometer to some target unitary transformation $V$ in
two ways. The first way is empirical. For each $V$, one optimizes the control
parameters in order to bring the experimental results of optical fields
measurements closer to the target ones [17]. This does not require building an
exact model of interferometer internal structure, but is impractical, since it
is necessary to measure all matrix elements at each step. For each new target
matrix $V$, the procedure must be repeated. Another way is to first build a
complete mathematical model of the interferometer based on the results of some
set of measurements [18]. Then, the search for the desired set of control
parameters is performed without any additional measurements. This procedure
faces some difficulties. Since large-scale linear optical devices can be
implemented in the integral design only, it is almost impossible to split the
large-scale reconstruction problem into several low-scale ones. It needs
additional optical outputs or modular optical circuits [19] that causes
additional losses and instability. Trying to characterize an integrated device
as a whole, one has to face the fact that the results of each measurement are
determined by a large number of internal device parameters. As a result, the
efficiency of the optimization procedure rapidly decreases with increasing
dimension $N$.

In this paper, we present a method for constructing a model of a programmable
integrated interferometer based on algebraic operations only. This approach
needs limited number of operations and is much more computationally efficient
than the one using numerical optimization. The method is based on the layered
interferometer architecture [15], which is described in Section 2. Further in
Section 3, we describe the procedure for reconstructing interferometer mixing
layers. As in [18], it requires the measurement of all matrix elements of the
complex transfer matrix $V$ (up to a global phase) for a given configuration
of phase layers. Common methods for measuring the transfer matrices of optical
interferometers using laser [20], thermal [21], and two-photon [22, 23] fields
make it possible to measure matrix elements up to input and output phases
only. A non-ambiguous measurement of the matrix $V$ can be performed using
homodyne detection of all the input and output coherent fields in (1) and
requires $O(N^{2})$ measurements [24]. In Section 4, we show that the proposed
algebraic approach allows one to predict the transformation matrix $V$ for an
arbitrary configuration of phase layers with an accuracy not lower than that
for the method based on numerical optimization. The mixing layers are
determined ambiguously, however, in the case of homogeneous losses, our method
makes it possible to uniquely determine the normalized transmission
coefficients of each individual layer. In addition, we show that the method is
highly robust to phase layers control errors.

##  2 Interferometer model

Below we consider the universal layered interferometer architecture consisting
of alternating mixing and phase layers [15, 18] (Fig. 1). The corresponding
transfer matrix has the form

$V=\Phi_{K+1}U_{K}\Phi_{K}\dots U_{2}\Phi_{2}U_{1}\Phi_{1}.$ (2)

Here $U_{k}$ and $\Phi_{k}$ are transfer matrices of the mixing and phase
layers, respectively. Each mixing layer is an $N$-channel transformation that
sufficiently distributes the signal of each input channel to all the output
channels. All matrices $U_{k}$ can be different. Ideally, they are unitary,
but this property is lost in the presence of linear losses. For the phase
layer, each channel is subjected to a controlled phase shift through a thermo-
or electro-optical effect. The corresponding transfer matrix has a diagonal
form:

$\Phi_{k}=\operatorname{diag}\quantity[\exp(i\varphi_{n}^{(k)})],\quad
n=1,\dots,N,\quad k=1,\dots,K+1.$ (3)

Expression (2) can also be used to describe any other interferometer
architecture (for example, the Reck [12] or Clements [13] architecture). In
this case, $U_{k}$ can be the transfer matrices of beam splitters acting on
two channels only, and the phase layers do not control all channels, but only
a few ones.

With the number of mixing layers $K=N$, the layered architecture makes it
possible to obtain an arbitrary transformation $V$ due to the proper choice of
phases $\varphi_{n}^{(k)}$. We mostly consider this particular case.

Figure 1: Schematic representation of the layered architecture of a
programmable interferometer. The number of mixing layers $K$ is equal to the
number of channels $N=5$.

The problem of constructing an interferometer model is reduced to determining
all transfer matrices $U_{k}$ and $\Phi_{k}$. The dependence of
$\varphi_{n}^{(k)}$ on the control parameters can be pre-calibrated even for
unknown $U_{k}$. For this, a simple $\sin\varphi$ dependence of the
interferometer output intensity on the value $\varphi$ of any single
controlled phase is used (all other phases are disabled) [2]. Below we assume
that all phase layers are pre-calibrated (later we will show that our
reconstruction method is robust to phase control errors). Reconstruction of
the matrices $U_{k}$ is complicated by the fact that there is no way to “turn
off” all mixing layers, except for one, and perform its measurement. However,
in the next section, we show that it is possible to obtain a matrix equation
that depends on a single mixing layer only by an appropriate choice of phase
layers configurations. This equation is reduced to an eigendecomposition
problem, which is solved algebraically.

##  3 Mixing layers reconstruction

In this section, we assume that we have full control over the phase layers.
Later in Section 4, we show that the proposed method is robust to control
errors. We also assume that none of the mixing layers has absolute losses in
any channel, so all the matrices $U_{1},\dots U_{K}$ are invertible. We also
do not restrict them to be unitary.

We describe the reconstructing algorithm in several steps. First, we consider
the ideal case when, for any configuration of the interferometer, we know
exactly its full transfer matrix $V$. Next, we consequentially take into
account the impossibility of experimental determination of the global phase of
$V$, the presence of errors in measurements and phase layers control, and the
special case of homogeneous losses inside the interferometer. Finally, we
present the pseudo-code for the complete algorithm. Its software
implementation is available online [25].

###  3.1 Ideal case

Consider two configurations of the interferometer phase layers. The first
configuration has all the phase layers disabled. The corresponding transfer
matrix is $V_{0}=U_{K}\dots U_{2}U_{1}$. Next, the phases
$\\{\varphi_{n},n=1,\dots,N\\}$ are set in the second phase layer:
$V_{1}=U_{K}\dots U_{2}\Phi U_{1}$, where
$\Phi=\operatorname{diag}[\exp(i\varphi_{n})]$. Consider the matrix
$A_{1}=V_{0}^{-1}V_{1}=U_{1}^{-1}\Phi U_{1}$. Since the matrix $\Phi$ is known
and diagonal, this equation is solved with respect to $U_{1}$ by the
diagonalization of $A_{1}$. By induction, one can iteratively obtain all
mixing matrices up to $U_{K-1}$. In particular, we enable only the $(k+1)$-th
phase layer to obtain the $k$-th mixing layer:

$V_{k}=U_{K}\dots U_{k+1}\Phi U_{k}\dots U_{2}U_{1}.$ (4)

Using $V_{k}$ and the results for all previous $k-1$ matrices we compute the
following matrix:

$A_{k}=U_{k-1}\dots U_{1}V_{0}^{-1}V_{k}U_{1}^{-1}\dots
U_{k-1}^{-1}=U_{k}^{-1}\Phi U_{k}.$ (5)

The following equation determines its eigendecomposition:

$A_{k}\vec{u}_{q}=\lambda_{q}\vec{u}_{q},\quad q=1,\dots,N.$ (6)

Since eigenvectors diagonalize $A_{k}$, we can use it to construct
$U_{k}^{-1}$. To do this, one must sort the solutions of (6) in such a way
that matches eigenvalues with the diagonal of $\Phi$. Let us define the
sorting indices $S=\\{s_{n}:\lambda_{s_{n}}=\exp(i\varphi_{n})\\}$. Then the
inverse transfer matrix has the form

$U_{k}^{-1}=\matrixquantity(\vec{u}_{s_{1}}&\vec{u}_{s_{2}}&\cdots&\vec{u}_{s_{N}}).$
(7)

This procedure goes successively for $k=1,\dots,K-1$. Finally, the transfer
matrix of the last mixing layer is $U_{K}=V_{0}U_{1}^{-1}U_{2}^{-1}\dots
U_{K-1}^{-1}$. The procedure does not involve the first and last phase layers.

Note that for the correct work of the described method it is necessary that
all $\lambda_{q}$ from (6) (so the phases $\varphi_{n}$) must be different.
Eigenvectors with equal eigenvalues cannot be uniquely defined, this leads to
errors in constructing the columns of $U_{k}^{-1}$.

###  3.2 Ambiguity

The described algorithm determines all mixing layers of the interferometer up
to following transformations:

$U_{1}\rightarrow\Lambda_{1}U_{1},\quad
U_{k}\rightarrow\Lambda_{k}U_{k}\Lambda_{k-1}^{-1}(1<k<K),\quad
U_{K}\rightarrow U_{K}\Lambda_{K-1}^{-1},$ (8)

where $\Lambda_{k}$ are arbitrary invertible diagonal matrices. Such ambiguity
is caused by the fact that the equation (6) is invariant under the
multiplication of the eigenvector by an arbitrary number. This, however, is
not the limitation of the proposed method, but the consequence of the
impossibility to address each individual mixing layer separately in the
experiment. Indeed, let us write the transfer matrix of the interferometer
with arbitrary control phases taking (8) into account:

$\displaystyle V$
$\displaystyle=\Phi_{K+1}(U_{K}\Lambda_{K-1}^{-1})\Phi_{K}\dots(\Lambda_{k}U_{k}\Lambda_{k-1}^{-1})\Phi_{k}\dots\Phi_{2}(\Lambda_{1}U_{1})\Phi_{1}$
(9)
$\displaystyle=\Phi_{K+1}U_{K}(\Lambda_{K-1}^{-1}\Lambda_{K-1})\Phi_{K}\dots(\Lambda_{k}^{-1}\Lambda_{k})\Phi_{k+1}U_{k}(\Lambda_{k-1}^{-1}\Lambda_{k-1})\Phi_{k}\dots(\Lambda_{1}^{-1}\Lambda_{1})\Phi_{1}U_{1}$
$\displaystyle=\Phi_{K+1}U_{K}\Phi_{K}\dots\Phi_{k+1}U_{k}\Phi_{k}\dots\Phi_{2}U_{1}\Phi_{1},$

which is equivalent to the initial transfer matrix. We have used the
commutativity of diagonal matrices $\Lambda_{k}$ and $\Phi_{k}$ in the second
line of the equation (9). Thus, transformations (8) cannot be observed in the
experiment, and they do not affect the ability of the model to predict the
transfer matrix for any configuration of phase layers. Physically this means
that we cannot distinguish the losses and constant phase shifts at the output
of one mixing layer from the losses and phase shifts at the input of the next
mixing layer (Fig. 2).

Figure 2: Illustration of the impossibility to distinguish the losses and
phase shifts at the output of one mixing layer from the losses and phase
shifts at the input of the next mixing layer. Boxes with symbols $\theta$ and
$\eta$ denote phase shifts and linear losses, respectively.

###  3.3 Effect of the global phase

The approach described above gives an exact solution, when the matrices
$V_{k}$ are exactly known. In practice, the elements of the transfer matrices
are measured experimentally and, therefore, differ from the original ones.
These differences are primarily related to the undetectable global phase. In
this regard, in equation (6), we obtain eigenvalues
$\lambda_{q}=\exp(i\delta_{q})$, where the set of phases $\\{\delta_{q}\\}$
differs from the original set $\\{\varphi_{n}\\}$ by some constant value
$\chi$. It makes the matching between $\\{\delta_{q}\\}$ and
$\\{\varphi_{n}\\}$ ambiguous and could violate the sorting of eigenvectors.
For example, a set of phases

$\varphi_{n}=\frac{2\pi(n-1)}{N},\quad n=1,\dots,N$ (10)

from Fig. 3a is invariant under $2\pi/N$ rotation and can’t be used to
determine $\chi$ in a unique way. This problem is solved by an appropriate
choice of controllable phases such that $\\{\varphi_{n}+\chi\ (\mathrm{mod}\
2\pi)\\}\neq\\{\varphi_{n}\\}$ for any non-trivial $\chi$. Consider the
example presented in Fig. 3b:

$\varphi_{n}=\frac{2\pi(n-1)}{N+1},\quad n=1,\dots,N.$ (11)

One can observe that this set is invariant under the full $2\pi$ rotations
only. Further, it is convenient to move from absolute phases to relative
phases between adjacent channels of the phase layer (channel $N$ is considered
adjacent to channel $1$):

$\Delta\varphi_{n}=\begin{cases}\varphi_{1}-\varphi_{N}\ (\mathrm{mod}\
2\pi),&n=1\\\ \varphi_{n}-\varphi_{n-1}\ (\mathrm{mod}\
2\pi),&n>1.\end{cases}$ (12)

It is easy to see that the relative phase for the first channel is twice the
others for the choice (11):

$\Delta\varphi_{1}=\frac{4\pi}{N+1},\quad\Delta\varphi_{n>1}=\frac{2\pi}{N+1}.$
(13)

We can then calculate the sorting indices of eigenvalues from (6) by relative
phases that are insensitive to the global phase $\chi$:

$S=\quantity{s_{n}:\Delta\delta_{s_{1}}=\Delta\varphi_{1};\;\delta_{s_{n>1}}=\delta_{s_{n-1}}+\Delta\varphi_{n}}.$
(14)

Figure 3: The distribution of phases on the unit circle. (a) Uniform
distribution of controllable phases on the single phase layer for $N=5$. If
the constant $\chi$ is a multiple of $2\pi/5$ then the transformation
$\varphi_{n}\rightarrow\varphi_{n}+\chi\ (\mathrm{mod}\ 2\pi)$ leads to a
similar distribution making it impossible to unambiguously compare the
eigenvalues $\lambda_{q}=\exp(i\delta_{q})$ with the phases $\varphi_{n}$. (b)
The choice of asymmetric phase distribution (11) distinguishes $\varphi_{1}$
from the others by the greatest distance to the previous phase. (c) The search
for the largest distance between adjacent phases of eigenvalues $\lambda_{q}$
allows one to determine the index $q$ corresponding to the phase $\varphi_{1}$
and retrieve the constant $\chi$. (d) Strong phase fluctuations of the
eigenvalues $\lambda_{q}$ at large $N$ (in this case $N=19$) violates the
sorting of the eigenvectors.

Systematic and statistical errors also distort the spectrum of the matrix
$A_{k}$, leading to more general eigenvalues:
$\lambda_{q}=r_{q}\exp(i\delta_{q})$, where $r_{q}\approx 1$,
$\delta_{q}\approx\varphi_{q}+\chi$. If the errors are not too large then we
can again use the above approach, but instead of an exact equality we look for
an approximate one:

$S=\quantity{s_{n}:s_{1}=\operatorname*{arg\,max}_{q}\Delta\delta_{q};\;s_{n>1}=\operatorname*{arg\,min}_{q}\quantity[\delta_{q}-\delta_{s_{n-1}}\
(\mathrm{mod}\ 2\pi)]}.$ (15)

Note that in practice this minimization is not necessary. The set (11) is
sorted in ascending order, so one could also sort all $\delta_{q}$ and then
find the maximum value of $\Delta\delta_{q}$. In view of (13), this value
corresponds to the index $s_{1}$. The remaining elements will follow it on the
unit circle in the anti-clockwise direction (Fig. 3c).

###  3.4 Noise handling

The described algorithm works well on small dimensions. However, as $N$
increases, the density of $\\{\varphi_{n}\\}$ values on the unit circle
increases. As a result, the eigenvalue sorting rule becomes more noise
sensitive (Fig. 3d). One can reduce this effect by activating only a few
phases at a time and calculating only a few columns of the matrix
$U_{k}^{-1}$. Denote by $C=\\{c_{m}\\}$ the set of indices of these columns
and select the following set of phases (Fig. 4a):

$\varphi_{n}=\begin{cases}0,&n\notin C,\\\ \frac{2\pi
m}{M+2},&n=c_{m},\end{cases}\quad n=1,\dots,N,$ (16)

where $M<N$ is the size of $C$. The set (16) contains $M+1$ unique values.
After measuring the interferometer transfer matrix $V_{k}$ and obtaining the
matrix $A_{k}$ from (5), its eigenvalues are sorted according to rule (15). It
turns out that the values $\delta_{s_{n}}$ for $1\leq n\leq N-M$ are close to
each other, and $\Delta\delta_{s_{1}}$ is maximal (Fig. 4b). The remaining $M$
values correspond to the desired eigenvectors (Fig. 4c).

Figure 4: An example of a partial reconstruction of the inverse matrix of the
$k$-th mixing layer. The dimension of the interferometer is $N=9$. $M=4$
matrix columns from 5 to 8 are estimated. (a) The $k$-th phase layer
configuration. (b) The phases of matrix $A_{k}$ eigenvalues are sorted on the
unit circle in ascending order. The first element is considered to be the one
that is separated from the previous one by the largest distance. The first
$N-M=5$ elements correspond to eigenvectors with indices $n\notin C$. The
remaining $M=4$ elements are eigenvectors with indices $n\in C$. (c) The
corresponding eigenvectors are set to be the columns of the matrix
$U_{k}^{-1}$.

The smaller the number $M$ of simultaneously estimated columns of
$U_{k}^{-1}$, the less densely the corresponding phase values can be located.
As a result it is possible to resolve the individual eigenvalues of the matrix
$A_{k}$ from (6) quite well even for large values of dimension $N$. In the
limit $M=1$, we reconstruct the columns of the matrix $U_{k}^{-1}$ one by one.
This, however, complicates the experiment: it is required to measure $\propto
N/M$ transfer matrices for each mixing layer. The total number of transfer
matrices used is $L=(K-1)\lceil N/M\rceil+1$.

Note that in the above algorithm, the phases $\delta_{q}$ of $A_{k}$
eigenvalues were used for the correct eigenvectors sorting only. At the same
time, the absolute values of $\delta_{q}$ contain information about the true
values of the controlled phases $\varphi_{n}$ at the phase layer. These values
can serve for more accurate calibration of the phase layer, taking its cross-
talk into account.

###  3.5 The case of homogeneous losses

The described procedure does not impose any constraints on resulting matrices
$U_{k}$. In certain cases, it may be known a priory that the interferometer
losses are homogeneous. Then each $U_{k}$ belongs to the set of unitary
matrices up to a constant factor $\eta_{0}$ describing the homogeneous losses.
Thus, one can enhance the accuracy by projecting each $U_{k}$ to this set. For
some non-unitary matrix $X$, the projection is given by

$\operatorname{proj}[X]=QW,\quad X=QSW^{\dagger}.$ (17)

Here matrices $Q$, $S$ and $W$ form the singular value decomposition of $X$.
The matrices $Q$ and $W$ are unitary, while $S$ is non-negative diagonal
matrix and equals to the identity matrix for unitary $X$. Within the unitary
model, the transformations $\Lambda_{k}$ from (8) are limited by arbitrary
phase shifts only (there are no linear loss blocks in Fig. 2). This ambiguity
does not affect the normalized transmission coefficients of mixing layers

$[T_{k}]_{mn}=\frac{\absolutevalue{[U_{k}]_{mn}}^{2}}{\frac{1}{N}\Tr(U_{k}^{\dagger}U_{k})}=\frac{1}{\eta_{0}^{2}}\absolutevalue{[U_{k}]_{mn}}^{2}.$
(18)

Thus, this coefficients are defined unambiguously. Below we consider the case
$\eta_{0}=1$ since homogeneous losses could be measured separately and
compensated by collecting more data.

###  3.6 Algorithm pseudo-code

Algorithms 1 summarizes the described procedure of mixing layers
reconstructions. It takes the following parameters as the input: $N$ – the
number of input and output channels of the interferometer; $K$ – the number of
mixing layers; $M$ – the maximal number of $U_{k}^{-1}$ columns, estimated at
once. The output of the algorithm are transfer matrices $U_{k}$
($k=1,\dots,K$) of all mixing layers.

The algorithm software implementation is available online [25].

Algorithm 1 Mixing layers reconstruction

Measure $V_{0}$ (all phase layers are disabled)

$V_{T}\leftarrow I$ $\triangleright$ The product of all currently estimated
$U_{k}$

for $k=1,\dots,K-1$ do

$U_{k}^{(\text{inv})}\leftarrow I$ $\triangleright$ Initial estimation

$j\leftarrow 0$ $\triangleright$ Current number of estimated columns

while $j<N$ do

$M^{\prime}=\text{min}(M,N-j)$

if $M^{\prime}=N$ then $\triangleright$ All columns of $U_{k}^{-1}$ are
estimated at once

$\varphi_{n}\leftarrow\frac{2\pi(n-1)}{N+1}$ ($n=1,\dots,N$)

else$\triangleright$ Only $M^{\prime}$ columns of $U_{k}^{-1}$ are estimated

$\varphi_{n}\leftarrow 0$ ($n=1,\dots,N$)

$\varphi_{j+m}\leftarrow\frac{2\pi m}{M^{\prime}+2}$ ($m=1,\dots,M^{\prime}$)

end if

Measure $V_{k}$ (phases of $(k+1)$-th phase layer are $\varphi_{n}$)

$A_{k}\leftarrow V_{T}V_{0}^{-1}V_{k}V_{T}^{-1}$

Solve $A_{k}\vec{u}_{q}=r_{q}e^{i\delta_{q}}\vec{u}_{q}$ for $\vec{u}_{q}$,
$r_{q}\geq 0$, $\delta_{q}\in[0,2\pi)$ ($q=1,\dots,N$)

Find sorting indices $S$ using (15)

$[U_{k}^{(\text{inv})}]_{n,j+m}\leftarrow[\vec{u}_{s_{N-M^{\prime}+m}}]_{n}$
($m=1,\dots,M^{\prime}$) $\triangleright$ Update $M^{\prime}$ columns of
$U_{k}^{-1}$ to match sorted eigenvectors of $A_{k}$

$j\leftarrow j+M^{\prime}$

end while

$U_{k}\leftarrow[U_{k}^{(\text{inv})}]^{-1}$

if Using unitary model then

$U_{k}\leftarrow\text{proj}[U_{k}]$ using (17)

end if

$V_{T}\leftarrow U_{k}V_{T}$

end for

$U_{K}\leftarrow V_{0}V_{T}^{-1}$

##  4 Numerical simulation

In this section, we analyze the efficiency of the proposed algorithm in the
presence of various kinds of errors. In addition, we compare it with the
optimization algorithm from [18].

###  4.1 Error models

We define random mixing layers in two steps. In the first step, we generate a
random unitary matrix

$U_{H}^{\gamma}=\operatorname{proj}[U_{H}+\gamma G],\quad
G_{mn}\sim\operatorname{Normal}(0,1)+i\operatorname{Normal}(0,1),\quad
m,n=1,\dots,N,$ (19)

where $\operatorname{proj}[\cdot]$ is the unitary projection operation (17),
and $\gamma$ characterizes the degree of distinction of $U_{H}^{\gamma}$ from
the $N$-dimensional Hadamard transform, proportional to the discrete Fourier
transform matrix (the reconstruction method is not restricted to this
particular type of mixing layers and could by applied to any invertable
matrix):

$[U_{H}]_{mn}=\frac{1}{\sqrt{N}}\exp(-i\frac{2\pi(m-1)(n-1)}{N}),\quad
m,n=1,\dots,N.$ (20)

Further we assume $\gamma=0.1$. In the second step, if necessary, we introduce
heterogeneous linear losses:

$U_{H}^{\gamma,\beta}=\sqrt{U_{H}^{\gamma}}R_{\beta}\sqrt{U_{H}^{\gamma}},\quad
R_{\beta}=\operatorname{diag}[1-\beta r_{n}],\quad
r_{n}\sim\operatorname{Uniform}(0,1),\quad n=1,\dots,N.$ (21)

In other words, we insert random linear losses of no more than $\beta$ into
the “middle” part of each random unitary transformation. Random matrices of
the form (21) are generated independently for each mixing layer to create a
true model of the programmable interferometer. We denote these matrices as
$U_{k}^{\text{true}}$.

We also simulate the errors of phase layers control. Let $\\{\varphi_{n}\\}$
be the desired set of phases in the $k$-th phase layer. The relative control
error of each phase is given by $\alpha$:

$\varphi_{n}^{\alpha}=\varphi_{n}(1+\alpha g_{n}),\quad
g_{n}\sim\operatorname{Normal}(0,1),\quad n=1,\dots,N.$ (22)

We simulate the cross-talk between the phases of a the single phase layer as
following:

$\varphi_{n}^{\alpha,w}=\sum_{m=1}^{N}{\varphi_{n}^{\alpha}\exp(-\frac{(n-m)^{2}}{2w^{2}})},\quad
n=1,\dots,N,$ (23)

where $w$ characterizes the cross-talk strength. The errors of the type (22)
are simulated independently for each phase layer each time the control phases
of the interferometer are changed. An important condition is the absence of
cross-talk between different phase layers. This condition is usually satisfied
in the experiment, since the corresponding control elements (electrodes or
heaters) are usually quite far apart on the optical chip.

We simulate the statistical noise of the transfer matrix measurements by a
random Gaussian error for each matrix element:

$V\rightarrow V+\varepsilon\sqrt{L}G,\quad
G_{mn}\sim\operatorname{Normal}(0,1)+i\operatorname{Normal}(0,1),\quad
m,n=1,\dots,N,$ (24)

where $\varepsilon$ is the noise level, and $L$ is the measurement protocol
size (the number of transfer matrices to measure). The $\varepsilon$ value
includes all the statistical errors in the transfer matrix measurement
procedure (e. g. described in [24]) and depends on the experiment details. The
correction $\sqrt{L}$ is based on the standard assumption that statistical
errors are inversely proportional to the square root of the acquisition time.
It allows to compare different methods under the conditions of a fixed data
acquisition time. For example, let Method A requires to measure $50$ transfer
matrices and Method B requires $100$. Then, if the experiment duration is
limited, one needs to reduce the measurement time twice for each matrix for
Method B (we neglect the phase layers configuration switching time). This
relatively increases the variance of matrix elements estimates by a factor of
$2$ for Method A. The linear error increases $\sqrt{2}$ times. The error (24)
is generated independently each time, when the algorithms need to obtain a new
transfer matrix.

###  4.2 Benchmarks

After the numerical experimental data generation we reconstruct the
interferometer model. Then we estimate the ability of the model to predict the
full transfer matrix given arbitrary phase layer configuration (each phase of
each phase layer is chosen randomly from $0$ to $2\pi$) by calculating the
fidelity between the predicted $V_{\text{pred}}$ and true $V_{\text{true}}$
transfer matrices of the entire interferometer (the normalized squared dot
product of the matrices):

$F=\frac{\absolutevalue{\Tr(V_{\text{true}}^{\dagger}V_{\text{pred}})}^{2}}{\Tr(V_{\text{true}}^{\dagger}V_{\text{true}})\Tr(V_{\text{pred}}^{\dagger}V_{\text{pred}})}.$
(25)

The quantity $\Delta F=1-F$, in turn, plays the role of the distance between
the matrices, which is insensitive to the value of the undetectable global
phase (similar to the distance between two quantum processes [26]). We
consider the $95$-th percentile $\Delta F_{95}$ in a test set of $1000$ random
phase layer configurations [27]. Note that the value of $\Delta F_{95}$ is
affected not only by the quality of the mixing layers reconstruction, but also
by the presence of phase layer control errors in $V_{\text{true}}$, which are
assumed to be zero in $V_{\text{pred}}$.

In Section 3.5, we pointed out that the method makes it possible to uniquely
determine the normalized transmission coefficients (18) of mixing layers.
Therefore, we also consider the measure

$\Delta T_{\max}=\max_{\begin{subarray}{c}k=1,\dots,K\\\
m,n=1,\dots,N\end{subarray}}{\absolutevalue{[T_{k}^{\text{true}}]_{mn}-[T_{k}^{\text{rec}}]_{mn}}}.$
(26)

It characterizes the maximal error of transmission coefficients estimation
over the all mixing layers. Here $T_{k}^{\text{true}}$ and
$T_{k}^{\text{rec}}$ are respectively the true and reconstructed transmission
matrices of the $k$-th mixing layer. Both $\Delta F_{95}$ and $\Delta
T_{\max}$ are insensitive to homogeneous losses. Considering $\Delta T_{\max}$
in addition to $\Delta F_{95}$ is useful for two reasons. First, $\Delta
F_{95}$ is an integral estimate that includes both errors in determining
mixing layers and errors in controlling phase layers. Since the latter remains
unknown, $\Delta F_{95}$ may be quite high even if each mixing layer is
determined accurately. Secondly, the determination of the transmission
coefficient is of independent interest for debugging the manufacturing
technology of integrated interferometers. The use of $\Delta T_{\max}$,
however, makes sense for the unitary model only (see Section 3.5). In the
presence of heterogeneous losses, the matrices $T_{k}^{\text{rec}}$ are
determined ambiguously and it makes no sense to compare them directly with
$T_{k}^{\text{true}}$.

###  4.3 Comparison with optimization method

Let us consider the case when the number of interferometer channels $N$ and
the number of mixing layers $K$ are 5 ($N=K=5$), and there are no linear
losses ($\beta=0$) and phase control errors ($\alpha=w=0$). We generate a
random programmable interferometer according to Section 4.1 and perform the
reconstruction of all its mixing layers according to Section 3. As in [18], we
use the model of unitary mixing layers. The initial data used in the algorithm
is also fed to the input of the algorithm from [18] with an open source code.
We use the version of this algorithm that uses the approximate prior knowledge
of mixing layers transfer matrices. Following the example from the library
files, the following algorithm parameters are chosen: mini-batch size is five,
the prior knowledge coefficient $\beta=0.24$. The number of training epochs
varies from $20$ to $1500$ depending on the dimension $N$ and the measurement
protocol (the value is high enough to reliably rich the accuracy saturation).

Figure 5: The errors in determining the transmission coefficients $\Delta
T_{\max}$ (a,c) and predicting the full transfer matrices $\Delta F_{95}$
(b,d). The average values over 100 independent trials are shown. Dimension
$N=5$, number of mixing layers $K=5$, no losses ($\beta=0$). The results for
our algorithm (solid lines) and for the algorithm from [18] (diamonds) are
shown. (a,b) Dependence on the measurement error $\varepsilon$; no phase
control errors ($\alpha=w=0$). (c,d) Dependence on the phase control error
$\alpha$; measurement error $\varepsilon=10^{-3}$, cross-talk strength
$w=0.1$.

In total, we considered 3 different measurement protocols:

  1. 1.

First, all phases are set to zero. Next, $20$ configurations are considered.
In each one a single non-zero phase is set equal to $2\pi/3$. Each of these
$20$ configurations allows one to define $M=1$ column of the inverse transfer
matrix of one mixing layer (see Section 3.4). The total number of
configurations is $L=21$.

  2. 2.

First, all phases are set to zero. Next, $4$ configurations are considered. In
each one the phases (11) are set on one of the phase layers. Each of these $4$
configurations allows one to determine all $M=5$ columns of the inverse
transfer matrix of one mixing layer (see Section 3.4). The total number of
configurations is $L=5$.

  3. 3.

A set of $L=300$ random interferometer configurations is considered (only for
the optimization algorithm from [18]).

Fig. 5 shows the simulation results. In the absence of phase layers control
errors, the dependencies of $\Delta T_{\max}$ (Fig. 5a) and $\Delta F_{95}$
(Fig. 5b) on the measurement error $\varepsilon$ turned out to be close,
regardless of the measurement protocol and the reconstruction method used. The
closeness is so high that plotting them in the same axes on a log scale make
the distinction barely seen. Approximately, we have observed $\Delta
T_{\max}\propto\varepsilon$ and $\Delta F_{95}\propto\varepsilon^{2}$. From
Fig. 5c it can be seen that the presence of phase layer control errors
$\alpha$ does not affect our method’s ability to estimate the mixing layer
transmission coefficients $[T_{k}]_{mn}$: it is limited by the interferometer
measurements accuracy only. This is because the algorithm uses the phases
values only to do the correct sorting. However, starting from a certain
critical value of the error $\alpha_{c}$, the accuracy of determining
$[T_{k}]_{mn}$ sharply decreases, since the sorting is violated. The value of
$\alpha_{c}$ can be increased by decreasing the number of unique phases at
each iteration (it is defined by the parameter $M$ of the algorithm). The
optimization algorithm from [18] relies on the particular phase values, so the
control errors impair the estimation of $[T_{k}]_{mn}$. Despite this, the
algorithm manages to train the entire interferometer as a whole with an error
$\Delta F_{95}$ not higher than that provided by our algorithm (Fig. 5d).
Approximately, we have observed $\Delta F_{95}\propto\alpha^{2}$.

Figure 6: Comparison of the computation time for different methods versus the
dimension $N$. The average values over 100 independent trials are shown. The
number of mixing layers $K=N$, no losses ($\beta=0$), no phase control errors
($\alpha=w=0$), transfer matrices measurement error $\varepsilon=10^{-5}$.

Despite the closeness of the results for optimization algorithm and our
algorithm, the latter performs much faster. To show this, consider two extreme
versions of our algorithm ($M=N$ and $M=1$). For the optimization algorithm
from [18], we use a set of 300 training phase layer configurations. Starting
from a certain value, changing this number does not affect the computation
speed significantly. The number of epochs is chosen to be minimal in order to
guarantee a correct result. On the range of considered dimensions $N$, the
presented reconstruction algorithm worked more than two orders of magnitude
faster (Fig. 6). We also considered the optimization algorithm in the “black
box” model (the absence of any prior knowledge about the mixing layers). In
this case, as expected, there was a rapid growth in computational complexity
with an increase in $N$ from two to six (as in [18], we also failed to obtain
the convergence of the optimization algorithm in the “black box” model for
$N>6$). Our algorithm, also based on the “black box” model, showed a much
slower increase in computation time. The independent trials of numerical
experiments were performed in parallel on 52 CPU cluster; Intel(R) Xeon(R) CPU
E5450 @ 3.00GHz.

###  4.4 Algorithm scalability

Above we have mentioned that the total number of required interferometer
configurations for our method is $L=(K-1)\lceil N/M\rceil+1$, where $N$ is the
number of input and output channels, $K$ is the number of mixing layers, and
$M\in[1;N]$ is the number of simultaneously estimated columns of the inverse
matrix of a single mixing layer. The time required to measure the transfer
matrix for each configuration is estimated as $O(N^{2})$ [24]. Therefore, for
$K=N$ the total measurement time scales as $O(N^{4}/M)$.

The reconstruction algorithm requires $O(L)$ algebraic operations including
matrix multiplication, taking the inverse matrix, and computing matrix
eigendecomposition. The complexity of these procedures does not exceed
$O(N^{3})$ [28]. For $K=N$ the total complexity of the entire procedure can be
estimated as $O(N^{5}/M)$.

Figure 7: (a) Scalability of the errors $\Delta T_{\max}$ in determining the
mixing layers transmission coefficients. The average values over 100
independent trials are shown. Number of mixing layers $K=N$, no losses
($\beta=0$), cross-talk strength $w=0.1$, transfer matrices measurement error
$\varepsilon=10^{-5}$. (b) Distribution of transmission coefficients errors
for various mixing layers within a single trial; confidence intervals show the
upper and lower quartiles over all $m,n=1,\dots,N$. Parameters: $N=K=25$,
$\beta=0$, $w=0.1$, $\alpha=10^{-3}$, $\varepsilon=10^{-5}$. (c) Scalability
of the errors $\Delta F_{95}$ in predicting the full transfer matrices. The
average values over 100 independent trials are shown. Parameters: $K=N$,
$M=N$, $\alpha=w=0$.

From these estimates, one can see the number $M$ to be an important parameter
that affects the algorithm speed. The value $M=N$ is optimal from the point of
view of measurement and computation speed. However, in this case, with
increasing $N$, the sensitivity to control errors of phase layers increases:
if the error $\alpha$ is higher than the critical value $\alpha_{c}$, the
mixing layers estimates accuracy sharply decreases. At the same time, for a
fixed $M$, the value of $\alpha_{c}$ changes slightly with increasing $N$
(Fig. 7a). Note also that for $\alpha<\alpha_{c}$ transmission coefficients
matrices $T_{k}$ are determined with approximately the same accuracy for all
mixing layers (Fig. 7b). That is, despite the iterative nature of the
algorithm (see eq. (5)), there is no effect of error accumulation.

The prediction error $\Delta F_{95}$ grows polynomially with increasing $N$
(Fig. 7c). This is primarily because at a fixed data acquisition time, the
size $L$ of the measurement protocol increases. As a result, the measurement
accuracy of each transfer matrix decreases (see eq. (24)), impairing the
mixing layers estimation accuracy. The figure shows that one can slightly
improve the accuracy by performing the projection of the mixing layers to the
set of unitary matrices (Section 3.5). This, however, does not affect the rate
of growth (about $\Delta F_{95}\propto N^{4}$). When the heterogeneous linear
losses are present, the rate of growth becomes higher (about $\Delta
F_{95}\propto N^{5.8}$ for $\beta=0.1$). This can be explained by the fact
that the condition number of the transfer matrices of the mixing layers
increases when heterogeneous losses occur. As a result, the matrix inversion
procedure becomes more sensitive to statistical fluctuations.

##  5 Conclusion

We have shown that the estimation of mixing layers in the layered architecture
of the programmable integrated optical interferometer is possible without any
use of numerical optimization. The proposed algorithm has a single parameter
$M$ –- the number of simultaneously estimated columns of a mixing layer
inverse transfer matrix. We have performed the numerical simulation and
compared the results with the results obtained using the optimization
algorithm. It has been found that under the equal conditions both methods
provide the same accuracy, regardless of the measurement protocol. At the same
time, the exclusively algebraic approach in our algorithm has shown a
significantly lower computation time.

Moreover, our algorithm uses the control phases values only to correctly sort
the columns of the adjacent mixing layer. This makes the method almost
insensitive to phase control errors up to some critical error value. This
critical value can be increased by choosing a smaller value of the parameter
$M$. Small $M$, however, increases the number of operations in the
reconstruction algorithm, so in practice it is worth choosing an optimal $M$
based on the experimental accuracy of phases control. We believe that
extending the algorithm to consider the absolute values of the phases can help
in determining the phase control errors. This is the subject of further
research.

Finally, note that any interferometer architecture (for example, the Reck [12]
or Clements [13] architecture) can be brought to the layered architecture
considered here by adding controllable phase shifters in some channels.
However, the possibility to directly apply our method to an arbitrary
architecture remains unclear and is also the subject of further research.

Acknowledgments This work was supported by Russian Science Foundation (RSF),
project No: 22-12-00263. We also thank I. Dyakonov for valuable discussions
and G. Struchalin for help in conducting computations.

Disclosures The authors declare no conflicts of interest.

Data Availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

## References

  * [1] P. Minzioni, C. Lacava, T. Tanabe _et al._ , “Roadmap on all-optical processing,” Journal of Optics 21, 063001 (2019).
  * [2] J. Carolan, C. Harrold, C. Sparrow, E. Martín-López, N. J. Russell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh, G. D. Marshall, M. G. Thompson, J. C. F. Matthews, T. Hashimoto, J. L. O’Brien, and A. Laing, “Universal linear optics,” Science 349, 711–716 (2015).
  * [3] N. C. Harris, J. Carolan, D. Bunandar, M. Prabhu, M. Hochberg, T. Baehr-Jones, M. L. Fanto, A. M. Smith, C. C. Tison, P. M. Alsing, and D. Englund, “Linear programmable nanophotonic processors,” Optica 5, 1623–1631 (2018).
  * [4] J. Wang, S. Paesani, Y. Ding, R. Santagati, P. Skrzypczyk, A. Salavrakos, J. Tura, R. Augusiak, L. Mančinska, D. Bacco, D. Bonneau, J. W. Silverstone, Q. Gong, A. Acín, K. Rottwitt, L. K. Oxenløwe, J. L. O’Brien, A. Laing, and M. G. Thompson, “Multidimensional quantum entanglement with large-scale integrated optics,” Science 360, 285–291 (2018).
  * [5] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, “Integrated photonic quantum technologies,” Nature Photonics 14, 273–284 (2020).
  * [6] I. Schwartz, D. Cogan, E. R. Schmidgall, Y. Don, L. Gantz, O. Kenneth, N. H. Lindner, and D. Gershoni, “Deterministic generation of a cluster state of entangled photons,” Science 354, 434–437 (2016).
  * [7] H.-S. Zhong, Y. Li, W. Li _et al._ , “12-photon entanglement and scalable scattershot boson sampling with optimal entangled-photon pairs from parametric down-conversion,” Physical review letters 121, 250505 (2018).
  * [8] W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J.-i. Yoshikawa, N. C. Menicucci, H. Yonezawa, and A. Furusawa, “Generation of time-domain-multiplexed two-dimensional cluster state,” Science 366, 373–376 (2019).
  * [9] R. Hamerly, L. Bernstein, A. Sludds, M. Soljačić, and D. Englund, “Large-scale optical neural networks based on photoelectric multiplication,” Physical Review X 9, 021032 (2019).
  * [10] G. Wetzstein, A. Ozcan, S. Gigan, S. Fan, D. Englund, M. Soljačić, C. Denz, D. A. Miller, and D. Psaltis, “Inference in artificial intelligence with deep optics and photonics,” Nature 588, 39–47 (2020).
  * [11] H. Zhang, M. Gu, X. D. Jiang, J. Thompson, H. Cai, S. Paesani, R. Santagati, A. Laing, Y. Zhang, M. H. Yung, Y. Z. Shi, F. K. Muhammad, G. Q. Lo, X. S. Luo, B. Dong, D. L. Kwong, L. C. Kwek, and A. Q. Liu, “An optical neural chip for implementing complex-valued neural network,” Nature communications 12, 1–11 (2021).
  * [12] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, “Experimental realization of any discrete unitary operator,” Physical review letters 73, 58 (1994).
  * [13] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, “Optimal design for universal multiport interferometers,” Optica 3, 1460–1465 (2016).
  * [14] S. A. Fldzhyan, M. Y. Saygin, and S. P. Kulik, “Optimal design of error-tolerant reprogrammable multiport interferometers,” Optics Letters 45, 2632–2635 (2020).
  * [15] M. Y. Saygin, I. Kondratyev, I. Dyakonov, S. Mironov, S. Straupe, and S. Kulik, “Robust architecture for programmable universal unitaries,” Physical Review Letters 124, 010501 (2020).
  * [16] S. Pai, B. Bartlett, O. Solgaard, and D. A. Miller, “Matrix optimization on universal unitary photonic devices,” Physical review applied 11, 064044 (2019).
  * [17] T. W. Hughes, M. Minkov, Y. Shi, and S. Fan, “Training of photonic neural networks through in situ backpropagation and gradient measurement,” Optica 5, 864–871 (2018).
  * [18] S. Kuzmin, I. Dyakonov, and S. Kulik, “Architecture agnostic algorithm for reconfigurable optical interferometer programming,” Optics Express 29, 38429–38440 (2021).
  * [19] P. L. Mennea, W. R. Clements, D. H. Smith, J. C. Gates, B. J. Metcalf, R. H. S. Bannerman, R. Burgwal, J. J. Renema, W. S. Kolthammer, I. A. Walmsley, and P. G. R. Smith, “Modular linear optical circuits,” Optica 5, 1087–1090 (2018).
  * [20] S. Rahimi-Keshari, M. A. Broome, R. Fickler, A. Fedrizzi, T. C. Ralph, and A. G. White, “Direct characterization of linear-optical networks,” Optics express 21, 13450–13458 (2013).
  * [21] K. Katamadze, G. Avosopiants, A. Romanova, Y. I. Bogdanov, and S. Kulik, “Linear optical circuits characterization by means of thermal field correlation measurement,” Laser Physics Letters 18, 075201 (2021).
  * [22] A. Laing and J. L. O’Brien, “Super-stable tomography of any linear optical device,” arXiv preprint arXiv:1208.2868 (2012).
  * [23] A. Peruzzo, A. Laing, A. Politi, T. Rudolph, and J. L. O’brien, “Multimode quantum interference of photons in multiport integrated devices,” Nature communications 2, 1–6 (2011).
  * [24] K. V. Jacob, A. E. Mirasola, S. Adhikari, and J. P. Dowling, “Direct characterization of linear and quadratically nonlinear optical systems,” Physical Review A 98, 052327 (2018).
  * [25] B. Bantysh, “The package for simulating and training integrated linear optical devices,” https://github.com/PQCLab/ILOptics.
  * [26] A. Gilchrist, N. K. Langford, and M. A. Nielsen, “Distance measures to compare real and ideal quantum processes,” Physical Review A 71, 062310 (2005).
  * [27] B. I. Bantysh, A. Y. Chernyavskiy, and Y. I. Bogdanov, “Quantum tomography benchmarking,” Quantum Information Processing 20, 1–20 (2021).
  * [28] J. W. Demmel, _Applied numerical linear algebra_ (Society for Industrial and Applied Mathematics, 1997).


# Matter-wave chaos with a cold atom in a standing-wave laser field

S.V. Prants<EMAIL_ADDRESS>tel.007-4232-312602, fax 007-4232-312573
Laboratory of Nonlinear Dynamical Systems, Pacific Oceanological Institute of
the Russian Academy of Sciences, Baltiiskaya St., 43, 690041 Vladivostok,
Russia

###### Abstract

Coherent motion of cold atoms in a standing-wave field is interpreted as a
propagation in two optical potentials. It is shown that the wave-packet
dynamics can be either regular or chaotic with transitions between these
potentials after passing field nodes. Manifestations of de Broglie-wave chaos
are found in the behavior of the momentum and position probabilities and the
Wigner function. The probability of those transitions depends on the ratio of
the squared detuning to the Doppler shift and is large in that range of the
parameters where the classical motion is shown to be chaotic.

###### keywords:

cold atom; matter wave; quantum chaos

###### PACS:

03.75.-b, 42.50.Md

##  1 Introduction

Cold atoms in a standing-wave laser field [1, 2, 3, 4, 5] are ideal objects
for studying fundamental principles of quantum physics, quantum-classical
correspondence, and quantum chaos. The proposal [6] to study atomic dynamics
in a far-detuned modulated standing wave made atomic optics a testing ground
for quantum chaos. A number of impressive experiments have been carried out in
accordance with this proposal [7, 8, 9]. Furthermore, cold atoms have been
recently used to study different phenomena in statistical physics including
ratchet effect of a directed transport of atoms in the absence of a net bias
[10, 11, 12] and Brillouin-like propagation modes in optical lattices [13].

New possibilities are opened if one works near the atom-field resonance where
the interaction between the internal and external atomic degrees of freedom is
intense. A number of nonlinear dynamical effects have been found [14, 15, 16]
with point-like atoms in a near-resonant rigid optical lattice: chaotic Rabi
oscillations, chaotic walking, dynamical fractals, Lévy flights, and anomalous
diffusion.

Dynamical chaos in classical mechanics is a special kind of random-like motion
without any noise and/or random parameters. It is characterized by exponential
sensitivity of trajectories in the phase space to small variations in initial
conditions and/or control parameters. Such sensitivity does not exist in
isolated quantum systems because their evolution is unitary, and there is no
well-defined notion of a quantum trajectory. Thus, there is a fundamental
problem of emergence of classical dynamical chaos from more profound quantum
mechanics which is known as quantum chaos problem and the related problem of
quantum-classical correspondence. In a more general context it is a problem of
wave chaos. It is clear now that quantum chaos, microwave, optical, and
acoustic chaos have much in common (see [17] for a review). The common
practice is to construct an analogue for a given wave object in a
semiclassical (ray) approximation and study its chaotic properties (if any) by
well-known methods of dynamical system theory. Then, it is necessary to solve
the corresponding linear wave equation in order to find manifestations of
classical chaos in the wave-field evolution in the same range of the control
parameters. If one succeeds in that a quantum-classical or wave-ray
correspondence is announced to be established.

In this paper we perform this program with cold two-level atoms in a one-
dimensional standing-wave laser field and show that coherent dynamics of the
atomic matter waves is really complicated in that range of the control
parameters where the corresponding classical point-like atomic motion can be
strictly characterized as a chaotic one. The effect is explained by a
proliferation of atomic wave packets at the nodes of the standing wave.

##  2 Order and chaos in dynamics of atomic wave packets in a laser standing
wave

The Hamiltonian of a two-level atom, moving in a one-dimensional classical
standing-wave laser field, can be written in the frame rotating with the laser
frequency $\omega_{f}$ as follows:

$\hat{H}=\frac{\hat{P}^{2}}{2m_{a}}+\frac{\hbar}{2}(\omega_{a}-\omega_{f})\hat{\sigma}_{z}-\hbar\Omega\left(\hat{\sigma}_{-}+\hat{\sigma}_{+}\right)\cos{k_{f}\hat{X}},$
(1)

where $\hat{\sigma}_{\pm,z}$ are the Pauli operators for the internal atomic
degrees of freedom, $\hat{X}$ and $\hat{P}$ are the atomic position and
momentum operators, $\omega_{a}$ and $\Omega$ are the atomic transition and
Rabi frequencies, respectively. We will work in the momentum representation
and expand the state vector as follows:

${|\Psi(t)\closeket}=\int[a(P,t){|2\closeket}+b(P,t){|1\closeket}]{|P\closeket}dP,$
(2)

where $a(P,t)$ and $b(P,t)$ are the probability amplitudes to find atom at
time $t$ with the momentum $P$ in the excited and ground states, respectively.
The Schrödinger equation for these amplitudes is

$\displaystyle i\dot{a}(p)$
$\displaystyle=\frac{1}{2}(\omega_{r}p^{2}-\Delta)a(p)-\frac{1}{2}[b(p+1)+b(p-1)],$
(3) $\displaystyle i\dot{b}(p)$
$\displaystyle=\frac{1}{2}(\omega_{r}p^{2}+\Delta)b(p)-\frac{1}{2}[a(p+1)+a(p-1)],$

where dot denotes differentiation with respect to dimensionless time
$\tau\equiv\Omega t$, and the atomic momentum $p$ is measured in units of the
photon momentum $\hbar k_{f}$. The normalized recoil frequency,
$\omega_{r}\equiv\hbar k_{f}^{2}/m_{a}\Omega$, and the atom-field detuning,
$\Delta\equiv(\omega_{f}-\omega_{a})/\Omega$, are the control parameters.

We will treat the wave-packet motion in the dressed-state basis [18]

${|+\closeket}_{\Delta}=\sin{\Theta}{|2\closeket}+\cos{\Theta}{|1\closeket},\
{|-\closeket}_{\Delta}=\cos{\Theta}{|2\closeket}-\sin{\Theta}{|1\closeket}.$
(4)

The dressed states are eigenstates of an atom at rest in a laser field with
the eigenvalues of the quasienergy $E_{\Delta}^{(\pm)}$ and the mixing angle
$\Theta$:

$E_{\Delta}^{(\pm)}=\pm\sqrt{\frac{\Delta^{2}}{4}+\cos^{2}{x}},\
\tan{\Theta}\equiv\frac{\Delta}{2\cos{x}}-\sqrt{\left(\frac{\Delta}{2\cos{x}}\right)^{2}+1}.$
(5)

The depth of the nonresonant optical potential is

$U_{\Delta}=\left|\sqrt{\frac{\Delta^{2}}{4}+1}-\frac{|\Delta|}{2}\right|.$
(6)

At $\Delta=0$, the depth changes abruptly from $1$ to $2$. The resonant
$E_{0}^{(\pm)}$ ($\Delta=0$) and non-resonant $E_{\Delta}^{(\pm)}$
($\Delta\not=0$) potentials of an atom in a standing-wave field are drawn in
Fig. 1. The ground atomic state can be written as a superposition of the
dressed states

${|1\closeket}={|+\closeket}_{\Delta}\frac{1}{\sqrt{1+\tan^{2}{\Theta}}}-{|-\closeket}_{\Delta}\frac{\tan{\Theta}}{\sqrt{1+\tan^{2}{\Theta}}}.$
(7)

In general, an atom in the ground state, placed initially at $x_{0}=0$, will
move along two trajectories simultaneously because it is situated
simultaneously at the top of $E_{\Delta}^{(+)}$ and the bottom of
$E_{\Delta}^{(-)}$ (see Fig.1).

In the dressed-state basis, the probability amplitudes to find the atom at the
point $x$ in the potentials $E_{\Delta}^{(+)}$ and $E_{\Delta}^{(-)}$ are,
respectively,

$C_{+}(x)=a(x)\sin{\Theta}+b(x)\cos{\Theta},\
C_{-}(x)=a(x)\cos{\Theta}-b(x)\sin{\Theta},$ (8)

where the amplitudes in the bare-state basis, $a(x)$ and $b(x)$, are computed
in the position representation with the help of the Fourier transform

$a(x)={\rm
const}\int_{-\infty}^{\infty}dp^{\prime}e^{ip^{\prime}x}a(p^{\prime}),\
b(x)={\rm
const}\int_{-\infty}^{\infty}dp^{\prime}e^{ip^{\prime}x}b(p^{\prime}).$ (9)

To clarify the character of the bipotential motion we write down the
Hamiltonian of the internal degrees of freedom in the basis
${|\pm\closeket}=({|1\closeket}\pm{|2\closeket})/\sqrt{2}$:

$\hat{H}_{\rm int}=\hat{\sigma}_{z}\cos x+\frac{\Delta}{2}\hat{\sigma}_{x}.$
(10)

Let us linearize the potential in the vicinity of a node of the standing wave
and estimate a small distance the atom makes when crossing the node as
follows: $\delta x\simeq\omega_{r}|p_{\rm node}|\tau$, where $|p_{\rm node}|$
is the mean atomic momentum near the node. The quantity $\omega_{r}|p_{\rm
node}|$ is a normalized Doppler shift for an atom moving with the momentum
$|p_{\rm node}|$, i.e., $\omega_{D}\equiv\omega_{r}|p_{\rm node}|\equiv
k_{f}|v_{\rm node}|/\Omega$. We arrive now at the famous Landau–Zener problem
[19] to find a probability of transition between the states (4) when the
energy difference varies linearly in time. In other words it is a a
probability of transition from one potential to another or from one trajectory
of motion to another. The asymptotic solution is

$P_{\rm LZ}=\exp\left(-\frac{\pi\Delta^{2}}{\omega_{D}}\right).$ (11)

There are three possibilities.

1\. $\Delta^{2}\gg\omega_{D}$. The transition probability is exponentially
small and an atom moves adiabatically along, in general, two trajectories
without any transitions between them.

2\. $\Delta^{2}\ll\omega_{D}$. The probability of a Landau-Zener transition is
close to unity, and an atom changes the potential each time upon crossing any
node, i.e., it moves mainly in the resonant potentials $E_{0}^{(\pm)}$.

3\. $\Delta^{2}\simeq\omega_{D}$. The probabilities to change the potential or
to remain in the same one upon crossing a node are of the same order. In this
regime one may expect strong complexification of the wave function.

At large and small detunings, the translational motion splits into two
independent motions in the potentials $E_{\Delta}^{(\pm)}$, and the wave-
packet motion is regular in the first two cases. In the third case, the motion
is complex because of a proliferation of wave packets at the nodes of the
standing wave. We call such a motion as ”a chaotic” one by the reasons which
will be clear in section 3. The switch between the regular and chaotic regimes
of atomic motion can be easily performed by changing the detuning.

Our normalization enables to change dimensionless values of the recoil
frequency $\omega_{r}$ and the detuning $\Delta$ varying the Rabi frequency
$\Omega$. Working with a cesium atom at the transition $6S_{1/2}$ – $6P_{3/2}$
($m_{a}=133$ a. u., $\lambda_{a}=852.1$ nm and $\nu_{\rm rec}\simeq 2$ KHz
[20]), we have $\omega_{r}=10^{-3}$ at $\Omega=1$ MHz. Let the atom is
initially prepared in the ground state as a Gaussian wave packet in the
momentum space with $p_{0}=55$. The probability to find the atom with the
momentum $p$ at time $\tau$ is ${\it
P}(p,\tau)=|a(p,\tau)|^{2}+|b(p,\tau)|^{2}$.

To illustrate the difference between the regular and chaotic regimes of the
wave-packet motion we take two values of the detuning, $\Delta=1$ and
$\Delta=0.2$. At $\Delta=1$, the motion is expected to be adiabatic and
regular because $\Delta^{2}\gg\omega_{D}=0.055$, and the Landau–Zener
probability, $P_{\rm LZ}$, is exponentially small. At $\Delta=0.2$, one
expects a much more complicated wave-packet motion with nonadiabatic
transitions between the two potentials at the field nodes because
$\Delta^{2}=0.04\simeq\omega_{D}$. Figure 2 shows the dependence of the mean
atomic momentum $<p>$ over a large time scale in those two cases. In both the
cases $<p>$ oscillates in a rather irregular way. The difference is that for
an adiabatically moving wave packet, which we refer as a regular motion, the
mean atomic momentum oscillates in a narrow range around $p_{0}=55$ (the upper
curve in the figure). Whereas the range of its oscillations for a wave packet
moving with nonadiabatic transitions at the field nodes, which we refer as a
chaotic motion, is much more broad (the lower curve in the figure). It is a
simple illustration of the two different regimes of the wave-packet
propagation in terms of the classical variable.

In Fig. 3a we plot the dependence of the momentum probability-density on time
at $\Delta=1$. The initial wave packet splits from the beginning to a few
components because the initial ground state is a superposition of the dressed
states (4). The initial kinetic energy is enough to perform a ballistic
motion. The momentum changes in a comparatively small range, from 40 to 70 of
the photon-momentum units. The packet does not split at the nodes of the
standing wave but it, on the contrary, recollects in the momentum space at the
nodes and spreads in between. However, this recollection smears out in course
of time.

At $\Delta=0.2$, the atomic ground state is the following superposition of the
dressed states: ${|1\closeket}\simeq
0.74{|+\closeket}_{\Delta}+0.66{|-\closeket}_{\Delta}$. The
${|+\closeket}_{\Delta}$-component of the initial wave packet, i.e., that one,
starting from the top of the potential $E_{\Delta}^{(+)}$, overcomes the
barriers of that potential and moves in the positive direction of the axis $x$
proliferating at the nodes. As to the ${|-\closeket}$-component with decreased
values of $p$, it will be trapped in the potential well performing
oscillations in the momentum and position spaces. The period of those
oscillations is about $T\simeq 280$ which is equal approximately to the period
of revivals of the Rabi oscillations for the population inversion. Figure 3b
illustrates the effect of simultaneous trapping and ballistic motion of the
atomic wave packet in the chaotic regime resulting in a broad momentum
distribution, from $p=-60$ to $p\simeq 80$.

To illustrate the nonadiabatic transitions from one potential to another and
their absence at the nodes more explicitly, we go to the position space and
compute the probabilities $|C_{\pm}(x,\tau)|^{2}$ (8) to be at the point $x$
at time $\tau$ in the potentials $E_{\Delta}^{(+)}$ and $E_{\Delta}^{(-)}$,
respectively. In Fig. 4 we plot the evolution of the probability density
$|C_{-}(x,\tau)|^{2}$ in the frame moving with the initial atomic velocity
$\omega_{r}p_{0}=0.055$. The slope straight lines in the figure mark positions
of the nodes in the moving frame. In the case of the regular motion at
$\Delta=1$ (Fig. 4a), no transitions happen when the atom crosses the nodes.
In the chaotic regime at $\Delta=0.2$, one observes visible changes in the
probability-density $|C_{-}(x)|^{2}$ exactly at the node lines (see Fig. 4b).
It means transitions from one trajectory to another at the field nodes that
should occur in a specific range of the control parameters if
$\Delta^{2}\simeq\omega_{D}$. This results in a proliferation of the
components of the wave packet at the nodes and, therefore, a complexification
of the wave function (see Fig. 3b).

The Wigner function can be used to visualize complexity of the wave function
in the chaotic regime of the atomic motion. We compute the evolution of the
Wigner function of the ground state in the momentum space

$W_{b}(p,x,\tau)={\rm
const}\int_{-\infty}^{\infty}dp^{\prime}e^{-ip^{\prime}x}a(p-p^{\prime}/2)a^{*}(p+p^{\prime}/2).$
(12)

This quantity gives a quasi-probability distribution corresponding to a
general quantum state (2). Figure 5 shows a contour plot of the Wigner
function (12) at two moments of time, $\tau=50$ and $\tau=200$, when the atom
moves in a regular way ($\Delta=1$). Figure 6 is a contour plot of this
function at the same times, but with an atom making nonadiabatic transitions
at the field nodes ($\Delta=0.2$). In the chaotic regime (Fig. 6b) we see a
dust-like distribution of nonzero values of the Wigner function at $\tau=200$
which occupy much more larger area in the phase space than the function for
the regular motion (Fig. 5b).

##  3 Quantum-classical correspondence

In this section we compare the results of the quantum treatment with those
obtained for the same problem in the semiclassical approximation when the
translational motion has been treated as a classical one [14, 15]. We must
compare quantum results for a single atomic wave packet not with a single
point-like atom but with an ensemble of point-like atoms. Dynamical chaos has
been found and analyzed in detail in Refs. [14, 15] in the semiclassical
approximation. Both the internal and external degrees of freedom of a two-
level atom in a standing wave field have been shown to be chaotic in a
specific range of values of the detuning $\Delta$, the recoil frequency
$\omega_{r}$, and the initial momentum $p_{0}$.

Coherent semiclassical evolution of a point-like two-level atom is governed by
the Hamilton-Schrödinger equations with the same normalization as in the
quantum case [15]

$\begin{gathered}\dot{x}=\omega_{r}p,\quad\dot{p}=-u\sin x,\quad\dot{u}=\Delta
v,\\\ \dot{v}=-\Delta u+2z\cos x,\quad\dot{z}=-2v\cos x,\end{gathered}$ (13)

where

$u\equiv 2\operatorname{Re}(a_{0}b_{0}^{*}),\
v\equiv-2\operatorname{Im}(a_{0}b_{0}^{*}),\ z\equiv|a_{0}|^{2}-|b_{0}|^{2}$
(14)

are the atomic-dipole components ($u$ and $v$) and population-inversion ($z$),
and $a_{0}$ and $b_{0}$ are the complex-valued probability amplitudes to find
the atom in the excited, ${|2\closeket}$, and ground, ${|1\closeket}$, states,
respectively. The system (13) has two integrals of motion, the total energy

$W\equiv\frac{\omega_{r}}{2}p^{2}-u\cos x-\frac{\Delta}{2}z,$ (15)

and the length of the Bloch vector, $u^{2}+v^{2}+z^{2}=1$. Equations (13)
constitute a nonlinear Hamiltonian autonomous system with two and half degrees
of freedom and two integrals of motion.

In classical mechanics there is a qualitative criterion of dynamical chaos,
the maximal Lyapunov exponent $\lambda$, which measures the mean rate of the
exponential divergence of initially closed trajectories in the phase space. In
Fig. 7 we plot this quantity, computed with semiclassical equations of motion
(13) vs the detuning $\Delta$ at $\omega_{r}=10^{-3}$ and $p_{0}=55$. At zero
detuning, the set of semiclassical equations acquires an additional integral
of motion and becomes integrable. The center-of-mass motion is regular at
small values of the detuning, $\Delta\ll 1$, and at large ones, $\Delta>0.8$.
Positive values of $\lambda$ at $0<\Delta<0.8$ characterize unstable motion.
Local instability of the center-of-mass motion produce chaotic motion of an
atom in a rigid standing wave without any modulation of its parameters in
difference from the situation with atoms in a periodically kicked optical
lattice [7, 8]. There is a range of initial conditions and the control
parameters where the center-of-mass motion in an absolutely deterministic
standing wave resembles a random walking. It means that a point-like atom
alternates between flying through the lattice, and being trapped in its wells
changing the direction of motion in a random-like way [15, 16].

In Fig. 8 we illustrate semiclassical chaos with the Poincaré mapping for a
number of atomic ballistic trajectories in the western ($u<0$) and eastern
($u>0$) hemispheres of the Bloch sphere $(u,v,z)$ on the plane $v-z$. One can
see a typical structure with regions of regular motion in the form of islands
and chains of islands filled by regular trajectories. The islands are imbedded
into a stochastic sea, and they are produced by nonlinear resonances of
different orders. Increasing the resolution of the mapping, one can see that
large islands are surrounded by islands of a smaller size each of which, in
turn, is surrounded by a chain of even more smaller islands, and so on.

In quantum mechanics there is no well-defined notion of a trajectory in the
phase space, the very phase space is not continuous due to the Heisenberg
uncertainty relation, and, hence, the Lyapunov exponents can not be computed
(however, see Ref. [21] where a notion of a Lyapunov exponent for quantum
dynamics has been discussed). The main result of this paper is the
establishment of the fact, that chaotically-like complexification of the wave
function, caused by nonadiabatic transitions at the field nodes, occurs
exactly at the same range of the control parameters where the semiclassical
dynamics has been shown to be chaotic in Refs. [14, 15]. It should be stressed
that quantum motion of a wave packet with nonadiabatic transitions between the
two optical potentials is compared with the center-of-mass motion of an
ensemble of atoms each of which moves in a single optical potentials. So, when
we say about a quantum-classical correspondence we mean a correspondence
between the wave function of a single quantum atom and the trajectories of the
ensemble of classical atoms with different values of the initial momentum
$p_{0}$ and at the other equal conditions.

##  4 Conclusion

We have studied coherent dynamics of cold atomic wave packets in a one-
dimensional standing-wave laser field. The problem has been considered in the
momentum representation and in the dressed-state basis where the motion of a
two-level atom was interpreted as a propagation in two optical potentials. The
character of that motion has been shown to depend strongly on the ratio of the
squared detuning, $\Delta^{2}$, to the normalized Doppler shift $\omega_{D}$.
In the regular regime, when $\Delta$ is comparatively large or small, wave
packets move in a simple way. The chaotic regime occurs if
$\Delta^{2}\simeq\omega_{D}$ when the probability for an atom to make
nonadiabatic transitions while crossing the nodes of the standing wave is
large. Atom in this regime of motion simultaneously moves ballistically and is
trapped in a well of the optical potential. This type of motion and
proliferation of wave packets at the nodes result in a complexification of the
wave function both in the momentum and position spaces manifesting itself in
the irregular behavior of the Wigner function.

Comparing the results of the quantum treatment with those obtained in the
semiclassical approximation, when the translational motion has been treated as
a classical one [14, 15], we have found that the wave-packet dynamics is
complicated exactly in that range of the atom-field detuning and recoil
frequency where the classical center-of-mass motion has been shown to be
chaotic in the sense of exponential sensitivity to small variations in initial
conditions or parameters.

As to possible practical applications of the results obtained we mention
atomic lithography to produce small-scale complex prints of cold atoms (see,
for example, beautiful experiments on coherent matter-wave manipulation [22,
23, 24]), new ways to manipulate atomic motion in optical lattices by varying
the atom-filed detuning and atomic ratchets with cold atoms.

## Acknowledgments

This work was supported by the Russian Foundation for Basic Research (project
no. 09-02-00358), the Integration grant from the Far-Eastern and Siberian
branches of the RAS, and the Program “Fundamental Problems of Nonlinear
Dynamics” of the RAS. I would like to thank L. Konkov for preparing Figs.6 and
7.

## References

  * [1] V.G. Minogin, V.S. Letokhov, Laser Light Pressure on Atoms, Gordon and Breach, New York, 1987.
  * [2] A.P. Kazantsev, G.I. Surdutovich, V.P. Yakovlev, Mechanical Action of Light on Atoms, Singapore, World Scientific, 1990.
  * [3] P. Meystre, Atom Optics, New York, Springer-Verlag, 2001\.
  * [4] W.P. Schleich, Quantum Optics in Phase Space, New York, Wiley, 2001.
  * [5] M.V. Fedorov, M.A. Efremov, V.P. Yakovlev, W.P. Schleich, JETP 97 (2003) 522 [Zh. Eksp. Teor. Fiz. 124 (2003) 578].
  * [6] R. Graham, M. Schlautmann, P. Zoller, Phys. Rev. A 45 (1992) R19 .
  * [7] M.G. Raizen, Adv. At. Mol. Opt. Phys. 41 (1999) 43. D.A. Steck, et al, Science 293 (2001) 274.
  * [8] W.K. Hensinger, N.R. Heckenberg, G.J. Milburn, H. Rubinsztein-Dunlop, J. Opt. B: Quantum Semiclass. Opt. 5 (2003) 83.
  * [9] M. Sadgrove, S. Wimberger, S. Parkins, and R. Leonhardt, Phys. Rev. Lett. 94, (2005) 174103.
  * [10] C. Mennerat-Robilliard et al, Phys. Rev. Lett. 82, (1999) 851.
  * [11] M. Schiavoni, L. Sanchez-Palencia, F. Renzoni, G. Grynberg, Phys. Rev. Lett. 90, (2003) 094101.
  * [12] G. G. Carlo, G. Benenti, G. Casati, S. Wimberger, O. Morsch, R. Manella, E. Arimondo, Phys. Rev. A 74, (2006) 033617.
  * [13] L. Sanchez-Palencia, F. R. Carminati, M. Schiavoni, F. Renzoni, G. Grynberg, Phys. Rev. Lett. 88, (2002) 133903.
  * [14] S.V. Prants, JETP Letters 75 (2002) 651. [Pis’ma ZhETF 75, 777 (2002)].
  * [15] V.Yu. Argonov, S.V. Prants, Phys. Rev. A 75 (2007) art. 063428.
  * [16] V.Yu. Argonov, S.V. Prants, Phys. Rev. A 78 (2008) art. 043413.
  * [17] D. Makarov, S. Prants, A. Virovlyansky, G. Zaslavsky, Ray and wave chaos in ocean acoustics, Singapore, World Scientific, 2010.
  * [18] C. Cohen-Tannoudji, J. Dupon-Roc, G. Grynberg, Atom-Photon Interaction, Weinheim, Wiley, 1998.
  * [19] L. Landau, Phys. Z. Sowjetunion 2 (1932) 46. C. Zener, Proc. R. Soc. London A 2 (1932) 137.
  * [20] C.S. Adams, M. Sigel, J. Mlynek, Phys. Rep. 240 (1994) 143.
  * [21] V.I. Man’ko, R.V. Mendes, Physica D 145 330 (2000) 330.
  * [22] G. Zabow, R.S. Conroy, M. G. Prentiss, Phys. Rev. Lett. 92, (2004) 180404.
  * [23] A. Zenesini, H. Lignier, D. Ciampini, O. Morsch, E. Arimondo, Phys. Rev. Lett. 102, (2009) 100403.
  * [24] S. Wu, A. Tonyushkin, M. G. Prentiss, Phys. Rev. Lett. 103, (2009) 034101.

Figure 1: The resonant $E_{0}^{(\pm)}$ ($\Delta=0$, dotted curves) and non-
resonant $E_{\Delta}^{(\pm)}$ ($\Delta\not=0$, solid curves) potentials of a
two-level atom in a standing-wave laser field. An atomic wave packet, centered
at $x_{0}=0$, and its initial evolution in the upper and lower potentials are
shown schematically. Figure 2: Mean atomic momentum $\left<p\right>$ vs time
in the regular ($\Delta=1$, the upper bold curve) and chaotic ($\Delta=0.2$)
regimes of motion.

Figure 3: Momentum probability-density distribution ${\it P}(p,\tau)$ vs time
in (a) the regular and (b) chaotic regimes of motion. The color codes the
corresponding values of ${\it P}(p)$.

Figure 4: The probability $|C_{-}(x)|^{2}$ to find the atom in the potential
$E_{\Delta}^{(-)}$ in the moving frame of reference in (a) the regular and (b)
chaotic regimes of motion. The slopes mark positions of the nodes in the
moving frame and the color codes the values of $|C_{-}(x)|^{2}$.

Figure 5: Contour plots of the Wigner function at (a) $\tau=50$ and (b)
$\tau=200$ for the regular regime of the atomic motion. Color online: red and
blue areas show positive and negative values of the Wigner function,
respectively.

Figure 6: The same as in Fig. 5 for the chaotic regime of motion. Figure 7:
Maximal Lyapunov exponent $\lambda$, computed with semiclassical equations of
motion (13), vs the detuning $\Delta$ at $\omega_{r}=10^{-3}$ and $p_{0}=55$.
Figure 8: Poincaré sections of the Bloch sphere illustrating the effect of
semiclassical chaos with point-like atoms at $\omega_{r}=10^{-5}$,
$\Delta=-0.05$ and the total atomic energy $W=36.45$. (a) $u<0$ (western Bloch
hemisphere), (b) $u>0$ (eastern Bloch hemisphere).


# Blind Omnidirectional Image Quality Assessment: Integrating Local Statistics
and Global Semantics

###### Abstract

Omnidirectional image quality assessment (OIQA) aims to predict the perceptual
quality of omnidirectional images that cover the whole 180$\times$360∘ viewing
range of the visual environment. Here we propose a blind/no-reference OIQA
method named S2 that bridges the gap between low-level statistics and high-
level semantics of omnidirectional images. Specifically, statistic and
semantic features are extracted in separate paths from multiple local
viewports and the hallucinated global omnidirectional image, respectively. A
quality regression along with a weighting process is then followed that maps
the extracted quality-aware features to a perceptual quality prediction.
Experimental results demonstrate that the proposed S2 method offers highly
competitive performance against state-of-the-art methods.

Index Terms—  Omnidirectional image, blind image quality assessment, low-level
statistics, high-level semantics

##  1 Introduction

The rapid recent advancement in virtual reality (VR) technologies makes it
possible to create immersive multimedia quality-of-experience (QoE) for end-
users. As a representative form of VR, omnidirectional content has
increasingly emerged in our daily life. To evaluate and optimize the
perceptual QoE of omnidirectional content, objective omnidirectional image
quality assessment (OIQA) models play a critical roles in the development of
modern VR systems.

In the literature, objective OIQA models have emerged that follow both full-
reference (FR) and no-reference (NR) frameworks. FR-OIQA models assume full
access to information of the reference image and are usually direct extensions
of traditional FR methods developed for regular rectangular 2D image quality
assessment (IQA). For example, based upon the peak signal-to-noise ratio
(PSNR), Yu et al. [1] propose the spherical PSNR (S-PSNR) algorithm, where
PSNR is calculated for uniformly distributed points on a sphere instead of
projected rectangular image. In [2], the weighted-to-spherically uniform PSNR
(WS-PSNR) method is presented, where a weighting map is created by considering
the stretched degree. Zakharchenko et al. [3] propose the Craster parabolic
projection PSNR (CPP-PSNR) approach, which maps the reference and distorted
omnidirectional images on the Craster parabolic projection followed by PSNR
computation.

Fig. 1: Perceptual cues in omnidirectional image quality assessment. Existing
models extract spatial information from various viewports and may obtain help
from global projected maps, whereas the proposed method combines local image
statistics and global semantic reconstruction.

Fig. 2: Framework of the proposed S2 method for blind OIQA.

NR-OIQA methods do not require access to the reference image and are more
desirable in many application scenarios. Existing NR-OIQA approaches can
generally be classified into two categories, depending on whether the
conventional hand-crafted or learned deep features are employed for quality
prediction. Multi-frequency information and local-global naturalness are
applied to develop the MFILGN model [4]. More recent models employ deep
convolutional neural networks (CNNs) or graph convolution networks (GCNs).
These models demonstrate promising performance, including the multi-channel
CNN for blind 360-degree image quality assessment (MC360IQA) [5], the viewport
oriented graph convolution network (VGCN) [6], and its variant named adaptive
hypergraph convolutional network (AHGCN) [7].

In a 360-degree viewing environment, e.g. using a head-mounted device, the
observer is not able to visualize the whole omnidirectional content
simultaneously, and thus an important step in the human subjective viewing
experience is to establish or reconstruct a sense of the global semantics by
browsing and integrating information from many viewports. During the course of
image quality assessment, such global semantics are integrated with local
observations on image fidelity, naturalness, and/or artifacts to produce an
overall quality evaluation. Motivated by this observation, we propose a
statistic and semantic oriented quality prediction framework named S2 for
blind OIQA as illustrated in Fig. 1, by integrating features extracted from
both low-level image statistics of multiple local viewports and high-level
semantics of the hallucinated global omnidirectional image. A quality
regression module is then leveraged to map the collection of the quality-
sensitive features extracted from the two separate paths to an overall
prediction of the subjective quality rating. Extensive experimental results
demonstrate that the proposed method is superior to many state-of-the-art
quality assessment models. In addition, we make some interesting observations
on the relationship between semantic confidence and image distortions, as well
as how the individual components affect the ultimate quality prediction
performance in ablation studies.

##  2 Proposed Method

The overall framework of the proposed S2 method is shown in Fig. 2, which
consists of a statistic path, a semantic path, and a final quality regression
step.

Since a variety of viewports are browsed by the viewers, we first convert the
distorted omnidirectional image (OI) to multiple viewports. Given each input
distorted OI denoted by $D$, we exploit the non-uniform viewport sampling
strategy [8, 9] and obtain $N$ viewports $V_{n},n=1,2,...,N$.

To capture the multi-scale characteristics of the human visual system [10], we
construct pyramid representations [11, 12] of multiple local viewports.
Specifically, multi-level Laplacian pyramids [13] are created by iterative
Gaussian filtering, down-sampling, and subtracting, resulting in Gaussian and
Laplacian pyramids in the same process. For a specific viewport $V_{n}$,
layers of the Gaussian pyramid are calculated as follows:

$\small G_{n}^{i}(x,y)=\left\\{\begin{array}[]{l}V_{n},i=1\\\
\sum\limits_{u=-2}^{2}{\sum\limits_{v=-2}^{2}k(u,v)G_{n}^{i-1}(2x+u,2y+v)},i>1\end{array}\right.,$
(1)

where $i$ is the layer index of the Gaussian pyramid, $x\in[0,X)$ and
$y\in[0,Y)$ are the pixel position indices in which $X$ and $Y$ are the image
dimensions, and $k(u,v)$ denotes the generating kernel that is typically
defined by the coefficients of a low pass filter such as a 2D Gaussian filter.

We then interpolate each layer of the Gaussian pyramid by:

$\small\hat{G}_{n}^{i}(x,y)=4\sum_{u=-2}^{2}\sum_{v=-2}^{2}k(u,v)G_{n}^{i}\left(\frac{u+x}{2},\frac{v+y}{2}\right).$
(2)

The residual between the current layer of the Gaussian pyramid and the
interpolation result from the next layer defines the current layer of the
Laplacian pyramid:

$\small L_{n}^{i}=G_{n}^{i}-\hat{G}_{n}^{i+1}.$ (3)

Since the computation of the $i$-th layer in the Laplacian pyramid requires
the $(i+1)$-th layer of the Gaussian pyramid, the number of layers in the
Laplacian pyramid is one less than that in the Gaussian pyramid.

To extract features from the Gaussian pyramid, we compute the default uniform
local binary pattern (LBP) descriptors, resulting in 59 statistics for each
Gaussian layer. When a 3-layer Gaussian pyramid is employed, this leads to 177
Gaussian pyramid features denoted by $f_{GP}$. For a Laplacian pyramid,
motivated by the success of natural scene statistics (NSS) in IQA research
[14, 15, 16], we extract mean subtracted and contrast normalized coefficients,
leading to 36 features for each layer. When a 2-layer Laplacian pyramid is
employed, this results in 72 Laplacian pyramid features denoted by $f_{LP}$.
The full statistic feature set $f_{st}$, one for each viewport, is obtained by
concatenating the statistical features extracted from the Gaussian and
Laplacian pyramids as:

$\small f_{st}=\left[f_{GP},f_{LP}\right].$ (4)

We employ the VGGNet trained on the large ImageNet dataset [17] as the
semantic feature extraction backbone, mainly for its simplicity and ability to
capture image distortion-related representations [18]. In [19], three
different structures of VGGNet have been proposed to balance between
complexity and accuracy, namely fast VGG (VGG-F), medium VGG (VGG-M), and slow
VGG (VGG-S). Each of them contains 5 convolutional (Conv) layers and 3 fully
connected (FC) layers. The first two FC layers have 4,096 neurons, while the
last one has 1,000 nodes indicating the 1,000 classes for image recognition.
In our current implementation, we select the deep features from the first FC
layer of VGG-M as our semantic feature set $f_{se}$:

$\small f_{se}=FC_{1}(D).$ (5)

To learn the mapping from features to quality labels, we feed the statistic
features and semantic features separately to support vector regression (SVR)
models [20], and denote the regressed statistic and semantic quality scores as
$Q_{st}$ and $Q_{se}$, respectively. The overall quality score is calculated
by a weighted average:

$\small Q_{\text{overall }}=wQ_{st}+(1-w)Q_{\text{se }},$ (6)

where $w$ is a weighting factor that determines the relative importance of the
statistic and semantic feature predictors.

##  3 Validation

###  3.1 Experimental Setup and Performance Comparison

We evaluate the proposed approach on the CVIQD subjective database [21], which
is so far a relatively large and widely adopted database containing both
omnidirectional images and their corresponding quality labels given by human
subjects. It consists of 16 original images and 528 distorted images produced
by three classic image or video coding technologies, namely JPEG, AVC, and
HEVC. The subjective quality ratings in the form of mean opinion score (MOS)
are rescaled to the range of [0, 100], for which a higher MOS represents
better perceptual image quality.

To compare the performance of various IQA models, we take Spearman Rank-Order
Correlation Coefficient (SROCC), Pearson Linear Correlation Coefficient (PLCC)
and Root Mean Squared Error (RMSE) as the evaluation criteria. Before
calculating the PLCC and RMSE, a 5-parameter logistic nonlinear fitting
approach [22] is implemented to map the predicted quality into the subjective
quality space.

The database is randomly divided into 80% data for training and the remaining
20% data for cross-validation. In order to relieve the uncertainty in
training/testing splitting, we repeat this random-splitting and cross-
validation process 100 times and report the median performance.

Table 1: Performance comparisons of objective models.

Types | Methods | SROCC | PLCC | RMSE
---|---|---|---|---
FR-IQA | PSNR | 0.6239 | 0.7008 | 9.9599
SSIM [23] | 0.8842 | 0.9002 | 6.0793
MS-SSIM [10] | 0.8222 | 0.8521 | 7.3072
FSIM [24] | 0.9152 | 0.9340 | 4.9864
DeepQA [25] | 0.9292 | 0.9375 | 4.8574
FR-OIQA | S-PSNR [1] | 0.6449 | 0.7083 | 9.8564
WS-PSNR [2] | 0.6107 | 0.6729 | 10.3283
CPP-PSNR [3] | 0.6265 | 0.6871 | 10.1448
NR-IQA | BRISQUE [26] | 0.8180 | 0.8376 | 7.6271
BMPRI [27] | 0.7470 | 0.7919 | 8.5258
DB-CNN [28] | 0.9308 | 0.9356 | 4.9311
NR-OIQA | MFILGN [4] | 0.9670 | 0.9751 | 3.1036
MC360IQA [5] | 0.9428 | 0.9429 | 4.6506
VGCN [6] | 0.9639 | 0.9651 | 3.6573
AHGCN [7] | 0.9623 | 0.9643 | 3.6990
Proposed S2 | 0.9710 | 0.9781 | 2.8945

The performance of the proposed algorithm is compared against state-of-the-art
quality assessment models, including five FR-IQA, three FR-OIQA, three NR-IQA,
and four NR-OIQA methods. The results are shown in TABLE 1, where we observe
that for FR-IQA metrics, the PSNR-based models are inferior to more advanced
approaches such as structural (SSIM, MS-SSIM, FSIM) and deep learning (DeepQA)
models. Somewhat surprisingly, the FR-OIQA methods do not help further improve
upon FR-IQA approaches. By contrast, the NR-OIQA models show significant
superiority over NR-IQA methods. This is likely due to their specific design
to capture the characteristics of omnidirectional images. Among all metrics
tested, the proposed S2 method demonstrates highly competitive performance.

Fig. 3: The relationship between semantic confidence and image distortions.
The first, second and third rows correspond to three distortion types (JPEG,
AVC and HEVC compression, respectively). The first, second and third columns
correspond to increasing distortion levels (low, medium and high,
respectively) for each distortion type.

Fig. 4: Performance results of ablation experiments.

###  3.2 Semantic Confidence Versus Image Distortion

Since the proposed method contains a semantic path, it is interesting to
observe the relationship between semantic confidence and image distortion. An
example of distorted omnidirectional images with different JPEG, AVC and HEVC
distortion levels is shown in Fig. 3, where from the first column to the third
column, we observe that as the degree of distortion increases, the semantic
confidence level decreases. This suggests that semantic information may be
highly related to perceptual image quality. It is also interesting to see that
the semantic confidence shows various sensitivities to different distortion
types. In particular, the drop in semantic confidence levels is much less in
more advanced image/video coding method HEVC than in earlier JPEG and AVC
encoders.

###  3.3 Ablation and Parameter Sensitivity Tests

We evaluate the contributions from the statistic and semantic paths by
ablation experiments, and the results are shown in Fig. 4, where GP1, GP2 and
GP3, respectively, represent the cases of using the first, second, and third
layers of Gaussian pyramid statistics only. GP denotes the case of using three
layers of Gaussian pyramid statistics. We find that the performance increases
gradually. Similarly, LP1 and LP2, respectively, correspond to the cases of
using the first and second layers of Laplacian pyramid statistics only, while
LP denotes the case of using 2-layer Laplacian pyramid statistics. The results
show that LP produces the best performance among the three. The cases of
adopting the statistic path and the semantic path only are denoted by St and
Se, respectively. It is observed that either path alone can achieve promising
quality prediction performance, but adopting both paths (i.e. the All case)
delivers the best performance. Relative speaking, the more dominant factor
seems to be the statistic path. This may not be surprising as the statistic
features come from different viewports directly visualized by human subject
while the global semantics offer complementary information for additional cues
in quality assessment.

Table 2: Performance comparisons for different viewport numbers in the
statistic path.

Numbers | SROCC | PLCC | RMSE
---|---|---|---
$6$ | 0.9684 | 0.9769 | 3.0083
$20$ | 0.9686 | 0.9777 | 2.9626
$80$ | 0.9683 | 0.9771 | 2.9501

Table 3: Performance comparisons for different neural network architectures in
the semantic path.

Architectures | SROCC | PLCC | RMSE
---|---|---|---
VGG-F | 0.9497 | 0.9537 | 4.2107
VGG-M | 0.9517 | 0.9576 | 4.0329
VGG-S | 0.9451 | 0.9486 | 4.4345

Because different parameter settings may be employed in the implementations of
the proposed framework, here we test the sensitivity of our model with regard
to various viewport numbers and semantic architectures. The results are
reported in TABLE 2 and TABLE 3, respectively. We can see that the proposed
model is insensitive to the viewport number. This allows us to reduce the
number of viewports (for example, 6) to alleviate the computational complexity
in real-world applications. The results also show that VGG-M outperforms the
other neural network architectures in the semantic path. The possible reason
may be that VGG-M achieves a preferable tradeoff between algorithm complexity
and accuracy, making it a desired option for deep semantic backbone.

##  4 Conclusion

We propose a novel S2 framework for blind omnidirectional image quality
assessment that integrates both local low-level statistic and global high-
level semantic features. Extensive experiments show that the proposed method
achieves state-of-the-art performance. Observations on the relationship
between semantic confidence and image distortion, and the ablation/sensitivity
tests offer additonal useful insights. Under the same framework, more advanced
models for statistic and semantic analysis may be employed in the future,
aiming for more accurate QoE assessment models that may help drive the
advancement of immersive multimedia systems.

## References

  * [1] Matt Yu, Haricharan Lakshman, and Bernd Girod,  “A framework to evaluate omnidirectional video coding schemes,”  in IEEE International Symposium on Mixed and Augmented Reality, 2015, pp. 31–36.
  * [2] Yule Sun, Ang Lu, and Lu Yu,  “Weighted-to-spherically-uniform quality evaluation for omnidirectional video,”  IEEE Signal Processing Letters, vol. 24, no. 9, pp. 1408–1412, 2017\.
  * [3] Vladyslav Zakharchenko, Kwang Pyo Choi, and Jeong Hoon Park,  “Quality metric for spherical panoramic video,”  in Optics and Photonics for Information Processing X. International Society for Optics and Photonics, 2016, vol. 9970, p. 99700C.
  * [4] Wei Zhou, Jiahua Xu, Qiuping Jiang, and Zhibo Chen,  “No-reference quality assessment for 360-degree images by analysis of multifrequency information and local-global naturalness,”  IEEE Transactions on Circuits and Systems for Video Technology, 2021\.
  * [5] Wei Sun, Xiongkuo Min, Guangtao Zhai, Ke Gu, Huiyu Duan, and Siwei Ma,  “MC360IQA: A multi-channel cnn for blind 360-degree image quality assessment,”  IEEE Journal of Selected Topics in Signal Processing, vol. 14, no. 1, pp. 64–77, 2019.
  * [6] Jiahua Xu, Wei Zhou, and Zhibo Chen,  “Blind omnidirectional image quality assessment with viewport oriented graph convolutional networks,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 31, no. 5, pp. 1724–1737, 2020.
  * [7] Jun Fu, Chen Hou, Wei Zhou, Jiahua Xu, and Zhibo Chen,  “Adaptive hypergraph convolutional network for no-reference 360-degree image quality assessment,”  in Proceedings of the 30th ACM International Conference on Multimedia, 2022, pp. 961–969.
  * [8] Jiahua Xu, Ziyuan Luo, Wei Zhou, Wenyuan Zhang, and Zhibo Chen,  “Quality assessment of stereoscopic 360-degree images from multi-viewports,”  in IEEE Picture Coding Symposium, 2019, pp. 1–5.
  * [9] Zhibo Chen, Jiahua Xu, Chaoyi Lin, and Wei Zhou,  “Stereoscopic omnidirectional image quality assessment based on predictive coding theory,”  IEEE Journal of Selected Topics in Signal Processing, vol. 14, no. 1, pp. 103–117, 2020.
  * [10] Zhou Wang, Eero P Simoncelli, and Alan C Bovik,  “Multiscale structural similarity for image quality assessment,”  in The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003. IEEE, 2003, vol. 2, pp. 1398–1402.
  * [11] Zhou Wang and Eero P Simoncelli,  “Reduced-reference image quality assessment using a wavelet-domain natural image statistic model,”  in Human vision and electronic imaging X. International Society for Optics and Photonics, 2005, vol. 5666, pp. 149–159.
  * [12] Valero Laparra, Johannes Ballé, Alexander Berardino, and Eero P Simoncelli,  “Perceptual image quality assessment using a normalized laplacian pyramid,”  Electronic Imaging, vol. 2016, no. 16, pp. 1–6, 2016.
  * [13] Peter J Burt and Edward H Adelson,  “The laplacian pyramid as a compact image code,”  in Readings in Computer Vision, pp. 671–679. Elsevier, 1987.
  * [14] Anish Mittal, Rajiv Soundararajan, and Alan C Bovik,  “Making a ”completely blind“ image quality analyzer,”  IEEE Signal Processing Letters, vol. 20, no. 3, pp. 209–212, 2012\.
  * [15] Yuming Fang, Kede Ma, Zhou Wang, Weisi Lin, Zhijun Fang, and Guangtao Zhai,  “No-reference quality assessment of contrast-distorted images based on natural scene statistics,”  IEEE Signal Processing Letters, vol. 22, no. 7, pp. 838–842, 2014\.
  * [16] Zhibo Chen, Wei Zhou, and Weiping Li,  “Blind stereoscopic video quality assessment: From depth perception to overall experience,”  IEEE Transactions on Image Processing, vol. 27, no. 2, pp. 721–734, 2017.
  * [17] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei,  “ImageNet: A large-scale hierarchical image database,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 248–255.
  * [18] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang,  “The unreasonable effectiveness of deep features as a perceptual metric,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 586–595.
  * [19] Ken Chatfield, Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman,  “Return of the devil in the details: Delving deep into convolutional nets,”  arXiv preprint arXiv:1405.3531, 2014.
  * [20] Chih-Chung Chang and Chih-Jen Lin,  “LIBSVM: a library for support vector machines,”  ACM Transactions on Intelligent Systems and Technology, vol. 2, no. 3, pp. 1–27, 2011.
  * [21] Wei Sun, Ke Gu, Siwei Ma, Wenhan Zhu, Ning Liu, and Guangtao Zhai,  “A large-scale compressed 360-degree spherical image database: From subjective quality evaluation to objective model comparison,”  in IEEE 20th International Workshop on Multimedia Signal Processing, 2018, pp. 1–6.
  * [22] Video Quality Experts Group et al.,  “Final report from the video quality experts group on the validation of objective models of video quality assessment, phase II,”  VQEG, 2003.
  * [23] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli,  “Image quality assessment: from error visibility to structural similarity,”  IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.
  * [24] Lin Zhang, Lei Zhang, Xuanqin Mou, and David Zhang,  “FSIM: A feature similarity index for image quality assessment,”  IEEE Transactions on Image Processing, vol. 20, no. 8, pp. 2378–2386, 2011.
  * [25] Jongyoo Kim and Sanghoon Lee,  “Deep learning of human visual sensitivity in image quality assessment framework,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 1676–1684.
  * [26] Anish Mittal, Anush Krishna Moorthy, and Alan Conrad Bovik,  “No-reference image quality assessment in the spatial domain,”  IEEE Transactions on Image Processing, vol. 21, no. 12, pp. 4695–4708, 2012.
  * [27] Xiongkuo Min, Guangtao Zhai, Ke Gu, Yutao Liu, and Xiaokang Yang,  “Blind image quality estimation via distortion aggravation,”  IEEE Transactions on Broadcasting, vol. 64, no. 2, pp. 508–517, 2018.
  * [28] Weixia Zhang, Kede Ma, Jia Yan, Dexiang Deng, and Zhou Wang,  “Blind image quality assessment using a deep bilinear convolutional neural network,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 30, no. 1, pp. 36–47, 2018.

111© 2023 IEEE. Personal use of this material is permitted. Permission from
IEEE must be obtained for all other uses, in any current or future media,
including reprinting/republishing this material for advertising or promotional
purposes, creating new collective works, for resale or redistribution to
servers or lists, or reuse of any copyrighted component of this work in other
works.