#  Polar prediction of natural videos

Pierre-Étienne H. Fiquet  Eero P. Simoncelli

###### Abstract

Observer motion and continuous deformations of objects and surfaces imbue
natural videos with distinct temporal structures, enabling partial prediction
of future frames from past ones. Conventional methods first estimate local
motion, or optic flow, and then use it to predict future frames by warping or
copying content. Here, we explore a more direct methodology, in which each
frame is mapped into a learned representation space where the structure of
temporal evolution is more readily accessible. Motivated by the geometry of
the Fourier shift theorem and its group-theoretic generalization, we formulate
a simple architecture that represents video frames in learned local polar
coordinates. Specifically, we construct networks in which pairs of
convolutional channel coefficients are treated as complex-valued, and are
optimized to evolve with slowly varying amplitudes and linearly advancing
phases. We train these models on next-frame prediction in natural videos, and
compare their performance with that of conventional methods using optic flow
as well as predictive neural networks. We find that the polar predictor
achieves better performance while remaining interpretable and fast, thereby
demonstrating the potential of a flow-free video processing methodology that
is trained end-to-end to predict natural video content.

unsupervised representation learning, next-frame video prediction, invariance
/ equivariance factorization, harmonic analysis, computational neuroscience

##  1 Introduction

One way to frame the fundamental problem of vision is that of representing the
signal in a form that is more useful for performing visual tasks, be they
estimation, recognition, or motor action. Perhaps the most general “task” is
that of temporal prediction, which has been proposed as a fundamental goal for
unsupervised learning of visual representations (Földiák, 1991). But previous
research along these lines has generally focused on estimating temporal
transformations rather than using them to predict: for example, extracting
slow features (Wiskott & Sejnowski, 2002), or finding sparse codes that have
slow amplitudes and phases (Cadieu & Olshausen, 2012).

In video processing and computer vision, a common strategy for temporal
prediction is to first estimate local translational motion, and to then
(assuming no acceleration) use this to warp and/or copy previous content to
predict the next frame. Such motion compensation is a fundamental component in
video compression schemes like MPEG (Wiegand et al., 2003). These video coding
standards are the result of decades of engineering efforts, and have enabled
reliable and efficient digital video communication that is now commonplace.
But motion estimation is a difficult nonlinear problem, and existing methods
fail in regions where temporal evolution is not translational and smooth: for
example, expanding or rotating motions, discontinuous motion at occlusion
boundaries, or mixtures of motion arising from semi-transparent surfaces
(e.g., viewing the world through a dirty pane of glass). In compression
schemes, these failures of motion estimation lead to prediction errors, which
must then be repaired by sending additional corrective bits.

Human perception does not seem to suffer from such failures - subjectively, we
can anticipate the time-evolution of visual input even in the vicinity of
these commonly occurring non-translational changes. In fact, those changes are
often highly informative, as they reveal object boundaries, and provide
ordinal depth and other information about the visual scene. This suggests that
the human visual system uses a different strategy, perhaps bypassing
altogether the estimation of local motion, to represent and predict evolving
visual input. Toward this end, and inspired by recent hypotheses that primate
visual representations support prediction by “straightening” the temporal
trajectories of naturally-occurring input (Hénaff et al., 2019), we formulate
an objective for learning an image representation that facilitates prediction
by linearizing the temporal trajectories of frames of natural video.

To motivate the separation of instantaneous spatial representation from
temporal prediction, we first consider the special case of rigidly translating
video content. When expressed in the frequency domain, translation corresponds
to linear phase advancement (section 2.1), and prediction of rigidly
translating content reduces to angular extrapolation (section 2.2). We
generalize this factorization using group representation theory (section 2.3),
and describe a neural network architecture that maps individual video frames
to a latent complex-valued representation. Within this latent space,
coefficients can be temporally predicted by phase advancement and then mapped
back to generate an estimated frame. The entire systems may then be trained
end-to-end to minimize next frame prediction errors (section 3). We report
training results of several such systems, and show that they produce
systematic improvements in predictive performance over both conventional
motion compensation methods, and direct predictive neural networks (section
4). Finally, we relate this approach to previous work (section 5) and discuss
its significance and implications (section 6).

##  2 Background

###  2.1 Base case: the Fourier shift theorem

Our approach is motivated by the well-known behavior of Fourier
representations with respect to signal translation (note that this elementary
example will later lead to our proposed generalization). Specifically, the
complex exponentials that make up the Fourier basis are the eigenfunctions of
the translation operator, and translation of inputs produces systematic phase
advances of frequency coefficients. Let $x\in\mathbb{R}^{N}$ be a discrete
signal indexed by spatial location $n\in[0,N-1]$, and let
$\widetilde{x}\in\mathbb{C}^{N}$ be its Fourier transform indexed by
$k\in[0,N-1]$. We write $x^{v}(n)=x(n-v)$, the translation of $x$ by $v$
modulo $N$ (ie. circular shift with period N). Defining $\phi=e^{i2\pi/N}$,
the primitive N-th root of unity, and $\mathcal{F}_{nk}=\phi^{nk}$, the
$N\times N$ Fourier matrix, we can express the Fourier shift theorem111Proof
by substituting $m=n-v$:  $\displaystyle\widetilde{x^{v}}(k)$
$\displaystyle=\sum_{n=0}^{N-1}\phi^{-kn}x(n-v)=\sum_{m=-v}^{N-1-v}\phi^{-kv}\phi^{-km}x(m)$
$\displaystyle=\phi^{-kv}\sum_{n=0}^{N-1}\phi^{-kn}x(n)=\phi^{-kv}\widetilde{x}(k).$
as: $x^{v}(n)=\frac{1}{N}\mathcal{F}D(v)\mathcal{F}^{*}x(n),$ where
$\mathcal{F}^{*}$ is the conjugate transpose of the Fourier matrix and
$D(v)=\textrm{diag}(\phi^{0},...,\phi^{-(n-1)v})$ is a diagonal matrix.

This relationship may be depicted in a compact diagram:

${\widetilde{x}(k)}$${\phi^{-kv}\widetilde{x}(k)}$${x(n)}$${x(n-v).}$$\scriptstyle{\rm
advance\;\;phase}$$\scriptstyle{\mathcal{F}}$$\scriptstyle{\mathcal{F^{*}}}$$\scriptstyle{\rm
shift}$ (1)

In the context of our goals, the diagram illustrates the point that
transforming to the frequency domain renders translation a “simpler”
operation: a phase advance is a rotation in the two dimensional (complex)
plane.

Figure 1:  Translation of a 1D signal consisting of a sum of two sinusoidal
components: $x(n,t)=\sin(2\pi(n-t))+\sin(2\pi 3(n-t))/2$. Lower left: three
snapshots of the signal as it translates. Lower right: In the high-dimensional
space representing the signal (each axis corresponding to the signal value at
one location), the temporal trajectory is highly curved. Shown is the
projection of the signal vector into the 3D space of the top three principal
components. Three colored points indicate the three snapshots in lower left
panel. Upper left: Fourier transform of the signal, showing complex-valued
coefficients as a function of frequency. In this representation the temporal
trajectory corresponds to linearly increasing phase of the two sinusoidal
components, each at a rate proportional to its frequency. Upper right: a polar
coordinate transform to amplitude and phase of each frequency component leads
to a representation that evolves along a straight line, and is thus readily
predictable (phases are unwrapped for display purposes).

###  2.2 Prediction via angular extrapolation

Now consider observations of a signal that translates at a constant velocity
over time, $x(n,t)=y(n-vt)$. Although the temporal evolution is easy to
describe, it traces a highly non-linear trajectory in the signal state space,
rendering prediction difficult (specifically, linear extrapolation fails). As
an example, Figure 1 shows a signal consisting of a sum of two sinusoidal
components. Transforming the signal to the frequency domain simplifies the
description. In particular, the translational motion now corresponds to
circular motion of the two (complex-valued) Fourier coefficients associated
with the constituent sinusoids.

The motion is further simplified by a polar coordinate transform to extract
phase and amplitude of each Fourier coefficient. Specifically, the motion is
now along a straight trajectory, with both phases advancing linearly (but at
different rates), and both amplitudes constant. Note that this is a geometric
property that holds for any rigidly translating signal, and offers a simple
means of predicting content over time. Indeed, we can use the shift property
(see section 2.1) on $x(n,t+1)=x^{v}(n,t)$ and observe that prediction is now
reduced to linear extrapolation of each coefficient’s phase. We have the three
step process:

$\displaystyle\widetilde{x}(k,t)$
$\displaystyle=\sum_{n=0}^{N-1}\phi^{-kn}x(n,t),$ (analyze)
$\displaystyle\widetilde{x}(k,t+1)$
$\displaystyle=\phi^{-kv}\widetilde{x}(k,t),$ (advance phase) $\displaystyle
x(n,t+1)$
$\displaystyle=\frac{1}{N}\sum_{k=0}^{N-1}\phi^{nk}\widetilde{x}(k,t+1).$
(synthesize)

Since we assumed that the motion from time $t$ to $t+1$ is identical to that
from time $t-1$ to $t$ (ie. no acceleration), the phase advance $kv$ can be
computed from the past two representations as
$kv=\angle\widetilde{x}(k,t)-\angle\widetilde{x}(k,t-1)$, where $\angle z$
indicates the phase of the complex number $z$. Thus, a polar coordinate
transformation in the frequency domain converts translational motion into
trajectories that are predictable via linear phase extrapolation.

###  2.3 Generalization: representing commutative Lie groups

Natural videos are replete with rich temporal transformations, such as
continuous deformations of objects and surfaces. Assuming that these
transformations can be described as groups, we will aim to learn their group
representation from data. To this end, we seek a parameterization that
generalizes beyond translation and the frequency domain. Remarkably, Fourier
analysis can be seen as a special case of the representation theory of compact
commutative Lie (ie. smooth) groups (Mackey, 1980).

In harmonic analysis, the celebrated Peter-Weyl Theorem (1927) establishes the
completeness of the irreducible representations of any compact continuous
group (an irreducible representation is a subspace that is invariant to group
action and that can not be further decomposed). Furthermore, it follows that
every compact Lie group admits a faithful (ie. injective) representation given
by an explicit complete orthogonal basis, constructed from finite-dimensional
irreducible representations (Hall, 2013). Accordingly, the action of a compact
Lie group can be expressed as a rotation within each irreducible
representation - thereby generalizing the Fourier shift property (an example
is the construction of steerable filters (Freeman et al., 1991) in the
computational vision literature).

In the case of compact commutative Lie groups, the irreducible representations
are one-dimensional and complex valued: they are pairs of real valued basis
functions. Therefore, the angular extrapolation mechanism described in the
previous section (2.2) can be employed for prediction in a much wider setting
than that of translational motion. We will rely on the parameterization
suggested by the representation theory of compact commutative Lie groups to
learn the harmonic basis functions of the transformations at play in natural
videos.

##  3 Learning to predict with angular extrapolation

To generalize beyond translation and the Fourier transform, we aim to learn a
representation of video frames that enables prediction via angular
extrapolation. Specifically, we focus on next frame prediction, and optimize
two parameterized mappings: one for the analysis and one for the synthesis
transform. This framework is illustrated in Figure 2, which provides a
generalization of the Fourier shift diagram (1).

Figure 2:  Unsupervised predictive representation learning framework. Each
video frame is transformed using a parametric mapping $f_{w}$, to an internal
representation consisting of pairs of coefficients arranged in spatial
channels. Predictions of individual complex coefficients at time $t+1$ are
computed by advancing the phase of the current coefficients by an amount equal
to the phase advance over the interval from $t-1$ to $t$. At each time step,
one such coefficient is depicted as a vector in two dimensions and the top
arrows indicates how they are combined (each vector corresponds to the complex
coefficient at a particular location within one channel pair). Predicted
frames are then generated by applying the parameterized inverse mapping
$g_{w}$ on the advanced coefficients. Forward and inverse mappings are jointly
trained to minimize mean squared prediction error between the predicted and
actual frame at time $t+1$.

###  3.1 Architecture and objective function

When focusing on a small region in an image sequence, the transformation
observed as time passes can be approximated as a local translation. That is to
say, in a spatial neighborhood around position $n$, $m\in N(n)$, we have:
$x(m,t+1)\approx x(m-v,t)$. We can use the decomposition described for global
rigid translation, replacing the Fourier transform with a local convolutional
operator (Fleet & Jepson, 1990), processing each spatial neighborhood of the
image independently and in parallel, and applying angular extrapolation to the
coefficients computed at each position.

We use the same weights for the encoding and decoding stages, that is to say
the analysis operator is the transpose of the synthesis operator (also true of
the Fourier transform and its inverse). Sharing these weights reduces the
number of parameters and simplifies interpretation of the learned solution.
This “polar predictor” (hereafter, PP) is consistent with the general scheme
described in figure 2 where $f_{w}$ is taken to be linear and convolutional,
and $g_{w}$ is its transpose. In practice, we assumed 64 convolutional
channels with filters of size $17\times 17$ pixels, with no additive
constants.

At every position in the image (spatial indices are omitted for clarity of
notation), each coefficient $y_{j}(t)$ is computed as an inner product between
the input $x(t)$ and the filter weights $w_{j}$ of each channel
$j\in\left[0,63\right]$ : $y_{j}(t)=w_{j}^{T}x(t)$. In order to obtain phases,
we combine coefficients in pairs, indexed by $k\in\left[0,31\right]$, which
can be written as single complex coefficient as:
$z_{k}(t)=y_{2k}(t)+iy_{2k+1}(t)\in\mathbb{C}$, and expressed in polar
coordinates as: $z_{k}(t)=a_{k}(t)e^{i\theta_{k}(t)}$. This polar coordinate
transformation is the only non-linear step used in the PP architecture, and
serves as a bivariate non-linear activation function, differing markedly from
the typical (pointwise) rectification operations found in convolutional neural
networks.

With this notation, linear phase extrapolation reduces to
$\hat{z}_{k}(t+1)=a_{k}(t)e^{i(\theta_{k}(t)+\Delta\theta_{k}(t))}$, where the
phase advance $\Delta\theta_{k}(t)$ is equal to the phase difference over the
interval from $t-1$ to $t$:
$\Delta\theta_{k}(t)=\theta_{k}(t)-\theta_{k}(t-1)$. The advanced coefficients
can be written in a more compact form, using complex arithmetic, as:

$\hat{z}_{k}(t+1)=\frac{z_{k}(t)^{2}\overline{z_{k}(t-1)}}{|z_{k}(t)||z_{k}(t-1)|},$
(2)

where $\overline{z}$ and $|z|$ respectively denote complex conjugation and
complex modulus of $z$. This formulation in terms of complex coefficients has
the benefit of handling phases implicitly, bypassing the discontinuities of
phase unwrapping and the instability of angular variables (phase is unstable
when amplitude is low). We find that such an indirect formulation of phase
processing is necessary for the stability of training, as previously noted in
the texture modeling literature (Portilla & Simoncelli, 2000). Finally, the
estimated next frame is generated by applying the transposed convolution
$g_{w}$ (with the same weights as $f_{w}$) to the advanced coefficients.

As a more substantial generalization of polar prediction, we use deep
convolutional neural networks to instantiate non-linear mappings for both the
encoder $f_{w}$ and the decoder $g_{w}$ (each with independent filters).
Specifically, the “deep polar predictor” (deepPP) operates by transforming two
frames of input into the encoding space, $z(t-1)=f_{w}(x(t-1))$ and
$z(t)=f_{w}(x(t))$, applying the polar prediction of equation 2 to this
encoded representation, and then decoding the next frame from this prediction,
$\hat{x}(t+1)=g_{w}(\hat{z}(t+1))$. While the PP model learns a linear
representation, the deepPP model is nonlinear, with potential to enhance
prediction by adapting to signal properties.

In order to isolate the effects of non-linearities from those of spatial
scale, we chose the number of layers and the kernel sizes of deepPP so that
the effective receptive field size was matched to that of the PP model.
Specifically, both the encoder and the decoder are composed of 4 convolutional
layers, each with 64 channels, and using filter kernels of size $5\times 5$
followed by half-wave rectification (ReLU).

For both PP and deepPP models, convolutional kernels $w$ are learned by
minimizing the average squared prediction error:

$\min_{w}\mathbb{E}_{t}{\left\lVert x(t+1)-\hat{x}(t+1)\right\rVert}_{2}^{2}.$

The computation of this prediction error is restricted to the center of the
image because moving content that enters from outside the video frame is
inherently unpredictable. Specifically, we trim a 17-pixel strip from each
side. Note that we only perform valid convolutions to avoid artificial
interference with prediction (zero-padding creates undesirable boundary
artifacts).

###  3.2 Comparison models

We compare our method to the traditional motion-compensated coding approach
that forms the core of inter-picture coding in well established compression
standards such as MPEG. Block matching is an essential component of these
standards, allowing the compression of video content by up to three orders of
magnitude with moderate loss of information. For each block in a frame,
typical coders search for the most similar spatially displaced block in the
previous frame (typically measured with MSE), and communicate the displacement
coordinates to allow prediction of frame content by translating blocks of the
(already transmitted) previous frame. We implemented a “diamond search”
algorithm (Zhu & Ma, 2000) operating on blocks of $8\times 8$ pixels, with a
maximal search distance of 8 pixels which balances accuracy of motion
estimates and speed of estimation (the search step is computationally
intensive). We use the estimated displacements to perform causal motion
compensation (cMC), using displacement vectors estimated from the previous two
observed frames ($x_{t-1}$ and $x_{t}$) to predict the next frame ($x_{t+1}$)
rather than the current one (as in MPEG).

To isolate the effects of the polar prediction, we also implemented a
predictor using linear extrapolation of the responses of a deep neural network
(deepL), with architecture identical to that of the deep polar predictor. That
is to say, we replace equation 2 by: $\hat{y}_{j}(t+1)=2y_{j}(t)-y_{j}(t-1)$,
which amounts to enforcing linear dynamics in the latent space of the non-
linear representation.

Finally, we implemented a more direct convolutional neural network predictor
(CNN), that maps two successive observed frames to an estimate of the next
frame (Mathieu et al., 2016). This predictor jointly transforms and predicts
visual signals without explicitly partitioning spatial content representation
and temporal feature extrapolation. For this, we used a CNN composed of 20
stages, each consisting of 64 channels, and computed with $3\times 3$ filters
without additive constants, followed by half-wave rectification. Note that,
unlike all other predictors, this model jointly processes pairs of frames to
generate predictions.

###  3.3 Datasets and training

To train, test and compare these models, we use the DAVIS dataset (Pont-Tuset
et al., 2017), which was originally designed as a benchmark for video object
segmentation. Image sequences in this dataset contain diverse motion of scenes
and objects (eg., with fixed or moving camera, and objects moving at different
speeds and directions), which make next frame prediction challenging. Each
clip is sampled at 25 frames per second, and is approximately 3 seconds long.
The set is subdivided into 60 training videos (4741 frames) and 30 test videos
(2591 frames).

We pre-processed the data, converting all frames to monochrome luminance
values, and scaling their range to the interval $[-1,1]$. Frames are cropped
to a $256\times 256$ central region, where most of the motion tends to occur,
and then spatially down-sampled to $128\times 128$ pixels. We assume the
temporal evolution of natural signals to be sufficiently and appropriately
diverse for training, and do not apply any additional data augmentation
procedures. We train on brief temporal segments containing 11 frames, which
allows for prediction of 9 frames, processing these in batches of size 4. We
train each model for one hundred epochs using the Adam optimizer (Kingma & Ba,
2015) with default parameters and a learning rate of $3\cdot 10^{-4}$. The
learning rate is halved at epochs 50, 60, 70, 80, 90, 100. We use batch
normalization before every half-wave rectification, rescaling by the standard
deviation of channel coefficients (but with no additive terms).

Similarly, we also trained on the larger UCF-101 dataset (Soomro et al.,
2012). This dataset, initially designed for action recognition, contains about
2.5 million frames, which amounts to over 27 hours of video data. Note that,
unlike the DAVIS dataset, the clips are only available in compressed video
formats and may contain motion artifacts (due to inter-frame coding). We used
the same pre-processing procedure, except that we reduced frames by directly
cropping a $128\times 128$ central region (without any down-sampling). We
employ the same training procedure, except that we only run training for 25
epochs.

##  4 Unsupervised representation learning

###  4.1 Recovery analysis

To experimentally validate our approach, we first verified that a the PP model
can robustly recover known symmetries in small synthetic datasets consisting
of translating or rotating image patches. For these experiments, we applied
encoding and decoding transforms to the entire patch (i.e., non-
convolutionally). Wen trained on translating image patches, the PP model
learned approximately sinusoidal filters, shifted in phase by $\pi/2$ \- i.e.,
a local Fourier transform. Similarly, when trained on rotating patches, the
learned filters represented circular harmonics. We also found that PP extracts
meaningful representations when multiple kinds of transformations are at play
(eg. mixing both translations and rotations), and when the transformation are
not perfectly translational (eg. translation with open boundary condition).
Learned filters for each of these cases are provided in Figure 5 in the
appendix.

###  4.2 Performance on Natural Videos

Table 1: Prediction error computed on the DAVIS and UCF-101 datasets. Values indicate average Mean Squared Error.  Algo. | DAVIS | UCF-101 | # param.
---|---|---|---
| train | test | train | test |
Copy | $0.064$ | $0.065$ | $0.0302$ | $0.0286$ |
cMC | $0.048$ | $0.049$ | $--$ | $0.0299$ |
deepL | $0.034$ | $0.037$ | $0.0220$ | $0.0217$ | $665,856$
CNN | $0.031$ | $0.035$ | 0.0210 | $0.0215$ | $666,496$
PP | $0.036$ | $0.035$ | $0.0245$ | $0.0229$ | $18,496$
deepPP | 0.028 | 0.032 | $0.0216$ | 0.0210 | $665,856$

(a)

(b)

(c)

(d)

Figure 3: Detailed performance comparison (in Root Mean Squared Error) of
predictive algorithms. Each point corresponds to a frame in the test set.
Vertical axes represent difference in performance: For points lying below the
horizontal axis, the method whose performance is plotted on the horizontal
axis achieves a lower RMSE than the comparison method. Black crosses indicate
average RMSE. Red point corresponds to example in Figure 4, green to Figure 8
and blue to Figure 9.

We summarize the main prediction results in Table 1. First, observe that the
predictive algorithms considered in this study perform significantly better
than the baseline obtained by simply copying the last frame. Second, the polar
predictor (PP) performs nearly as well as the convolutional neural network
(CNN) in terms of test mean squared error on DAVIS. This demonstrates the
remarkable power of the polar predictor: the PP model has roughly 30 times
fewer parameters and uses a single non-linearity, while the CNN is composed of
20 non-linear layers. Finally, observe that deepPP achieves the lowest mean
squared error, notably outperforming the deepL model which uses linear
extrapolation in an otherwise identical architecture. It also outperforms the
CNN on the UCF-101 test dataset while remaining significantly simpler. Thus,
the prediction task benefits substantially from use of a fixed nonlinear
transformation to polar coordinates.

While average performance values provide a compact summary, it is also
informative to examine the distribution of prediction errors on individual
frames from the test set. Figure 3 shows pairwise comparison of the predictive
algorithms for each frame in the DAVIS dataset. To make the contrast more
apparent, we display performance difference on the vertical axis. Note that
while the models have been optimized to reduce mean squared error, we show
root mean squared error (RMSE) in order to facilitate visual inspection of the
results (the concavity of the square root spreads out small differences). We
see that (a) The polar predictor representation systematically outperforms
causal motion compensation, especially on difficult examples. (b) The polar
predictor outperforms the CNN on the bulk of easy to medium cases but this
advantage is reversed for harder examples. (c) The deep polar predictor
outperforms the single layer polar predictor overall, indicating that the non-
linearity in representation can help. (d) The deep polar predictor clearly
outperforms the deep linear predictor, revealing the strong benefit of using a
polar extrapolation mechanism over a linear one.

###  4.3 Learned filters

In order to better understand these results, we visualized the learned PP
filters trained on the DAVIS dataset and observed that the learned filter are
selective for orientation and spatial frequency, and that that they tile the
frequency domain. Filters in each pair have a similar frequency preference,
and are related by a 90 degrees phase shift (see in Figure 6(a) and 6(b) in
the Appendix). This relationship is analogous to that of sines and cosines and
is consistent with the structure of the angular extrapolation described in
equation 2.

###  4.4 Examples

Figure 4: A typical example image sequence from the DAVIS test set. The first
three frames on the top row display the unprocessed images, and last five
frames show the respective prediction for each method (with their shorthand
above). The bottom row displays error maps computed as the difference between
the target image $x(t+1)$ and each predicted next frame on the corresponding
position in the first row. Images, predictions and error maps are all shown on
the same scale.

Consider a set of example videos, chosen to illustrate behaviors of the
methods being compared. In Figure 4, we see a wall, its shadow and their sharp
boundaries against a grass background as the camera moves. Both PP and deepPP
generate good results, cMC produces a sharp prediction at the expense of
significant blocking artifacts, and both the CNN and deepL tend toward
excessive blurring. Here, the cMC is is significantly sharper than the others,
but introduces substantial artifacts. Again, the PP methods produce sharper
results than either the CNN or deepL methods. A few additional informative
examples are displayed in the appendix (see Figure 8, 9, 10).

##  5 Related work

Our method is conceptually related to sparse coding with complex-valued
coefficients (Cadieu & Olshausen, 2012) in that it factorizes natural videos
into form and motion. But it differs in a number of important ways: (1) sparse
coding focuses on representing, not predicting, the signal; (2) we do not
promote sparsity of either amplitude or phase components; (3) finally, the
discontinuity arising from selection of sparse subsets of coefficients seems
at odds with the representation of continuous group actions, while our
explicit mapping into polar coefficients aims for a smooth and continuous
parameterization of the transformations that occur in natural videos. Several
other studies have aimed to learn representations that decompose signal
identity and attribute (ie. a what-where, or invariance-equivariance
factorization). In particular learning linearized features from video was
explored using a heuristic extrapolation mechanism (Goroshin et al., 2015).
The authors developed specialized “soft max-pooling” and “soft argmax-pooling”
modules and tested them on the small NORB dataset. A related approach aimed at
finding video representations which decompose content and pose in order to
enable prediction (Hsieh et al., 2018). This work explicitly identifies
spatial components that are easier to predict in the moving MNIST and bouncing
balls datasets. More sophisticated architectures have been developed to tackle
the challenge of natural video prediction. In particular, a recurrent
instantiation of the predictive coding theory (Rao & Ballard, 1999) introduced
a stacked convolutional LSTM architecture (Lotter et al., 2017). In contrast,
our framework scales to prediction of natural videos while remaining simple:
we rely on principles of signal processing and representation theory to employ
a polar non-linearity (and we describe an effective and stable
implementation), but we do not explicitly model the stochastic nature of the
video prediction problem.

Our method is also related to work that adopts a Lie group formalism in
representation learning. Since the seminal work that proposed learning Lie
group generators from dynamic signals (Rao & Ruderman, 1998), the polar
parametrization was explored in (Cohen & Welling, 2014) to identify
irreducible representations in a synthetic dataset. The continuous group
formalism has also been combined with sparse coding (Chen et al., 2018; Chau
et al., 2020) to model natural images as points on a latent manifold. More
recently, bispectral neural networks (Sanborn et al., 2022) have been shown to
learn image representations invariant to a given global transformation (in
particular cyclic translation and rotation of MNIST digits). In contrast to
the coding approach, our formulation relies on a prediction objective to
jointly discover and exploit the symmetries implicit in data. In order to
scale to natural video data, where multiple unknown and noisy transformations
are play, we developed a convolutional approach that adapts to the local
structure of transformations. This formulation can represent a very large
family of local symmetries (including diffeomorphisms and non-smooth fields of
local translations). This generality comes at the cost of precisely
identifying what groups of transformations are captured by the learned
representation.

Finally, in the fluid mechanics literature, the Koopman operator approach
(Mezić, 2005) has been used to lift a system from its original state-space to
a higher dimensional representation space where its dynamics can be linearized
- a dynamical analog of the well known kernel trick. This formalism has
spurred a line of work in machine learning that relies on autoencoders to
learn coordinate systems that approximately linearize dynamics (Lusch et al.,
2018; Azencot et al., 2020). In this perspective, our work can also be
interpreted as learning the spectral properties of an abstract Koopman
operator operating on video data, specifically estimating its complex
eigenvectors. Our approach makes an inertial assumption and does not require
an auxiliary network to compute velocities. Moreover it relies on a
convolutional approach and is able to predict raw videos (which tend to
contain richer structure than typical fluid flows).

##  6 Discussion

We have presented a simple self-supervised representation-learning framework
based on next-frame prediction. It unveils the temporal structure of natural
videos using local polar coordinates. Our approach jointly discovers and
exploits the local symmetries present in the temporal evolution of image
sequences, in particular the spatio-temporal redundancies due to local
deformation of image content. We assumed that spatial processing and temporal
extrapolation can be partitioned into i) the learned parameterized mappings,
one that extract pairs of local features from individual frames and one that
generates a frame from the coefficients; and ii) a fixed angular extrapolation
mechanism that advances coefficients and embodies an inertial hypothesis (ie.
evolving content will continue to evolve in the same way). Our empirical
results demonstrate that these assumptions, far from being too limiting,
correspond well to the structure of natural videos and provide a natural
representation thereof.

Specifically, we used the the polar coordinate transformation as a bivariate
non-linear activation function acting on pairs of coefficients in the
representation. Predictions in this representation were computed by phase
advancement, which was implemented implicitly (Eq. 2). Compared to linear
extrapolation, angular extrapolation achieved higher prediction accuracy on
natural video. Using terminology from group theory, our polar models factorize
signals into an invariant part, which is stable in time, and an equivariant
part, which evolves linearly. This choice of prediction mechanism, motivated
by principles of signal processing and harmonic analysis, acts as a structural
prior. Although the conventional deep convolutional network (CNN) considered
here could in principle have discovered this solution, it failed to do so
(within the constraints of our architecture). The polar predictor, on the
other hand, is well-matched to the task, and achieves a good solution using
only a fraction of the number of parameters. It is optimized on a mean squared
error objective, without any other additional regularization - which
facilitates interpretability. This exemplifies a fundamental theme in
computational vision and machine learning: when possible, let the
representation do the analysis.

Our approach to prediction has the advantage of being motion-informed while
not relying on explicit motion estimation. Because it is not constrained to
assigning a single motion vector at every location and instead represents a
distribution of phases, this method bypasses known difficulties of motion
estimation in handling non-translational motions and outperforms a
conventional causal motion compensated algorithm. In the era of GPU computing,
it admits a very fast implementation that has potential for applications in
video compression. Moreover, the polar predictor takes the form of a
predictive auto-encoder that associates a latent representation vector to each
frame. This representation may prove useful for other tasks like object
categorization, segmentation, or estimation of heading direction for a moving
observer. Several natural extensions of the work presented here can be further
explored: (i) treating the angular extrapolation prediction mechanism as a
more general building block that is cascaded in a multi-layer architecture;
(ii) optimizing representation layers deeper in the hierarchy to make
predictions at longer timescales; (iii) measuring prediction error directly in
the representation domain, while avoiding representation collapse - such a
local objective function would allow a potential connection with biological
neural architecture and to human visual perception. (iv) examining and
interpreting what is learned in the deepPP model, especially around occlusion
boundaries (which is not invertible, and therefore not a group action).

## References

  * Azencot et al. (2020) Azencot, O., Erichson, N. B., Lin, V., and Mahoney, M.  Forecasting sequential data using consistent koopman autoencoders.  In _International Conference on Machine Learning_ , 2020.
  * Cadieu & Olshausen (2012) Cadieu, C. F. and Olshausen, B. A.  Learning intermediate-level representations of form and motion from natural movies.  _Neural computation_ , 2012.
  * Chau et al. (2020) Chau, H. Y., Qiu, F., Chen, Y., and Olshausen, B.  Disentangling images with lie group transformations and sparse coding.  _arXiv preprint arXiv:2012.12071_ , 2020.
  * Chen et al. (2018) Chen, Y., Paiton, D., and Olshausen, B.  The sparse manifold transform.  _Advances in neural information processing systems_ , 31, 2018.
  * Cohen & Welling (2014) Cohen, T. and Welling, M.  Learning the irreducible representations of commutative lie groups.  In _International Conference on Machine Learning_. PMLR, 2014.
  * Fleet & Jepson (1990) Fleet, D. J. and Jepson, A. D.  Computation of component image velocity from local phase information.  _International journal of computer vision_ , 1990.
  * Földiák (1991) Földiák, P.  Learning invariance from transformation sequences.  _Neural Computation_ , 1991.
  * Freeman et al. (1991) Freeman, W. T., Adelson, E. H., et al.  The design and use of steerable filters.  _IEEE Transactions on Pattern analysis and machine intelligence_ , 1991.
  * Goroshin et al. (2015) Goroshin, R., Mathieu, M. F., and LeCun, Y.  Learning to linearize under uncertainty.  In _Advances in Neural Information Processing Systems_ , 2015.
  * Hall (2013) Hall, B. C.  Lie groups, lie algebras, and representations.  In _Quantum Theory for Mathematicians_. Springer, 2013.
  * Hénaff et al. (2019) Hénaff, O. J., Goris, R. L., and Simoncelli, E. P.  Perceptual straightening of natural videos.  _Nature neuroscience_ , 2019.
  * Hsieh et al. (2018) Hsieh, J.-T., Liu, B., Huang, D.-A., Fei-Fei, L. F., and Niebles, J. C.  Learning to decompose and disentangle representations for video prediction.  _Advances in neural information processing systems_ , 2018.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _ICLR_ , 2015.
  * Lotter et al. (2017) Lotter, W., Kreiman, G., and Cox, D.  Deep predictive coding networks for video prediction and unsupervised learning.  _International Conference on Learning Representations_ , 2017.
  * Lusch et al. (2018) Lusch, B., Kutz, J. N., and Brunton, S. L.  Deep learning for universal linear embeddings of nonlinear dynamics.  _Nature communications_ , 2018.
  * Mackey (1980) Mackey, G. W.  Harmonic analysis as the exploitation of symmetry–a historical survey.  _Bulletin (New Series) of the American Mathematical Society_ , 3(1. P1):543–698, 1980.
  * Mathieu et al. (2016) Mathieu, M., Couprie, C., and LeCun, Y.  Deep multi-scale video prediction beyond mean square error.  _ICLR_ , 2016.
  * Mezić (2005) Mezić, I.  Spectral properties of dynamical systems, model reduction and decompositions.  _Nonlinear Dynamics_ , 2005.
  * Pont-Tuset et al. (2017) Pont-Tuset, J., Perazzi, F., Caelles, S., Arbeláez, P., Sorkine-Hornung, A., and Van Gool, L.  The 2017 davis challenge on video object segmentation.  _arXiv:1704.00675_ , 2017.
  * Portilla & Simoncelli (2000) Portilla, J. and Simoncelli, E. P.  A parametric texture model based on joint statistics of complex wavelet coefficients.  _International journal of computer vision_ , 2000.
  * Rao & Ruderman (1998) Rao, R. and Ruderman, D.  Learning lie groups for invariant visual perception.  _Advances in neural information processing systems_ , 1998.
  * Rao & Ballard (1999) Rao, R. P. and Ballard, D. H.  Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.  _Nature neuroscience_ , 1999.
  * Sanborn et al. (2022) Sanborn, S., Shewmake, C., Olshausen, B., and Hillar, C.  Bispectral neural networks.  _arXiv preprint arXiv:2209.03416_ , 2022.
  * Soomro et al. (2012) Soomro, K., Zamir, A. R., and Shah, M.  Ucf101: A dataset of 101 human actions classes from videos in the wild.  _arXiv preprint arXiv:1212.0402_ , 2012.
  * Wang et al. (2004) Wang, Z., Bovik, A. C., Sheikh, H. R., and Simoncelli, E. P.  Image quality assessment: from error visibility to structural similarity.  _IEEE transactions on image processing_ , 2004.
  * Wiegand et al. (2003) Wiegand, T., Sullivan, G. J., Bjontegaard, G., and Luthra, A.  Overview of the h. 264/avc video coding standard.  _IEEE Transactions on circuits and systems for video technology_ , 2003.
  * Wiskott & Sejnowski (2002) Wiskott, L. and Sejnowski, T. J.  Slow feature analysis: Unsupervised learning of invariances.  _Neural computation_ , 2002.
  * Zhu & Ma (2000) Zhu, S. and Ma, K.-K.  A new diamond search algorithm for fast block-matching motion estimation.  _IEEE transactions on Image Processing_ , 2000.

##  Appendix A Appendix

(a) Translation (cyclic boundary)

(b) Rotation

(c) Translation and rotation

(d) Translation (open boundary)

Figure 5: Filters of polar predictor networks trained to predict small
synthetic sequences. We randomly select 100 image patches of size $16\times
16$ from the DAVIS dataset and generate training data by manually transforming
them - applying translations or rotations. We verify that PP recovers the
known harmonic functions: Fourier modes for translation (panel 5(a)), and disk
harmonics for rotation (panel 5(b)). To show that the recovery of harmonics is
robust, we design two additional synthetic datasets. i) the combination of
translational and rotational sequences. In this case, PP learns filters that
correspond to either group, suggesting that our approach can generalize to
situations with more than than one group at play (panel 5(c)); ii) generalized
translation sequences: spatially sliding a square window on a large image (ie.
new content creeps in and falls off at boundaries), instead of using cyclic
boundary condition (ie. content wraps around the edges). In this case, PP
learns localized Fourier-like modes (panel 5(d)), indicating that approximate
group actions still provide meaningful training signal - although it is much
more noisy. In each panel, the 32 pairs of filters are sorted by their norm.
Notice that some of the filters are not structured and generally miss high
frequency harmonics. This is due to the spectral properties of the datasets,
which have more power at lower frequencies, and to the discretization of the
transformations.

(a) spatial domain filters

(b) corresponding Fourier amplitude spectra

Figure 6: Filters of a polar predictor trained to predict natural videos from
the DAVIS dataset. The 32 pairs of convolutional filters are sorted by their
norm and their amplitude spectrum is displayed at corresponding locations on
the right panel. Observe that the filters are selective for orientation and
spatial frequency, tile the frequency spectrum, and form quadrature pairs.
Notice that some of the learned filters do not exactly conform to the
idealized description just given. Table 2: Prediction error computed on the
DAVIS dataset. Values indicate average PSNR and SSIM. All methods obtain
relatively low structural similarity scores (Wang et al., 2004), a perceptual
measure of similarity that equals one when both images are identical. This
indicates that prediction on this dataset is quite challenging. Algorithm
---
metric | set | Copy | cMC | deepL | CNN | PP | deepPP
PSNR $\uparrow$ | train | 21.32 | 23.82 | 23.18 | 23.84 | 24.49 | 24.52
| test | 20.06 | 22.37 | 22.30 | 22.82 | 23.46 | 23.35
SSIM $\uparrow$ | train | 0.52 | 0.64 | 0.58 | 0.62 | 0.65 | 0.64
| test | 0.50 | 0.62 | 0.55 | 0.59 | 0.63 | 0.61
Figure 7: Another comparison - see caption of Fig. 3. The deep polar predictor
generally outperforms the CNN over the whole range of difficulties. Figure 8:
Another example image sequence with nonrigid motions/deformations - see
caption of Fig. 4. As the biker advances to the right, the camera tracks and
leads its displacement. Figure 9: Another example image sequence - see caption
of Fig. 4. A video with one portion not moving. In these regions the cMC is
most effective, and PP and deepPP outperform CNN and deepL. Figure 10: Another
image sequence - see caption of Fig. 4. As expected, in presence of large non-
rigid motion, with texture and object deformations, prediction is very
difficult and none of the tested methods performs well.


* [40] E. de Klerk and M. Laurent.  A survey of semidefinite programming approaches to the generalized problem of moments and their error analysis. In World Women in Mathematics 2018, Araujo C., Benkart G., Praeger C., Tanbay B. (eds). Association for Women in Mathematics Series, Springer, Cham, vol. 20, pages 17–56, 2019.
  * [41] M. Korda.  Stability and performance verification of dynamical systems controlled by neural networks: Algorithms and complexity.  IEEE Control Systems Letters, 6:3265–3270, 2022.
  * [42] M. Korda and C. N. Jones.  Stability and performance verification of optimization-based controllers.  Automatica, 78:34–45, 2017.
  * [43] J. B. Lasserre.  Global optimization with polynomials and the problem of moments.  SIAM J. Optim., 11:796–817, 2001.
  * [44] J. B. Lasserre.  Convergent SDP-relaxations in polynomial optimization with sparsity.  SIAM Journal on Optimization, 17(3):822–843, 2006.
  * [45] J. B. Lasserre.  A semidefinite programming approach to the generalized problem of moments.  Mathematical Programming, 112:65–92, 2008.
  * [46] J. B. Lasserre.  Moments, Positive Polynomials and Their Applications.  Imperial College Press, 2009.
  * [47] J. B. Lasserre.  An Introduction to Polynomial and Semi-Algebraic Optimization,  Cambridge University Press, 2015.
  * [48] J. B. Lasserre and Y. Emin.  Semidefinite relaxations for Lebesgue and Gaussian measures of unions of basic semialgebraic sets.  Mathematics of Operations Research, 44(4):1477–1493, 2019.
  * [49] M. Laurent.  Sums of squares, moment matrices and optimization over polynomials. In Emerging Applications of Algebraic Geometry, vol. 149 of IMA Volumes in Mathematics and its Applications, M. Putinar and S. Sullivant (eds.), Springer, pages 157–270, 2009.
  * [50] V. Magron, M. Forets, and D. Henrion.  Semidefinite approximations of invariant measures for polynomial systems.  Discrete & Continuous Dynamical Systems - B, 24(12):6745–6770, 2019.
  * [51] V. Magron and J. Wang.  Sparse Polynomial Optimization: Theory and Practice. Series on Optimization and Its Applications, World Scientific Press, 2023.
  * [52] D. K. Molzahn, J. T. Holzer, B. C. Lesieutre, and C. L. DeMarco.  Implementation of a large-scale optimal power flow solver based on semidefinite programming.  IEEE Transactions on Power Systems, 28(4):3987–3998, 2013.
  * [53] J. Nie.  The $\mathcal{A}$-truncated $K$-moment problem.  Foundations of Computational Mathematics, 14(6):1243–1276, 2014\.
  * [54] J. Nie.  Symmetric tensor nuclear norms.  SIAM Journal on Applied Algebra and Geometry, 1(1):599–625, 2017\.
  * [55] B. Reznick.  Extremal PSD forms with few terms.  Duke Mathematical Journal, 45(2):363–374, 1978.
  * [56] C. Riener, T. Theobald, L. J. Andrén, and J. B. Lasserre.  Exploiting symmetries in SDP-relaxations for polynomial optimization.  Mathematics of Operations Research, 38(1):122–141, 2013.
  * [57] L. Salce and P. Zanardo.  Completely positive matrices and positivity of least squares solutions.  Linear Algebra and its Applications, 178:201–216, 1993.
  * [58] N. Shaked-Monderer, I. M. Bomze, F. Jarre, and W. Schachinger.  On the cp-rank and minimal cp factorizations of a completely positive matrix.  SIAM Journal on Matrix Analysis and Applications, 34(2):355–368, 2013.
  * [59] J. Sliwak, E. D. Andersen, M. F. Anjos, L. Létocart, and E. Traversi.  A clique merging algorithm to solve semidefinite relaxations of optimal power flow problems.  IEEE Transactions on Power Systems, 36(2):1641–1644, 2021.
  * [60] J. Sliwak, M. F. Anjos, L. Létocart, and E. Traversi.  A semidefinite optimization-based branch-and-bound algorithm for several reactive optimal power flow problems, arXiv:2103.13648, 2021.
  * [61] G. Tang and P. Shah.  Guaranteed tensor decomposition: A moment approach.  In F. Bach and D. Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, PMLR 37:1491–1500, 2015.
  * [62] S. A. Vavasis.  On the complexity of nonnegative matrix factorization.  SIAM Journal on Optimization, 20(3):1364–1377, 2010.
  * [63] H. Waki, S. Kim, M. Kojima, and M. Muramatsu.  Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity.  SIAM Journal on Optimization, 17(1):218–242, 2006.
  * [64] J. Wang and V. Magron.  A second order cone characterization for sums of nonnegative circuits.  In International Symposium on Symbolic and Algebraic Computation (ISSAC ’20), July 20–23, 2020, Kalamata, Greece, pages 450–457. ACM, New York, NY, USA.
  * [65] J. Wang, V. Magron, and J. B. Lasserre.  Chordal-TSSOS: A moment-SOS hierarchy that exploits term sparsity with chordal extension.  SIAM Journal on Optimization, 31(1):114–141, 2021.
  * [66] J. Wang, V. Magron, and J. B. Lasserre.  TSSOS: A moment-SOS hierarchy that exploits term sparsity.  SIAM Journal on Optimization, 31(1):30–58, 2021.
  * [67] J. Wang, V. Magron, J. B. Lasserre, and N. H. A. Mai. CS-TSSOS: Correlative and term sparsity for large-scale polynomial optimization. ACM Transactions on Mathematical Software, 48(4):1–26, 2022.
  * [68] S. Xiang and S. Xiang.  Notes on completely positive matrices.  Linear Algebra and its Applications, 271(1):273–282, 1998.
  * [69] Y. Zheng and G. Fantuzzi.  Sum-of-squares chordal decomposition of polynomial matrix inequalities.  Mathematical Programming, 197:71–108, 2023.
  * [70] Y Zheng, G. Fantuzzi and A. Papachristodoulou.  Chordal and factor-width decompositions for scalable semidefinite and polynomial optimization.  Annual Reviews in Control, 52:243–279, 2021.


# Quality of Control in the Tavis-Cummings-Hubbard Model

Raffael Düll1, Alexei Kulagin2, Wanshun Lee2,
Yuri Ozhigov2,3, Miao Huei-huei2, Zheng Keli2,
1\. Technical University Munich, Munich, Germany.
2\. Moscow State University of M.V.Lomonosov, Faculty VMK
Moscow 119992, Vorobyovy Gory, GSP-2, 2nd hum. corpus, SKI.
3\. K.A.Valiev Institute of physics and technology RAS, Moscow, Russia
Nakhimovsky Prospekt 32-a

###### Abstract

The quality of controlling a system of optical cavities in the Tavis-Cummings-
Hubbard (TCH) model is estimated with the examples of quantum gates, quantum
walks on graphs, and of the detection of singlet states. This type of control
of complex systems is important for quantum computing, for the optical
interpretation of mechanical movements, and for quantum cryptography, where
singlet states of photons and charges play an essential role. It has been
found that the main reason for the decrease of the control quality in the THC
model is due to the finite width of the atomic spectral lines, which is itself
related to the time energy uncertainty relation. This paper evaluates the
quality of a CSign-type quantum gate based on asynchronous atomic excitations
and on the optical interpretation of the motion of a free particle.

PACS: 03.65, 87.10
Keywords: Optical field, Tavis-Cummings-Hubbard model, Quantum gate,
Decoherence

##  1 Introduction

The Tavis-Cummings-Hubbard model TCH ([1], [2]) in quantum electrodynamics is
well suited to simulate the dynamics of small atomic ensembles that interact
with light. The high flexibility of the model allows us to explore both
single-photon quantum scenarios ([3], [4]) and multi-photon scenarios ([5],
[6]). In this work, we will analyze and evaluate dynamic quantum scenarios
with one or two photons, which are suitable to simulate the dynamics of a
broad range of quantum systems. In many cases, it is much easier to manipulate
the dynamics of photons than of charges. We will call this technique optical
interpretation and demonstrate it on three examples: the motion of a free
massive particle, the entanglement of a quantum gate and the optical selection
of dark states in atomic ensembles. These tasks are important to develop
optimal technology for quantum computers and for quantum cryptography
applications. The limitations of all three scenarios within the framework of
the TCH model will be discussed both with theoretical considerations and with
computer simulations.

Even standard quantum key distribution protocols that operate on individual
photons and do not use entanglement, such as B84 or BB92, provide a
fundamentally new level of secrecy compared to classical protocols. However,
using singlet states of both photons ans massive particles provide additional
advantages as in protocol AK-2017 ([7]). The key challenge in this protocol is
to prepare singlet photon (or atoms) states in the Bell state of the form
$\frac{1}{\sqrt{2}}(|01\rangle-10\rangle)$.

The preparation of singlet atomic states and of optical imitations is limited
by the technical capabilities of the physical implementation of the TCH model.
The quality of quantum gates for computations depends on the theoretical
boundaries of the TCH model itself.

The Tavis-Cummings Hamiltonian describes a system of $s$ identical two-level
atoms in an optical resonator. The transition energy between the two levels
matches with high accuracy the energy $\hbar\omega$ of a photon trapped in the
optical cavity. The interaction energies between the atoms and the field
within the cavity is denoted by $g_{j},\ j=1,2,...,s$. If
$g_{j}/\hbar\omega\ll 1$, the rotating wave approximation can be used and the
Hamiltonian takes the form:

$H_{TC}=\hbar\omega(a^{+}a+\sum\limits_{j=1}^{s}\sigma^{+}_{j}\sigma_{j})+a^{+}\bar{\sigma}+a\bar{\sigma}^{+}$
(1)

where $\bar{\sigma}=\sum\limits_{j=1}^{s}g_{j}\sigma_{j}$. The field operators
for photon creation $a^{+}$ and annihilation $a$ depend on the given cavity,
so an additional subscript $q$ is provided to identify the respective cavity.
Similarly, the atomic operators $\sigma_{j}^{+}$ and $s_{j}$ denote the
excitation and relaxation of the $j$-th atom in the cavity.

Cavities can be connected with atoms tuned to the same frequency $\omega$ to
form a network of $c$ cavities in which photons can move between two cavities
$q$ and $p$ with an amplitude $\nu_{q,p}=r_{q,p}e^{i\phi_{q,p}}$. The
parameters $r_{q,p}$ and $\phi_{q,p}$ depend linearly on the thickness of and
length of the optical fiber connecting the involved cavities. Such a system is
described by the Hamiltonian:

$H_{TCH}=\sum\limits_{q=1}^{c}H^{q}_{TC}+\sum\limits_{1\leq q<p\leq
c}\nu_{q,p}a_{q}^{+}a_{p}+\bar{\nu}_{q,p}a_{q}a_{p}^{+}$ (2)

Note that it is always possible to choose the waveguide length such that the
phase of the transition amplitudes between the cavities is zero, thus the
amplitude $\nu_{q,p}$ remains real and only its absolute value may change.

##  2 Optical interpretation of the dynamics of a free massive particle

For now, we consider a system of connected optical cavities without any atoms,
such that one photon is located in the quantum state $|\Psi(x,t)\rangle$,
where $x$ is the index of the cavity and $t$ is the current time. The task is
to select waveguides in a way that the dynamics of this state can be
approximated by a free particle in a one-dimensional space, whose coordinate
is proportional to the cavity index $x$.

This problem is the simplest example of a continuous-time quantum walk (CTQW)
[8]. This kind of dynamics allows for far-reaching generalizations, such as
the introduction of a magnetic field as suggested by [9], where the motion of
a free particle was analyzed in a two-dimensional space. We will compare the
representation of the dynamics of a one-dimensional free particle in the JCH
model without atoms with a simple analytical expression of the Feynman kernel.

All physical quantities are represented as qubits, for which they are defined
as arithmetic progressions on a finite set of values. After an appropriating
linear transformation ${\cal D}$ on the coordinate $x$, which is equivalent to
choosing a new unit of length, we can assume that the definition domain of the
wave function $[x_{0},x_{N}]$ coincides with the interval $[0,1]$. For
$N=2^{n}$ quantities, we then obtain $x_{q}=q/N$ for $\ q=0,1,...,N-1$ and we
approximate the coordinate $x$ with accuracy $1/N$ by the by the binary string
$|e_{1}e_{2}...e_{n}\rangle$ where the binary digits $e_{k}\in\\{0,1\\}$ come
from the expansion $q=\sum\limits_{k=1}^{n}2^{n-k}e_{k}$.

Any of such strings represents a basis state of a system with $n$ qubits,
which can then be used as the discrete representation of the wave function.
This representation does not have a physical dimension, since only the
transition operator ${\cal D}$ from the continuous function $\Psi(x)$ to the
qubit representation and the basis states $|j\rangle$ are dimensional. The
coordinate of the variable $a$ is $a/\sqrt{N}$ and lies within the interval
$[0,\sqrt{N}]$ and its unit could be the Planck constant in the proper system
of units. The term $c/\sqrt{N}$ has then to be associated with a momentum, and
it seems natural to assume that it belongs to the interval
$[-\sqrt{N}/2,\sqrt{N}/2]$ because a particle in $[0,\sqrt{N}]$ can move in
both directions. Therefore, the momentum has to be equal to
$\sqrt{N}(c/N-1/2)$.

The discrete form of the momentum operator is the $N$-dimensional Hermitian
operator

$p_{discr}=QFT\sqrt{N}(x_{discr}-I/2)QFT^{-1}=A^{-1}QFT\sqrt{N}x_{discr}QFT^{-1}\
A$

with the diagonal operator $A=diag(e^{\pi ia})_{a=0,1,...,N-1}$. Its
eigenvectors are given by $A^{-1}QFT^{-1}|a\rangle$ in (3) and are associated
to the respective eigenvalues $\sqrt{N}(a-1/2)$ for all $\
a=0,1/N,...,(N-1)/N$.

The dynamics of a free particle with nonzero mass is determined by the
Hamiltonian $H_{free}=p^{2}/2m$, which is expressed in the coordinate basis as
$H_{free}=QFTdiag(p^{2}/2m)QFT^{-1}$. The Quantum Fourier Transforms are given
by:

$\begin{array}[]{ll}\quad
QFT:|c\rangle&\rightarrow\frac{1}{\sqrt{N}}\sum\limits_{a=0}^{N-1}e^{-\frac{2\pi
iac}{N}}|a\rangle\\\
QFT^{-1}:|a\rangle&\rightarrow\frac{1}{\sqrt{N}}\sum\limits_{c=0}^{N-1}e^{\frac{2\pi
iac}{N}}|c\rangle\end{array}$ (3)

Each value of the coordinates $x_{q}=q/N$ for the indices $q=0,1,...,N-1$ is
associated to one optical cavity without atoms, which are connected by optical
fibers such that the transition amplitude $\nu_{q,p}$ is equal to the term
$\langle q|H_{free}|p\rangle$ dependent on the Hamiltonian $H_{free}$ of a
free particle. Such a network of $N$ connected cavities stands for the motion
of a free particle along a straight line.

The quality criterion is determined by the Feynman kernel of a free particle
of the form:

$K(x,t)=At^{-\frac{1}{2}}e^{\frac{imx^{2}}{\hbar t}};$ (4)

Its real part is shown in figure 1, А.

Figure 1: A.  General aspect of the real part of the Feynman kernel $K(x,t)$
for a free particle with nonzero mass, which is located at the initial time
$t=0$ at the initial position $x=0$. B. Real part of the amplitude of the
photon position in the optical network simulating a free particle. The phase
shift is due to the choice of the initial photon position ($x=0$) in the
cavity system at the point with coordinate $\sqrt{N}/2$. The scale along the
axes has arbitrary units.

Figure 2 shows how the phase and the transition amplitude between cavities
depend on the distance between the respective points on the straight line
along which the particles propagate. The length of the waveguide is
proportional to the real distance from the initial particle position and its
thickness is inversely proportional to this distance. The calculation was
preformed on a personal computer and the number of optical cavities did not
exceed a couple of hundreds.

As it can be seen in Figure 1, the dynamics of a massive particle can only be
reproduced with low quality, such an optical interpretation can thus only
simulate the general aspect of the dynamics. The advantage of the model is
that it uses the same basic means as quantum gates, therefore the free
movement of massive particles in space can be easily included in a quantum
computer based on the TCH model, without the need to solve the Cauchy problem
of the Schrödinger equation, as in the standard approach for quantum computers
[10].

Figure 2: Dependence of the phase and the transition amplitude on the distance
between points along the propagation line for the Hamiltonian JCH, which
simulates the propagation of a free particle

##  3 Quality of an entaglement gate on optical cavities

The entanglement gate plays a key role for quantum computing (see [11],
[12],[13]) because, together with one-qubit gates, it allows to implement
every imaginable quantum circuit as described in [14].

Photonic systems are one of the most common ways to construct such gates, as
for instance with the well-known KLM-scheme [15] or by improving [16] with
quantum teleportation [17] and by switching on atomic excitation [18]. The
results of experiments on optical cavities in [19] or [20] indicate the
practical feasibility of such gates.

The CNOT gate is a common example of entanglement gates, it acts on a pair of
qubits $|x,y\rangle$ according to the rule $CNOT:\
|x,y\rangle\rightarrow|x,x\oplus y\rangle$. In [20], it was implemented using
the vibrational degrees of freedom of an atom. It is equivalent to the well-
known universal gate CSign, which acts according to the rule $CSign:\
|x,y\rangle\rightarrow(-1)^{xy}|x,y\rangle$, in the sense that one of these
gates can be obtained from the other with the help of simply realizable one-
qubit gates: $CNOT=w_{y}\ CSign\ w_{y}$, where the Hadamard gate
$w_{y}=(\sigma_{x}+\sigma_{y})/\sqrt{2}$ is applied to the second, controlled
qubit of the operator argument. The CSign gate inverts the phase of the basis
state $|11\rangle$. The closely related operator $coCSign:\
|x,y\rangle\rightarrow(-1)^{x(y\oplus 1)}|x,y\rangle$ can easily be reduced to
the CSign operator with the help of a one-qubit gate:
$coCSign=\sigma_{x}(y)CSign\ \sigma_{x}(y)$.

The implementation of the coCSign gate is proposed in [22]. It partially uses
the idea of [21], where the states of dual-rail photons are used as qubits,
however the scheme in [22] is much simpler. The present paper discusses the
quality of this implementation of the gate.

The focus of this paper is therefore on the coCSign gate $coCSign:\
|x,y\rangle\rightarrow(-1)^{(x\oplus 1)y}|x,y\rangle$. Since it is related to
$CSign=\sigma_{x}(x)\ coCSign\ \sigma_{x}(x)$ by single qubit-gates, which are
implemented as linear optical devices, the presented estimate also applies to
the $CNOT$ gate.

##  4 Layout of the coCSign gate

We now consider a cavity with one atom inside whose energy gap between the
ground state $|0\rangle$ and the excited state $|1\rangle$ exactly matches the
energy $\hbar\omega$ of the confined photon in the cavity. Under the
assumption of the rotating wave approximation, the Hamiltonian takes the form:

$H=H_{JC}=H_{0}+H_{int};\ H_{0}=\hbar\omega
a^{+}a+\hbar\omega\sigma^{+}\sigma,\ H_{int}=g(\alpha^{+}\sigma+a\sigma^{+}),$
(5)

The basic states of the atom and of the field can be written as
$|n\rangle_{ph}|m\rangle_{at}$, where $n=0,1,2,...$ is the number of photons
in the cavity and $m=0,1$ is the level of atomic excitation. In our setup,
there are in total 2 photons, so $n$ can take the values $n=0,1,2$. To
implement the gate, a term $H_{jump}$ is added to the Hamiltonian above to
represent the transition of a photon from the cavity $i$ to $j$ and vice
versa.

$H_{jump}=\nu(a_{i}a^{+}_{j}+a_{j}a^{+}_{i}),$ (6)

Let $\tau_{1}=\pi\hbar/g,\ \tau_{2}=\pi\hbar/g\sqrt{2}$ be the periods of the
Rabi oscillations for the total energies $\hbar\omega$ and $2\hbar\omega$. The
operator $U_{t}=e^{-\frac{i}{\hbar}Ht}$, that describes the evolution of the
system for most of the time, depends on the total energy in the cavity. If
this energy is equal to $\hbar\omega$ in the basis states
$|\phi_{0}\rangle=|1\rangle_{ph}|0\rangle_{at}$ and
$|\phi_{1}\rangle=|0\rangle_{ph}|1\rangle_{at}$, we obtain:

$U_{\tau_{1}/2}=-i\sigma_{x},\ U_{\tau_{1}}=-I,\ U_{2\tau_{1}}=I,$ (7)

In this equation, $\sigma_{x}$ is the Pauli matrix. A similar relation holds
for $\tau_{2}$ for the total energy of the cavity $2\hbar\omega$. The phase of
the state is relevant here, since states with different numbers of photons in
the cavities may be in superposition.

The transition of a photon from a cavity $i$ to another cavity $j$ is carried
out by the simultaneous switching on of the Pockels cells, which translates
into the addition of the term $H_{jump}$ to the main interaction $H_{int}$. In
the absence of atoms, it leads to the same dynamics as in Rabi oscillations,
but with a period $\tau_{jump}=\pi\hbar/\nu$. In the ideal case, we can assume
that $\nu\gg g$, which means that it is possible to move a photon between two
cavities without being affected by the atom at all. The paper [21] describes
the difficulty to implement this behavior in an experiment, however, we will
consider this problem as a technical issue.

Because of the incommensurability of the Rabi oscillation periods $\tau_{1}$
and $\tau_{2}$, we can choose two natural numbers $n_{1}$ and $n_{2}$ such
that the following approximate equality holds.

$2n_{2}\tau_{2}\approx 2n_{1}\tau_{1}+\frac{\tau_{1}}{2}$ (8)

The nonlinear shift required to implement the $coCSing$ gate is obtained with
exactly this difference in the period.

The basis state $|0\rangle$ is created in our model as the state of the
optical cavity $|0\rangle_{ph}|1\rangle_{at}$, and the basis state $|1\rangle$
as $|1\rangle_{ph}|0\rangle_{at}$. Thus, the two-qubit state $|01\rangle$ on
which the phase shift should act, represents the state
$|01\rangle_{ph}|10\rangle_{at}$, where the first photonic qubit belongs to
the cavity $x$ and the second to the cavity $y$. It is to be noted that after
a time $\tau_{1}/2$, the null and zero states change place with a phase
advance of $-\pi/2$.

The sequence of operations to implement the $coCSign$ gate is shown in Figure
3 and the involved cavities are shown in Figure 4. In the beginning, an almost
instantaneous photon exchange between the auxiliary and the $x$ cavity is
induced. After a delay $\tau_{1}/2$, the same photon exchange happens between
the auxiliary and the $y$ cavities. Then, the system evolves without any
exchanges for the time $2n_{2}\tau_{2}$, before the same sequence of photon
exchanges occurs again between the auxiliary and respectively the $x$ and $y$
cavities with a time gap $\tau_{1}/2$ in between. With this choice of photon
displacement times, it follows that at these moments there are exactly one or
zero photons in a cavity. Therefore, switching on the Pockels cells during a
short time $\delta\tau=\pi\hbar/2\nu\ll\tau_{1}$ results in exactly the
desired photon movements.

Figure 3: Operation sequence of the implementation of the coCSign gate:
$|x,y\rangle\rightarrow(-1)^{(x\oplus 1)y}|x,y\rangle$ on asynchronous atomic
excitations in optical cavities, split into 7 time segments, and short photon
transition times $\delta\tau=\tau_{jump}/2\ll\tau_{1(2)}$. The expected
behavior of the coCSign gate is observed after an additional $\tau_{1}/2$ at
the end of diagram Figure 4: The three cavities ($x,y$ and auxiliary) involved
in the coCSign gate

In the work [22], it was proved that this scheme ideally gives an
implementation of the coCSign gate. In the article [23], it was calculated
that even if the numbers $n_{1}$ and $n_{2}$ are below 100, and if all other
decoherence factors are ignored, the gate accuracy exceeds 90%.

The interaction force between the atom and the field in the cavity takes the
form:

$g=\sqrt{\hbar\omega/V}d\ E(x),$ (9)

where $V$ is the effective volume of the cavity, $d$ is the dipole moment of
the transition between the ground an the perturbed states and $E(x)$ describes
the spatial arrangement of the atom in the cavity and is equal to
$E(x)=sin(\pi x/L)$ with $L$ the length of the cavity. To ensure the
confinement of the photon in the cavity, $L$ has to be chosen such that
$L=n\lambda/2$ is a multiple of the photon wavelength $\lambda$. In
experiments, $n=1$ is often chosen to decrease the effective volume of the
cavity, which makes it possible to obtain dozens of Rabi oscillations (see,
for example [19]).

The physical limitation on the gate quality follows from the energy-time
uncertainty relation for photons. The parameter $\nu$ for the photon
transition from one cavity to another cannot be chosen too large because a
short transition time automatically means a large uncertainty in its energy. A
large error in the energy results in an altered wavelength and if the
wavelength becomes too different from twice the cavity length, the photon can
escape the cavity.

For instance, the photon frequency in experiments with the rubidium atom is
approximately $10^{10}\ sec^{-1}$ and the upper bound for the possible
uncertainty in frequency is $10^{9}\ sec^{-1}$ (the real bound is however much
lower). Taking into account the uncertainty relation $\delta\omega\ \delta
t\approx 1$, a lower estimate for the photon transition time window could be
given by $\delta\tau\approx 10^{-9}\ sec$. Since the period of the Rabi
oscillations is approximately $\tau\approx 10^{-6}\ sec$, we obtain the bounds
$10^{-9}\leq\delta\tau\ll 10^{-6}$ for the transition time, which means that a
single execution of the gate will induce an error of a least than $10^{-3}$.

Figure 5: Left side: exchange of a photon between the auxiliary and $x$
cavities for the ground-state atoms, with the parameter
$\alpha=\tau_{1}/2\sqrt{2\pi}\sigma$. Right side: Evolution of the segment
with length $2n_{2}\tau_{2}$.

Let us suppose, that the intensity $\nu$ of any photon transfer between
cavities from equation (6) is described by the Gaussian

$\nu(t)=\alpha e^{-\frac{(t-t_{id})^{2}}{2\sigma^{2}}}.$ (10)

The parameter $\sigma$ needs to be chosen $\sigma\ll\tau_{1,2}$ to fulfill the
aforementioned condition for the gate implementation. In this case, the
transition of a photon from the cavity $x$ to the auxiliary cavity and back is
regular for the ground states of an atom and an example of it is shown in the
left part of Figure 5. The evolution of the state in the central segment of
duration $2n_{2}\tau_{2}$ is shown in the right part.

To assess the quality of the gate, we consider the initial state
$|\psi_{0}\rangle=\frac{1}{2}(|00\rangle+|01\rangle+|10\rangle+|11\rangle)$.
First, we assume that the maximum intensity $\alpha$ has a fixed value. In
this case, the accuracy can be evaluated by the trace and modular distances
between the ideal state $\rho_{id}=|\Psi_{id}\rangle\langle\Psi_{id}|$ and the
real state $\rho=|\Psi\rangle\langle\Psi|$ of the system. We define the
differences $\delta\rho=\rho-\rho_{id}$ and
$\delta\Psi=|\Psi\rangle-|\Psi_{id}\rangle$ to estimate the errors:

$d_{tr}=tr\sqrt{\delta\rho^{+}\delta\rho}\qquad\
d_{mod}=\delta\Psi^{+}\delta\Psi.$ (11)

The results of the numerical simulations is shown in Figure 6.

Figure 6: Trace and modular distances between the ideal and the decoherent
results of the coCSign gate or the initial state $|\psi_{0}\rangle$ and for
various combinations of $n_{1}$ and $n_{2}$.

Figure 7 shows the modular distance between the ideal and the simulated,
decoherent gate for a fixed $\sigma$ and variable $\alpha$ and $\sigma$. The
observed periodicity follows from the possibility of multiple photon
transitions between the cavities during one activation of the Pockels cells
for increasing $\alpha$.

Figure 7: A. Modular distance between the ideal and decoherent states after
the application of the coCSign gate on different initial states in function of
$\alpha$ for the fixed values $\sigma=0.5,g=10^{-3},\hbar=1,\omega=1$. B.
Modular error for variable $\alpha$ and $\sigma$.

##  5 Dark states of two-level atomic ensembles

The complexity of the model can be further increased if we consider several
atoms in one cavity. A pure quantum effect then occurs: dark states of atomic
ensembles.

In a dark state, an atom cannot emit any light, although its total energy is
not zero. For two-level atoms, such states are obtained as the linear
combinations of tensor products of singlets. With the RWA approximation, the
ground states of atom parts are also allowed in the product. The Bell state
$|s\rangle=\frac{1}{\sqrt{2}}(|01\rangle-|10\rangle)$ has the remarkable
property that for any unitary operator $U$, its tensor square does not change
the singlet state: $U^{\otimes 2}|s\rangle=|s\rangle$. The quantum protocol to
distribute the secret key AK-2017 is based on this singlet property [7].

The preparation of such singlet states of photons is therefore a key issue for
quantum cryptographic protocols, and the quality assessment of this
preparation plays a practical role. Singlets can be generalized to ensembles
of multi-level atoms by taking the determinant of their base states. For
instance, the multi-singlet for a three-level atom takes the form
$|D_{3}\rangle=\frac{1}{\sqrt{6}}\sum\limits_{\pi\in
S_{3}}(-1)^{\sigma(\pi)}|\pi(1)\pi(2)\pi(3)\rangle$, where the sum is taken
over all permutations of three atoms.

The simplest way to obtain dark singlet atomic states is through optical
selection. One photon is launched into a cavity with an atomic ensemble and a
vacuum state of the field. Then, its escape time from the cavity is measured,
and after many repetitions, the average photon emission time can be determined
very accurately by virtue of the central limit theorem. If the initial state
of the atomic ensemble is a superposition of singlet products (for a diatomic
ensemble, of a single singlet), the photon escapes much faster because it
cannot interact with the atoms. On the other hand, if the ensemble is not in a
singlet state, there is interaction between the photon and the atoms, which
causes a delay in its exit.

The optical selection makes it also possible to estimate the degree of
darkness of an atomic ensemble, as illustrated by Figure 8.

Figure 8: Distribution density of the photon emission time from the cavity for
one dark and two different types of light states of an atomic ensemble with
three $Rb^{85}$ atoms (Figure taken from [24]).

In [24], it was shown that the quality of optical selection is reduced due to
the quality of the photon detectors, where the error of any kind is less then
3%. The accuracy of the selection itself is much higher, since this systematic
error does not affect the result obtained by the comparison of two alternative
hypotheses of dark or light states in the ensemble. The random error decreases
rapidly with the increase in the number of experiments.

##  6 Conclusion

The quality assessment of quantum control with the Tavis-cummings-Hubbard
model was carried out on the standard problem of a free quantum particle, on
the main controlled two-qubit gate in quantum computing and on the problem of
obtaining dark states of two-level atoms, which is important for quantum
cryptography. Decoherence occurs because of the incommensurability of the Rabi
oscillation periods and because of the broadening of spectral lines for an
atom in a cavity due to the energy-time uncertainty relation. These effects,
in combination with the finite response time of the Pockels cell, may
significantly reduce the accuracy of the gate. The latter factor seems to be
only a technical issue, and was not further investigated in this work. The
incommensurability induces an error of the order of one percent, which is a
very good result for the gate. However, this is not a definite assessment of
the feasible gate quality. The broadening of spectral lines was taken into
account too roughly to draw certain conclusions about its contribution to the
error and needs a more detailed analysis. Together with the technical
limitations on the response rate of the Pockels cell, it imposes significant
restrictions on the achievable accuracy of the main entanglement gate, and
this regardless of its specific shape (CNOT, CSign or coCSign), and yields a
serious obstacle to scale the Feynman scheme to a quantum computed based on
gates.

It was also shown that for dark states of atomic ensembles, quantum operations
can be achieved with a significantly higher accuracy, even with relatively
imprecise photon detectors. This is facilitated by statistical methods to
choose between two alternative hypotheses, however these are only measurable
operations. The use of statistical methods in quantum computing means elements
of classical parallelization, similar to the KLM scheme [15], which should be
part of any implementation of a quantum computer.

Dark states can theoretically persist for a very long time, since the atomic
ensemble in such a state does not interact with the resonator field.
Therefore, the development of the concept of quantum computing from the
traditional Feynman scheme to the use of dark states promises new weapons in
the fight against decoherence.

The analyzed examples demonstrate the great importance of the Tavis-Cummings-
Hubbard model in theoretical estimates of the quality of quantum operations.
This model could be the basis for the creation of an operating system for a
quantum computer, where the central part of the computing technology bases on
the first principles of quantum computing.

##  7 Acknowledgments

This work was supported by the Russian Foundation for Basic Research, grant
a-18-01-00695.

## References

  * [1] Jaynes E.T., Cummings F.W., Comparison of quantum and semiclassical radiation theories with application to the beam maser, Proc. IEEE 51 (1): 89–109, (1963). doi:10.1109/PROC.1963.1664
  * [2] Michael Thomas Tavis, A Study of an N Molecule Quantized-Radiation-Field Hamiltonian, Dissertation, https://arxiv.org/abs/1206.0078.
  * [3] M. I. Makin, Jared H. Cole, Charles D. Hill, Andrew D. Greentree, and Lloyd C. L. Hollenberg, Time evolution of the one-dimensional Jaynes-Cummings-Hubbard Hamiltonian, Phys. Rev. A 80, 043842 – Published 30 October 2009.
  * [4] Guilherme M. A. Almeida and Andre M. C. Souza, Quantum transport with coupled cavities on an Apollonian network, Phys. Rev. A 87, 033804 – Published 8 March 2013.
  * [5] Shifeng Cui, F. Hébert, B. Grémaud, V. G. Rousseau, Wenan Guo, and G. G. Batrouni, Two-photon Rabi-Hubbard and Jaynes-Cummings-Hubbard models: Photon-pair superradiance, Mott insulator, and normal phases, Phys. Rev. A 100, 033608 – Published 9 September 2019
  * [6] L. Garbe, I. L. Egusquiza, E. Solano, C. Ciuti, T. Coudreau, P. Milman, and S. Felicetti, Superradiant phase transition in the ultrastrong-coupling regime of the two-photon Dicke model, Phys. Rev. A 95, 053854 – Published 23 May 2017
  * [7] Матвеев Е.А., Игошина С.Е., Карманов А.А., Квантовый фазовый переход как основа для практической реализации ATF-технологии связи, Программная инженерия, 2019, том 20, N 7-8, стр. 305-310.
  * [8] A. Ambainis, Quantum walks and their algorithmic applications, International Journal of Quantum Information, 1:507-518, 2003.
  * [9] Luca Razzoli, Matteo G.A. Paris, Paolo Bordone, Continuous-time quantum walks on planar lattices and the role of the magnetic field, Phys. Rev. A 101, 032336 (2020).
  * [10] C.Zalka, Simulating quantum systems on a quantum computer, Proceedings of The Royal Society A 454(1969):313-322, January 1998.
  * [11] Richard P. Feynman, Simulating Physics with Computers, International Journal of Theoretical Physics, VoL 21, Nos. 6/7, 1982, pp. 467-488.
  * [12] L.Grover, A fast quantum mechanical algorithm for database search, Proceedings, 28th Annual ACM Symposium on the Theory of Computing (STOC), May 1996, pages 212-219. Proceedings, Melville, NY, 2006, vol. 810, electronic version: xxx.lanl.gov, quant-ph/0610052.
  * [13] P.Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, SIAM J.Sci.Statist.Comput. 26 (1997) 1484.
  * [14] A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J. Smolin, H. Weinfurter, Elementary gates for quantum computation Phys.Rev. A52 (1995) 3457.
  * [15] Knill, E., Laflamme, R., Milburn, G. J. (2001). "A scheme for efficient quantum computation with linear optics". Nature. Nature Publishing Group. 409 (6816): 46–52
  * [16] Gottesman, D., Chuang, I. L. (1999-11-25). "Demonstrating the viability of universal quantum computation using teleportation and single-qubit operations". Nature. 402 (6760): 390–393
  * [17] Bennett, Charles H.; Brassard, Gilles; Crépeau, Claude; Jozsa, Richard; Peres, Asher; Wootters, William K. (1993-03-29). "Teleporting an unknown quantum state via dual classical and Einstein-Podolsky-Rosen channels". Physical Review Letters. 70 (13): 1895–1899.
  * [18] Popescu, S., Knill-Laflamme-Milburn Quantum Computation with Bosonic Atoms, PRL 99, 130503 (2007).
  * [19] G. Rempe, H. Walther, and N. Klein. Observation of quantum collapse and revival in a one-atom maser (англ.) // Phys. Rev. Lett. : journal. — 1987. — Vol. 58, no. 4. — P. 353.
  * [20] C. Monroe, D. M. Meekhof, B. E. King, W. M. Itano, and D. J. Wineland, Demonstration of a Fundamental Quantum Logic Gate, Phys. Rev. Lett. 75, 4714 (1995).
  * [21] Azuma H., Quantum computation with the Jaynes-Cummings model, Prog. Theor. Phys. 126, 369-385 (2011).
  * [22] Y.I.Ozhigov, Quantum gates on asynchronous atomic excitations, Quantum Electronics (2020), 50 (10):947-950, DOI:10.1070/QEL17320.
  * [23] Ladunov, V., Ozhigov, Y., Skovoroda, N., Computer simulation of quantum effects in Tavis-Cummings model and its applications, SPIE Proceedings, vol. 10224, International Conference on Micro- and Nano-Electronics 2016; 102242X (2016) https://doi.org/10.1117/12.2267190
  * [24] A.V.Kulagin, Y.I.Ozhigov, Optical selection of dark states (accepted for publication in the journal ’’Computational Mathematics and modeling’’)