# Entanglement in the Quantum Game of Life

Peter-Maximilian Ney  Theoretical Physics, Saarland University, D-66123
Saarbrücken, Germany  Simone Notarnicola  Dipartimento di Fisica e Astronomia
“G. Galilei”, via Marzolo 8, I-35131, Padova, Italy  Padua Quantum
Technologies Research Center, Università degli Studi di Padova.  INFN, Sezione
di Padova, via Marzolo 8, I-35131, Padova, Italy  Simone Montangero
Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8, I-35131,
Padova, Italy  Padua Quantum Technologies Research Center, Università degli
Studi di Padova.  INFN, Sezione di Padova, via Marzolo 8, I-35131, Padova,
Italy  Giovanna Morigi  Theoretical Physics, Saarland University, D-66123
Saarbrücken, Germany

###### Abstract

We investigate the quantum dynamics of a spin chain that implements a quantum
analog of Conway’s game of life. We solve the time-dependent Schrödinger
equation starting with initial separable states and analyse the evolution of
quantum correlations across the lattice. We report examples of evolutions
leading to all-entangled chains and/or to time oscillating entangling
structures and characterize them by means of entanglement and network
measures. The quantum patterns result to be quite different from the classical
ones, even in the dynamics of local observables. A peculiar instance is a
structure behaving as the quantum analog of a blinker, but that has no
classical counterpart.

##  I Introduction

The Game of Life (GoL) is a celebrated cellular automaton originally proposed
by Conway and since then used as a paradigm to study selforganization Gardner
(1970). Its dynamics takes place on a network according to simple rules, it
depends on the initial state and it is Turing-complete Schulman and Seiden
(1978); Rendell (2002). The implementation of the GoL in the quantum world has
been object of increasing interest as a possible architecture for quantum
computation P.Arrighi , for variational quantum optimization protocols and
quantum state engineering Brennen and Williams (2003); Wintermantel et al.
(2020). Moreover, it can be used to investigate the role of quantum mechanics
in the emergence of physical complexity Wolfram (1983); Hillberry et al.
(2020). A central issue of quantum cellular automata is the characterization
of the spreading of quantum information and of entanglement. Recent works
analysed the spreading of entanglement generated by quantum cellular automata
Brennen and Williams (2003); Hillberry et al. (2020), in Ref. Piroli and Cirac
(2020) an explicit relation is identified between quantum cellular automata
and tensor networks, showing that the entanglement entropy satisfies an area
law bound that can be generated at each dynamical step of the quantum cellular
automata.

In this work we analyse the buildup of quantum correlations in the quantum
game of life introduced in Ref. Bleh et al. (2012): a Hamiltonian on a one-
dimensional lattice simulating a reversible cellular automaton implementing
the classical $F_{12}$-rule, which swaps the state of a given site whenever
either two or three sites of the nearest or next-nearest neihbours are alive
Wolfram (1983). We consider initial separable states with a well defined
classical analog and determine their evolution by integrating the Schrödinger
equation, both by methods based on ordinary differential equations as well as
by means of tensor networks simulations Schollwöck (2011); Montangero (2018).
We compare the classical and the quantum dynamics obtained by evolving the
same initial uncorrelated state and analyse the corresponding quantum
correlations by means of quantum-information theoretical measures. We show
that the emerging quantum patterns exhibit qualitatively different
characteristics from their classical counterpart in the behaviour of both
local and non-local observables. One peculiar structure is the quantum
blinker, namely a quantum recurrent structure that in the quantum GOL has no
classical counterpart.

This manuscript is organized as follows. In Sec. II we introduce the model and
the observables we use in order to characterize the dynamics. In Sec. III we
present the results of the numerical integration of the Schrödinger equation
for different initial states. An analytical model for the semiclassical analog
of classical recurrent structures is derived in Sec. IV. Finally, in Sec. V
the conclusions are drawn and outlooks are discussed.

##  II A quantum game of life

We consider a chain of $L$ quantum spins $1/2$. We denote the states of the
spin $j$ by $\\{\ket{0}_{j},\ket{1}_{j}\\}$, which are here chosen to be the
eigenstates of the Pauli operator $\hat{\sigma}_{z}^{j}$. We denote by
$\hat{b}_{j}$ and $\hat{b}_{j}^{\dagger}$ the lowering and rising operators,
which act on the basis states according to the relations
$\hat{b}_{j}\ket{0}=0$, $\hat{b}_{j}^{\dagger}\ket{0}=\ket{1}$, and
$\hat{b}_{j}^{\dagger}\ket{1}=0$, $\hat{b}_{j}\ket{1}=\ket{0}$, and such that
$\hat{\sigma}_{z}^{j}=\hat{b}_{j}\hat{b}_{j}^{\dagger}-\hat{b}_{j}^{\dagger}\hat{b}_{j}$.
We note that an eigenstate of the operator
$\sum_{j=1}^{L}\hat{\sigma}_{z}^{j}$ is a sequence of $L$ bits.
In what follows, $\ket{0}_{j}$ symbolizes the dead and $\ket{1}_{j}$ the alive
cell, respectively. Let $\ket{\psi}_{0}$ be the initial state of the chain.
Its evolution is governed by the Schrödinger equation

${\rm i}\partial_{t}\ket{\psi}_{t}=\hat{H}\ket{\psi}_{t}\,,$ (1)

where $t$ is a continuous dimensionless variable ($\hbar=1$). The dynamics is
unitary, the state at time $t$ is $\ket{\psi}_{t}=\hat{U}(t)\ket{\psi}_{0}$
where $\hat{U}(t)=\exp(-{\rm i}\hat{H}t)$ is the time evolution operator.
Following Ref. Bleh et al. (2012), the model Hamiltonian reads

$\displaystyle\hat{H}\coloneqq\sum_{i=3}^{L-2}\hat{\mathcal{S}}_{i}\quantity(\hat{\mathcal{N}}_{i}^{(2)}+\hat{\mathcal{N}}_{i}^{(3)})\,.$
(2)

The evolution it describes realises a continuous-time, quantum mechanical
analog of the classical $F_{12}$-rule, where the state of a cell is swapped if
and only if two or three cells among the nearest and next-nearest neighbours
are alive (namely, in state $\ket{1}$). Indeed,
$\hat{\mathcal{S}}_{i}=\hat{b}_{i}+\hat{b}_{i}^{\dagger}$ and the operators
$\hat{\mathcal{N}}_{i}^{(2),(3)}$ are defined to be different from zero over
the set of states where only two (three) of the four nearest and next-nearest
neighboring cells about $i$ are alive. They are the sums of tensor products of
projection operators $\hat{\bar{n}}_{j}\equiv\ket{0}_{j}\bra{0}$ and
$\hat{n}_{\ell}\equiv\ket{1}_{\ell}\bra{1}$ for the nearest and next-nearest
neighbouring sites $j$ and $\ell$ of site $i$ with eigenvalues is $0,1$. They
take the form

$\displaystyle\hat{\mathcal{N}}_{i}^{(2)}=\phantom{+}$
$\displaystyle\hat{\bar{n}}_{i-2}\hat{\bar{n}}_{i-1}\hat{n}_{i+1}\hat{n}_{i+2}+\hat{\bar{n}}_{i-2}\hat{n}_{i-1}\hat{\bar{n}}_{i+1}\hat{n}_{i+2}$
$\displaystyle+$
$\displaystyle\hat{\bar{n}}_{i-2}\hat{n}_{i-1}\hat{n}_{i+1}\hat{\bar{n}}_{i+2}+\hat{n}_{i-2}\hat{\bar{n}}_{i-1}\hat{\bar{n}}_{i+1}\hat{n}_{i+2}$
$\displaystyle+$
$\displaystyle\hat{n}_{i-2}\hat{\bar{n}}_{i-1}\hat{n}_{i+1}\hat{\bar{n}}_{i+2}+\hat{n}_{i-2}\hat{n}_{i-1}\hat{\bar{n}}_{i+1}\hat{\bar{n}}_{i+2}$
$\displaystyle\hat{\mathcal{N}}_{i}^{(3)}=\phantom{+}$
$\displaystyle\hat{\bar{n}}_{i-2}\hat{n}_{i-1}\hat{n}_{i+1}\hat{n}_{i+2}+\hat{n}_{i-2}\hat{\bar{n}}_{i-1}\hat{n}_{i+1}\hat{n}_{i+2}$
$\displaystyle+$
$\displaystyle\hat{n}_{i-2}\hat{n}_{i-1}\hat{\bar{n}}_{i+1}\hat{n}_{i+2}+\hat{n}_{i-2}\hat{n}_{i-1}\hat{n}_{i+1}\hat{\bar{n}}_{i+2}\,,$

their eigenvalues are $0,1$.

We integrate numerically the dynamics for chains of finite size and for
sufficiently long times. The resulting states and configurations are
significantly affected by the boundary conditions. Here, we implement open
boundary conditions, such that the state of cells $j=1,\,2$ and $j=L-1,\,L$ is
constant, as visible from Hamiltonian (2). The initial states we consider are
separable and are eigenstates of the operator
$\sum_{j=1}^{L}\hat{\sigma}_{z}^{j}$, that we denote below as Fock states and
have a well defined classical analog.

We remark that Eq. (2) describes a unitary, and thus reversible, continuous-
time dynamics. The classical GOL, instead, evolves according to a discrete-
time, irreversible map fixed by a rule (the $F_{12}-$rule in our case). Here,
in order to compare the dynamics of the two models we compare the unitary
evolution with the classical analog. The classical dynamics can also be
derived by considering a modified stroboscopic dynamics of the quantum system
in which we set the duration of the time intervals to $\pi/2$. In fact, this
is the time needed for a single spin in the state $\ket{0}$ or $\ket{1}$ to be
completely swapped under the action of the operator $\hat{\mathcal{S}}_{i}$.
At the beginning of each time step, moreover, we iteratively measure and fix
the values of the projectors $\hat{\mathcal{N}}_{i}^{(2),(3)}$ and
consequently let all the sites freely evolve or remain in their initial state.
By doing so, if the initial state is a Fock state, then the state at each step
is a Fock state, and the sequence is the same as the one generated by the
classical $F_{12}$-rule.

###  II.1 Local observables

We first discuss local observables, corresponding to measurements of operators
$\hat{n}_{i}$ and $\hat{\bar{n}}_{i}$. The evolution of the distribution of
the alive cells is given by the expectation value

$\displaystyle n_{i}(t)=\expectationvalue{\hat{n}_{i}}\,,$ (3)

and by its discretized description,

$\displaystyle\mathcal{D}_{i}(t)$ $\displaystyle\coloneqq\left\lceil
n_{i}(t)-0.5\right\rceil\,,$ (4)

which is either $1$ or $0$ depending on whether $n_{i}(t)>0.5$ or
$n_{i}(t)\leq 0.5$, respectively. According to this discretized description,
the fraction of alive sites (density) is described by the quantity

$\displaystyle\rho(t)$ $\displaystyle\coloneqq
L^{-1}\sum_{i}\mathcal{D}_{i}(t)\,.$ (5)

The formation of clusters of alive cells with size $\ell$ (with $0\leq\ell\leq
L$) is revealed by the clustering function

$\displaystyle C(\ell,t)$
$\displaystyle\coloneqq\smash{\sum_{i=1}^{L-\ell+1}}\,\hphantom{\cdot}\Theta_{P}\quantity(\mathcal{D}_{i-1}(t)=0)$
(6)
$\displaystyle\hphantom{\coloneqq\sum_{i=1}^{L-\ell+1}}\cdot\Theta_{P}(\mathcal{D}_{j}(t)=1\,\forall\,j\in\quantity{i,\ldots,i+\ell-1})$
$\displaystyle\hphantom{\coloneqq\sum_{i=1}^{L-\ell+1}}\cdot\Theta_{P}\quantity(\mathcal{D}_{i+\ell}(t)=0)\,,$

where here $\Theta_{P}$ denotes the Heaviside function for propositions
$\mathcal{A}$, such that $\Theta_{P}(\mathcal{A})=0$ if $\mathcal{A}$ is false
while $\Theta_{P}(\mathcal{A})=1$ otherwise. It is useful to consider also the
clustering function for dead cells:

$\displaystyle\bar{C}(\ell,t)$
$\displaystyle\coloneqq\smash{\sum_{i=1}^{L-\ell+1}}\,\hphantom{\cdot}\Theta_{P}\quantity(\mathcal{D}_{i-1}(t)=1)$
(7)
$\displaystyle\hphantom{\coloneqq\sum_{i=1}^{L-\ell+1}}\cdot\Theta_{P}(\mathcal{D}_{j}(t)=0\,\forall\,j\in\quantity{i,\ldots,i+\ell-1})$
$\displaystyle\hphantom{\coloneqq\sum_{i=1}^{L-\ell+1}}\cdot\Theta_{P}\quantity(\mathcal{D}_{i+\ell}(t)=1)\,.$

The number of different (living) cluster sizes present in the system at time
$t$ is quantified by the diversity Schulman and Seiden (1978):

$\displaystyle\Delta(t)$
$\displaystyle\coloneqq\sum_{\ell=1}^{L}\Theta_{0}(C(\ell,t))\,,$ (8)

that is defined using the Heaviside function $\Theta_{0}(x)$ Schulman and
Seiden (1978). Furthermore, we will also use the so-called improved diversity,
which counts over the average of dead and alive cells,

$\displaystyle\Delta^{\text{impr.}}(t)$
$\displaystyle\coloneqq\frac{1}{2}\sum_{\ell=1}^{L}\Theta_{0}(C(\ell,t))+\Theta_{0}(\bar{C}(\ell,t))\,.$
(9)

The diversity is a measure of complexity, in the sense that its increase is
associated with an increase of the complexity Grassberger (1986); Gros (2015).

###  II.2 Quantum correlations

Quantum (non-classical) correlations are studied by means of the entanglement
entropy $S(t)$ obtained by partitioning the system into a subset $A$ of cells
and the rest of the lattice, which we denote by $B$ Hillberry et al. (2020):

$\displaystyle S(t)$
$\displaystyle=-\Tr\\{\hat{\rho}_{A}(t)\log_{2}\quantity(\hat{\rho}_{A}(t))\\}\,,$
(10)

where the density matrix of the subset $A$ is obtained by tracing out the
degrees of freedom of $B$, namely,
$\hat{\rho}_{A}(t)=\Tr_{B}\\{\ket{\psi}_{t}\bra{\psi}\\}$, and $\Tr_{B}$
denotes the trace over a basis of subspace $B$.

We characterize quantum correlations by means of different partitions of the
lattice: We determine the single-point entanglement entropy $S_{i}(t)$ when
$A$ consists of the single site $i$. The two-point entanglement entropy
$S_{i,j}(t)$ is found when the subsystem $A$ consists of two sites $i,\,j$,
which are not necessarily neighbours Vargas (2016); Valdez et al. (2017).
These two quantities allow us to determine the mutual information,

$\displaystyle\mathcal{I}_{i,j}(t)$
$\displaystyle\coloneqq\frac{1}{2}\quantity(S_{i}(t)+S_{j}(t)-S_{i,j})(t)\mbox{\quad
with\quad}\mathcal{I}_{i,i}\equiv 0\,,$ (11)

which quantifies correlations between two sites Cover and Thomas (2006). The
matrix $\mathcal{I}$ can be interpreted as the adjacency matrix of a weighted
network. In order to quantify the entanglement between a pair of sites, we
also use the concurrence between sites $i$ and $j$ Horodecki et al. (2009).
Given the two-sites density matrix $\rho$ and its complex conjugate $\rho^{*}$
expressed in the $\sigma_{z}$ basis, one defines the matrix
$\tilde{\rho}=(\sigma_{y}\otimes\sigma_{y})\rho^{*}(\sigma_{y}\otimes\sigma_{y})$,
the concurrence is

$\mathcal{C}_{ij}:=\mathrm{max}\\{0,\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4}\\}\,,$
(12)

where $\lambda_{\mu}$ are the eigenvalues of the matrix
$R=\sqrt{\sqrt{\rho}\tilde{\rho}\sqrt{\rho}}$ put in descending order. In our
analysis we will use the averaged concurrence $\mathcal{C}(d)$ obtained by
averaging over all the pairs of sites with distance $d$. Moreover, following
Ref. Hillberry et al. (2020), one can determine typical network measures, such
as the density of the network, which gives the fraction of existing links
(correlations) between two nodes (cells) over all possible links:

$\displaystyle D$
$\displaystyle\coloneqq\frac{1}{L(L-1)}\sum_{i,j=1}^{L}\mathcal{I}_{i,j}\,.$
(13)

Further insight can be extracted from the average disparity $Y$ and from the
clustering coefficient $C$: The first quantifies the similarity of the network
to a backbone Serrano et al. (2009); Hillberry et al. (2020), while the second
is a measure of the transitivity in a network Newman (2003). They are defined
respectively as

$\displaystyle Y$ $\displaystyle\coloneqq
L^{-1}\sum_{i=1}^{L}\frac{1}{\quantity(\sum_{k=1}^{L}\mathcal{I}_{i,k})^{2}}\sum_{j=1}^{L}\mathcal{I}_{i,j}^{2}\,,$
(14) $\displaystyle C$
$\displaystyle\coloneqq\frac{\Tr\quantity(\mathcal{I}^{3})}{\sum_{i=1}^{L}\sum_{j=1,\,j\neq
i}^{L}\quantity[\mathcal{I}^{2}]_{i,j}}\,.$ (15)

We note that density, clustering, and disparity are here calculated using the
mutual information, they thus provide insight into the behaviour of
entanglement across the system, namely, a network whose nodes are the spins
and whose edges are dynamically determined by Eq. (1). Finally, entanglement
across the lattice is quantified by means of the bond entropy $S_{\rm
Bond}(t,j)$, which is found from Eq. (10) for a bipartition of the system at
the bond between sites $j$ and $j+1$.

##  III Quantum dynamics

In this section we discuss the dynamics of three types of initial states: (i)
spatially localized states, where the initial clusters have unit size, (ii) a
cluster of a dozen cells, and (iii) an arbitrary sequence of 0 and 1 with
varying density of alive cells, which we denote by random Fock state. The
dynamics of these three kinds of states have been first analysed in Ref. Bleh
et al. (2012). Here, we extend those studies by reporting the corresponding
dynamics of quantum correlations and thereby elucidating the emergence of
complexity with respect to the classical game of life. Let us remind that all
the initial states we consider are separable Fock states. Therefore, they can
be represented as a sequence of 0 and 1 and have a well defined classical
analog.

In the following, for lattices up to 20 sites we numerically integrate the
Schrödinger equation, Eq. (1), by means of fourth-order ordinary differential
equations solver (ODE). For larger lattice sizes we resort to Matrix Product
States (MPS) Schollwöck (2011); Silvi et al. (2019) and the Time Dependent
Variational Principle (TDVP) Haegeman et al. (2016) from the library _OpenMPS_
OpenMPS Documentation contributors (2020); Jaschke et al. (2018); Wall and
Carr (2012).

###  III.1 Quantum blinkers

Figure 1: (color online) Time evolution of the initial state
$\ket{\psi(0)}=\ket{0}^{\otimes 13}\otimes\ket{101}\otimes\ket{0}^{\otimes
16}$ (the quantum blinker) according to (a) the classical $F_{12}$ rule and
(b,c,d) the Schrödinger equation, Eq. (1). Suplot (a) shows the alive cells as
a function of time according to the classical GOL. Subplots (b) and (c)
display respectively the expectation value $n_{i}$ over the quantum state, Eq.
(3), and the corresponding discretized population $\mathcal{D}_{i}(t)$, Eq.
(4). Subplot (d) shows the surface plot of the clustering function
$\mathcal{C}_{\ell}(t)$, Eq. (6) as a function of time $t$ and cluster size
$\ell$. The time evolution is performed over a lattice of $L=32$ sites using
TDVP with bond dimension $m=32$ and time step $\Delta t=0.01$. Note that in
the classical dynamics (subplot (a)) the evolution occurs over a discrete time
grid of step $\pi/2$.

We consider a spatially localized state with a string $\ket{1,0,1}$ embedded
in a sequence of zeros. Figure 1(a) displays its evolution for the classical
model, in which the distribution of alive cells spreads from the center of the
lattice with an approximate diffusion rate of half lattice cell per time step.
The quantum dynamics, see subplot (b), is instead characterized by dynamical
localization of the cell distributions about the initial positions of the
alive cells. The behaviour becomes more manifest in the dynamics of the
corresponding discretized population, visible in subplot (c). We here observe
collapses and revivals of the alive cells. These features are well captured by
the clustering function, (d), where the dynamics alternate no clusters with
two clusters of alive cells and of unit size. This kind of dynamics has been
first reported in Ref. Bleh et al. (2012) and has been called“quantum blinker”
for the recurrent behavior reminiscent of the blinkers of Conway’s game of
life. It is possible to relate the dynamics of the local population $n_{i}$
with the single-site entropy. We consider a single blinker on a chain with
$L=11$ and evolve it as shown in Fig. 2. In Fig. 2(a) we show the dynamics of
the local occupation $n_{i}$ while in Fig. 2(b) we plot its single-point
entropy. By looking, in particular, at the central site entropy one can notice
an oscillating behavior. We note that the revivals of the initial structure of
the blinker coincide in time with the minima in the single-site entropy,
showing that the entanglement of the chain oscillates and it is minimal when
the revivals occur. This behavior shares analogies with those due to dynamical
symmetries in closed and open quantum systems. Medenjak et al. (2020); Buca et
al. (2019). As visible from comparison with Fig. 1(a), the recurrent features
of this dynamics emerge only in the quantum regime, and are clearly a result
of the interplay between interactions and quantum interference.

Figure 2: (color online) Time evolution of the initial state
$\ket{\psi(0)}=\ket{0}^{\otimes 4}\otimes\ket{101}\otimes\ket{0}^{\otimes 4}$
(the quantum blinker) according to the Schrödinger equation, Eq. (1). Suplot
(a) shows the expectation value $n_{i}$ over the quantum state, Eq. (3), and
the corresponding (b) singe site entropy $S_{i}$. Panel (c) shows the average
concurrence $C(d)$, Eq. 12 The time evolution is performed over a lattice of
$L=11$.

In order to study the spreading of nonlocal correlations we also compute the
central bond entropy by taking a single blinker as initial state and
simulating its dynamics at $L=20$. The corresponding dynamics is displayed in
Fig. 3(a). After a transient, the central bond entropy increases linearly and
then starts performing oscillations about a slowly increasing mean value.
These oscillations are out-of-phase with respect to the oscillations of the
density pattern. Figure 3(b) reports the bond entropy at $t=100$ and as a
function of the lattice sites: the bond entropy is maximum when partitioning
the chain at the center (which is here the center of the blinker) and
decreases linearly when approaching the edges, suggesting that entanglement
follows a volume-law at sufficiently long times. This conjecture is supported
by the scaling of the central bond entropy, computed with tensor networks, as
a function the bond dimension $m$: the entanglement values converge to the
ones of the ODE only for bond dimension $m=256$ (dark brown line), which
corresponds to considering the full Hilbert space. The growth of the
entanglement and the emerging volume law are consistent with the results of
Ref. Schachenmayer et al. (2013) about quench dynamics in short ranged spin
Hamiltonians, where it is shown that the entanglement grows linearly in time
until a saturation value, proportional to the partition size, is reached. This
saturation value is due to the counter propagation of entangled quasi
particles along the lattice. The dynamics we report, moreover, can be compared
with the scenario considered in Ref. Piroli and Cirac (2020), in which a
single quantum cellular automaton (QCA) application is discussed. In fact, we
can cast our problem in terms of the dynamics of a continuous-time QCA over an
infinitesimal time. In these terms our findings are consistent with the area-
law entanglement bound of Ref. Piroli and Cirac (2020).

Network measures provide further information on complexity and insights on the
entanglement properties of the system Doggen et al. (2021). Figure 4(a)
displays the dynamics of the density of entangled spin pairs. The mean value
of this quantity slowly grows. The oscillations have relatively large
amplitude and show that collapses and revivals of the density are accompanied
by oscillations of the quantum correlations. This is a further manifestation
of the quantum dynamics underlying this phenomenon. These oscillations are
also visible in the clustering (b) and in the disparity (c). The relative
amplitude is visibly smaller. The behaviour of clustering reflects the one
visible of the clustering function in Fig. 1(d): the initial increase is
associated with the formation of a large cluster, that then decays into
smaller clusters. The slow asymptotic growth of the clustering can be
attributed to the spreading of correlations across the chain. The disparity,
(c), seems instead to settle about a constant mean value, suggesting that the
backbone of the network has reached a stationary structure about which it
oscillates. The behavior of the network measures can also be related to the
dynamics of the average concurrence $\mathcal{C}(d)$ defined in Eq. (12).
Indeed, the oscillations in the concurrence at $d=1$, shown in Figure 2(c)
resemble those of the density and the clustering highlighted before. The fact
that the concurrence vanishes while the bipartition bond entropy grows
linearly in time, suggests a multipartite entanglement structure.

Figure 3: (color online) Quantum dynamics of correlations across the lattice
with size $L=20$ and for the initial state $\ket{\psi(0)}=\ket{0}^{\otimes
9}\otimes\ket{101}\otimes\ket{0}^{\otimes 8}$ (quantum blinker). (a) Central
bond entropy $S_{\rm Bond}\quantity(t,\frac{L}{2})$, Eq. (10), as a function
of time and (b) bond entropy $S_{\rm Bond}(t,j)$ as a function of the cell $j$
and evaluated at time 100. The curves are determined using TDVP and bond
dimensions $m=32,64,128,256$ (darker colors correspond to higher bond
dimension values, see legenda). The data linked by the dashed line correspond
to a simulation using ordinary differential equations (ODE). The time step of
the evolution is $\Delta t=0.01$. Figure 4: (color online) Characterization of
the dynamics of a quantum blinker using network measures. (a) Density $D(t)$,
Eq. (13), (b) clustering $C(t)$, Eq. (15), (c) disparity $Y(t)$, Eq. (14). The
parameters and legenda are the same as in Fig. 3.

Quantum blinkers are reminiscent of discrete breathers Flach and Willis (1998)
and are a manifestation of quantum interference. We now report the numerical
study of the interaction of an initial state composed by two spatially
separated clusters of cells $\ket{101}$. For infinite distances between the
two clusters, we would observe two quantum blinkers. We plot their dynamics in
Fig. 5 for a finite distance. In particular, Fig. 5(a) shows the dynamics of
the local occupation expectation value $n_{i}$. The surface plot of the
density in Fig. 5(b) shows that blinking structures interfere with each other.
This effective interaction tends to destroy the collapse and revivals of the
discretized population. Nevertheless, the discretized density, Fig. 5(c),
remains localized about the position of the individual blinkers. The
stationary bond entropy, shown in Fig. 6, exhibits a maximum when the bond is
taken at half distance between the two structures, suggesting that the
blinkers become entangled with one another. A similar behaviour has been
reported for other Hamiltonian dynamics simulating quantum cellular automata
Hillberry et al. (2020); Farrelly (2020); Arrighi (2019). Further studies
shall analyse the dependence of these features on the intermediate distance
between the blinkers and on their number across the lattice.

Figure 5: (color online) Dynamics of the interaction between two quantum
blinkers. Subplot (a) is the surface plot of the density of each site as a
function of time for the dynamics of the double-quantum-blinker (“two quantum
blinkers”) state
$\ket{\psi(0)}=\ket{0000}\otimes\ket{101}\otimes\ket{00000}\otimes\ket{101}\otimes\ket{0000}$
and (b) the corresponding discretized population. Subplot (c) gives the
evolution of the clustering function for the two quantum blinkers. The size of
the lattice is $L=19$. The time step is $\Delta t=0.01$. Figure 6: (color
online) Dynamics of the central bond entropy (a) and equilibrium bond entropy
as a function of the lattice sites (b) for one (yellow, light curve, $L=11$)
and two (brown, dark curve, $L=19$) quantum blinkers. The bond entropy for a
single quantum blinker is shown for comparison.

###  III.2 Evolution of a cluster of alive cells.

Figure 7: (color online) Characterization of the time evolution of the initial
state $\ket{\psi(0)}=\ket{0}^{\otimes 4}\otimes\ket{1}^{\otimes
24}\otimes\ket{0}^{\otimes 4}$. The subplots and parameters are the same as
Fig. 1, except that here for the TVDP simulations the bond dimension is
$m=256$.

We now consider a different initial state, composed by a single cluster of
alive cells with size $\sim L$. The dynamics of the local density for the
classical and quantum model are displayed in Fig. 7. The classical dynamics
leads to the propagation of patterns of dead cells from the edges, that seem
to follow a classical trajectory. This leads to a fragmentation of the initial
cluster. After 6 time steps, where the two trajectories starting from the
edges meet at the chain center, the initial cluster of alive cells is splitted
into two large ones. On the contrary, in the quantum dynamics the size of the
cluster decreases as a function of time, such that after a time evolution
corresponding to 6$\pi/2$ the clustering function vanishes. On this time scale
the density reaches a uniform value across the lattice which is below $1/2$,
and thus smaller than the initial average density. The dynamics of the bond
entropy is shown in Fig. 8(a) for $L=16$ and Fig. 8(c) for $L=32$:
Entanglement grows across the lattice until reaching a saturation value. The
spatial profile at long times is shown in Figs. 8(b),(d) and displays a volume
law, as for the quantum blinker. Notice that the entropy spatial profile has
not converged for $L=32$ at bond dimension $m=256$ and $t\geq 4\pi/2$.

Figure 8: (color online) Quantum dynamics of correlations across the lattice
with size $L=16$ (top row) and $L=32$ (bottom row) and for the initial state
$\ket{\psi(0)}=\ket{0}^{\otimes 4}\otimes\ket{1}^{\otimes
12}\otimes\ket{0}^{\otimes 4}$ (top row) and $\ket{\psi(0)}=\ket{0}^{\otimes
4}\otimes\ket{1}^{\otimes 24}\otimes\ket{0}^{\otimes 4}$ (bottom row). (a-c)
Central bond entropy $S_{\rm Bond}\quantity(t,\frac{L}{2})$, Eq. (10), as a
function of time and (b-d) bond entropy $S_{\rm Bond}(t,j)$ as a function of
the cell $j$ . The curves are determined using TDVP and bond dimensions
$m=32,64,128,256$ (darker colors correspond to higher bond dimension values,
see legenda). The data linked by the dashed line correspond to a simulation
using ordinary differential equations. The time step of the evolution is
$\Delta t=0.01$. Figure 9: (color online) Time evolution of the network
measures for the initial state $\ket{\psi(0)}=\ket{0}^{\otimes
4}\otimes\ket{1}^{\otimes 24}\otimes\ket{0}^{\otimes 4}$. The subplots show
(a) the density $D(t)$, Eq. (13), (b) the clustering $C(t)$, Eq. (15), and (c)
the disparity $Y(t)$, Eq. (14). The parameters and legenda are the same as in
Fig. 8 for $L=32$.

The density and clustering measures, Figs. 9(a) and (b), exhibit a fast
increase of correlations, which then quickly decay. After about 6 time steps,
both density and clustering seem to have reached a stationary value which is
small but finite. The behaviour of the disparity, subplot (d), shows that
after an initial transient correlations uniformly spread across the chain.

###  III.3 Evolution of a random Fock state

Figure 10: (color online) Characterization of the time evolution of an initial
random Fock state with density $\rho(0)=0.25$. The subplots and parameters are
the same as Fig. 1, except that here for the TVDP simulations the bond
dimension is $m=64$ and the time step is $\Delta t=0.1$.

The states we discussed so far were characterized by a symmetry in the density
distribution about the center of the lattice. We now study the evolution of an
initial state which is an arbitrary sequence of $0$ and $1$ states along the
lattice, and which we denote by random Fock state. Figure 10 displays an
example of the dynamics of the density distribution for one specific
realisation. The comparison between the classical and the quantum dynamics
exhibits qualitative differences. The quantum dynamics leads to a uniform
average density across the density accompanied by fluctuations in the clusters
dynamics. The average size of the clusters is visibly smaller than in the
classical case.

The long-time behavior is studied by averaging over the dynamics for different
initial states with the same initial density $\rho(0)$ of alive cells. At
sufficiently long times the density of alive cells, Eq. (5), the diversity,
Eq. (8), and the improved diversity, Eq. (9), tend to a constant value. In
Fig. 11 we compare the stationary behaviour for the classical and quantum
models as a function of the initial density of the initial Fock states. We
observe interesting qualitative differences. The classical dynamics predicts
that the final density is a monotonic function of $\rho(0)$ and exceeds
$\rho(0)$ for $\rho(0)\lesssim 0.6$. In the quantum dynamics, on the contrary,
the final density tends to decrease again for $\rho(0)\gtrsim 0.7$. Moreover,
the rate of growth at small initial densities is smaller than the classical
one. The final density is maximal for $\rho(0)$ in the interval about
$[0.6,0.7]$ and exceeds the classical value: quantum dynamics here leads to a
larger number of alive cells than the classical game of life. The diversity
shows that the quantum structures at large initial densities are more complex
than the classical counterparts. The quantum improved diversity, in
particular, increases monotonically with the initial density and behaves
essentially as the quantum diversity, indicating that the clusters of dead
cells are sparse (and of small size). The results obtained taking a larger
lattice (bottom row), are consistent with these observations.

Figure 11: (color online) Characterization of the asymptotic dynamics of a
random Fock state as a function of the initial density $\rho(0)$ and for a
chain with $L=16$ (top) and $L=32$ sites (bottom). Subplots (a,d) display the
equilibrium density $\rho_{\text{equi.}}$ (Eq. (5)), (b,e) display the
equilibrium diversity $\Delta_{\text{equi.}}$ (Eq. (8)), (c,f) the equilibrium
improved diversity $\Delta^{\text{impr.}}_{\text{equi.}}$, Eq. (9)). The black
symbols indicate the results of the classical dynamics, the coloured symbols
the results of TDVP simulations with bond dimensions $m=8,16,32,64$ (darker
colors correspond to higher bond dimension values, see legenda in subplot (e))
and time step $\Delta t=0.1$ (see legenda in subplot (b)). The equilibrium
value is extracted by taking the time average over the time interval where the
numerical results seem to have reached a stationary value. In the quantum case
the interval is $[25,30]$, for the classical evolution it is $[83,100]$. Top
row: The orange symbols are the results of the quantum dynamics integrated
with a standard ODE algorithm with time step $\Delta t=0.01$. In the classical
simulation, we analysed all $65\,536$ different initial configurations. Bottom
row: Each point of the classical simulation is obtained by averaging over
about $5000$ initial states. In the quantum case, each point is determined by
evolving $32$ different initial Fock states (we took the same set of states
when varying the bond dimension $m$).

##  IV Semiclassical Model for classical recurrent configurations

The appearance of quasi-recurrent configurations with the features of the
blinkers of the classical GOL is one remarkable feature of the quantum GOL. In
this section we show an example where quantization of a classical recurrent
structure does not generate quantum recurrent structures. For this purpose we
consider a semiclassical model of the classical rules for the evolution that
combines a generic classically recurrent pattern with the analog unitary
evolution and the superposition of states. In order to identify the quantum
analog of a classically recurrent pattern, let us first recall the classical
$F_{12}$-rule: the state of a given site changes if and only if two or three
sites among the nearest and next-nearest neighbours are alive. In a chain of
length $L$, there are $2^{L}$ possible classical configurations. These
configurations can be classified in terms of their _connectivity_ to other
configurations.

We now loosely employ Fock states in order to denote a sequence of 0 and 1
across the lattice. Let $\ket{\psi_{0}}$ be a Fock state belonging to a
classical recurring configuration with periodicity $n$, and be
$\ket{\psi_{1}},\ldots,\ket{\psi_{n-1}}$ the sequence of states obtained by
applying the classical $F_{12}$ rule, such that
$\ket{\psi_{n}}=\ket{\psi_{0}}$ and $n>1$. Starting from this sequence, we now
construct a corresponding quantum dynamics assuming that it solely couples the
states of this sequence. We note that such dynamics will be implemented by a
specific Hamiltonian, which in the subspace spanned by the states
$\ket{\psi_{1}},\ldots,\ket{\psi_{n-1}}$ shall take the form

$\displaystyle\hat{H}^{(j)}$
$\displaystyle=J\sum_{k=0}^{n-1}|\psi_{k}^{(j)}\rangle\langle\psi_{k+1}^{(j)}|+|\psi_{k+1}^{(j)}\rangle\langle\psi_{k}^{(j)}|\,,$
(16)

where $|\psi_{n}^{(j)}\rangle\coloneqq|\psi_{0}^{(j)}\rangle$. This is
equivalent to the dynamics of a quantum particle in a ring lattice of $n$
sites, where each site is one “classical” configuration of the recurrent
sequence. The eigenstates can be found by using the properties of the
circulant matrices, which we review in the appendix, and take the form

$|\tilde{\psi}_{\ell}\rangle=\sum_{k=0}^{n-1}{\rm
e}^{2\pi\mathrm{i}k\ell/n}\ket{\psi_{k}}/\sqrt{n}$ (17)

with eigenvalue

$E_{\ell}=2J\cos\varphi_{\ell}=2J\cos(2\pi\ell/n)\,.$

Note that these states can be formally derived as the discrete Fourier
Transform of states localized in space. Periodic dynamics can be solely
realised by superpositions of eigenstates whose energy differences are
commensurate, which cannot be realised by means of this construction. This
necessarily implies that the quantum dynamics we constructed from a
classically recurrent structure does not exhibit periodicity in the spatial
density distribution.

##  V Conclusions

The dynamics of a spin lattice, that realises a quantum analog of Conway’s
game of life, exhibits qualitative differences from the corresponding
classical version, even when the initial state is equivalent to a classical
string of bits. These differences are already visible in the dynamics of local
observables. One striking feature is the appearance of metastable structures,
which are reminiscent of classical breathers but have no classical analogon.

We identify several interesting questions that will be object of future
investigations. First, the study of the single-site entropy and the
concurrence suggest a highly nontrivial structure of the entanglement.
Moreover we observe that the dynamics tend to establish entanglement across
the lattice following a volume-law scaling, for which it would be interesting
to investigate whether the predictions of the so-called page curve are
satisfied Fujita et al. (2018); Bianchi and Donà (2019). This property is a
signature of the emerging complexity of the quantum game of life, which calls
for further studies with larger lattice sizes and for its experimental
implementations Wintermantel et al. (2020); Farrelly (2020).

###### Acknowledgements.

We thank Frederic Folz, Rebecca Kraus, Ludger Santen, Tom Schmit and Pietro
Silvi for stimulating discussions and helpful comments. SM acknowledges the
late David Vargas, a brilliant master’s student he had the pleasure to
collaborate with during his stay at Colorado School of Mines, for insightful
discussions that partially inspired this work. This work was funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID
429529648 TRR 306 QuCoLiMa (”Quantum Cooperativity of Light and Matter”) and
by the Priority Program SPP 1929 ”GiRyd” (Giant Interactions in Rydberg
Systems). We acknowledge the support of the German Ministry of Education and
Research (BMBF) via the QuantERA projects NAQUAS and QTFLAG. We acknowledge
support from the EU Horizon 2020 program via the Quantum Flagship PASQUANS
project, and the Italian PRIN 2017. Project NAQUAS and QTFLAG have received
funding from the QuantERA ERA-NET Cofund in Quantum Technologies implemented
within the European Union’s Horizon 2020 program.

##  Appendix A Eigenvectors of a circulant matrix

The matrix corresponding to Hamiltonian (16) is _circulant_ Gray (2009).
Consider a generic circulant matrix of size $n\times n$:

with $c_{j}$ is real, in our case $c_{i}=\delta_{i,1}+\delta_{i,n-1}$. Then
the eigenvalue problem reads

$\displaystyle\sum_{k=0}^{m-1}c_{n+k-m}y_{k}+\sum_{k=m}^{n-1}c_{k-m}y_{k}$
$\displaystyle=Ey_{m}$ (18)

for $m=0,1,\ldots,n-1$, where $E\in\mathbb{C}$ is the eigenvalue, which is
real since $H$ is hermitian, and $y_{m}$ are the components of the
corresponding eigenvector
$\mathbf{y}=\quantity(y_{0},y_{1},\ldots,y_{n-1})^{\intercal}$.

We make the ansatz $y_{k}=\rho^{k}$ and obtain

$\displaystyle\sum_{\alpha=n-m}^{n-1}c_{\alpha}\rho^{\alpha-n}+\sum_{\alpha=0}^{n-1-m}c_{\alpha}\rho^{\alpha}$
$\displaystyle=E.$ (19)

where we divided both sides by $\rho^{m}$ and performed an index shift. We can
choose
$\rho^{-n}=1\quad\Leftrightarrow\quad\rho_{m}=\exp(-2\pi\mathrm{i}\frac{m}{n})$
with $m\in\mathbb{Z}$ such that
$\sum_{\alpha=0}^{n-1}c_{\alpha}\rho^{\alpha}=E$. The corresponding
eigenvector takes the form

$\displaystyle\quantity(\rho^{0},\rho^{1},\ldots,\rho^{n-1})^{\intercal}\,,$
(20)

as one can prove by inserting the eigenvalue-eigenvector pair in Eq. (18):

$\displaystyle\sum_{\alpha=0}^{n-1}c_{\alpha}\rho^{\alpha+m}$
$\displaystyle=\sum_{\alpha=0}^{n-1}c_{\alpha}\rho^{\alpha}\rho_{m}\,.$ (21)

Here, we used $\rho^{m}=\rho^{m+n}=\rho^{m-n}$ for $m=0,1,\ldots,n-1$.
Inserting the values of $c_{j}$ we obtain

$\displaystyle E_{m}$
$\displaystyle=\sum_{\alpha=0}^{n-1}c_{\alpha}\rho^{\alpha}_{m}=2\cos(2\pi\frac{m}{n})\,,$
(22)

with the corresponding eigenvector

$\displaystyle\mathbf{y}_{m}$ $\displaystyle=\matrixquantity(1\\\
\exp(-2\pi\mathrm{i}\frac{m}{n})\\\ \exp(-2\pi\mathrm{i}\frac{2m}{n})\\\
\vdots\\\ \exp(-2\pi\mathrm{i}\frac{(n-1)m}{n}))\eqqcolon\sqrt{n}\ket{y^{m}}.$
(23)

## References

  * Gardner (1970) M. Gardner, Scientific American 223, 120 (1970).
  * Schulman and Seiden (1978) L. S. Schulman and P. E. Seiden, Journal of Statistical Physics 19, 293 (1978), ISSN 1572-9613, URL https://doi.org/10.1007/BF01011727.
  * Rendell (2002) P. Rendell, _Turing Universality of the Game of Life. In: Adamatzky A. (eds) Collision-Based Computing_ (Springer (London), 2002), ISBN 978-1-4471-0129-1, URL https://link.springer.com/chapter/10.1007%2F978-1-4471-0129-1_18.
  * (4) P.Arrighi, Natural Computing 18, 885 (????), URL https://link.springer.com/article/10.1007%2Fs11047-019-09762-6.
  * Brennen and Williams (2003) G. K. Brennen and J. E. Williams, Phys. Rev. A 68, 042311 (2003), URL https://link.aps.org/doi/10.1103/PhysRevA.68.042311.
  * Wintermantel et al. (2020) T. M. Wintermantel, Y. Wang, G. Lochead, S. Shevate, G. K. Brennen, and S. Whitlock, Phys. Rev. Lett. 124, 070503 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.070503.
  * Wolfram (1983) S. Wolfram, Rev. Mod. Phys. 55, 601 (1983), URL https://link.aps.org/doi/10.1103/RevModPhys.55.601.
  * Hillberry et al. (2020) L. E. Hillberry, M. T. Jones, D. L. Vargas, P. Rall, N. Y. Halpern, N. Bao, S. Notarnicola, S. Montangero, and L. D. Carr, _Entangled quantum cellular automata, physical complexity, and Goldilocks rules_ (2020), eprint 2005.01763.
  * Piroli and Cirac (2020) L. Piroli and J. I. Cirac, Phys. Rev. Lett. 125, 190402 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.190402.
  * Bleh et al. (2012) D. Bleh, T. Calarco, and S. Montangero, EPL (Europhysics Letters) 97, 20012 (2012), URL https://doi.org/10.1209/0295-5075/97/20012.
  * Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96–192 (2011), ISSN 0003-4916, URL http://dx.doi.org/10.1016/j.aop.2010.09.012.
  * Montangero (2018) S. Montangero, _Introduction to Tensor Network Methods_ (Springer International Publishing, Cham, 2018), ISBN 978-3-030-01408-7, URL http://link.springer.com/10.1007/978-3-030-01409-4.
  * Grassberger (1986) P. Grassberger, Int. J. Theor. Phys. 25, 907?938 (1986), URL https://link.springer.com/article/10.1007%2FBF00668821.
  * Gros (2015) C. Gros, _Complex and Adaptive Dynamical Systems_ (Springer International Publishing, 2015), ISBN 9783319162652, URL http://dx.doi.org/10.1007/978-3-319-16265-2.
  * Vargas (2016) D. L. Vargas, _Quantum Complexity: Quantum Mutual Information, Complex Networks, and Emergent Phenomena in Quantum Cellular Automata_ (2016), master’s Thesis.
  * Valdez et al. (2017) M. A. Valdez, D. Jaschke, D. L. Vargas, and L. D. Carr, Physical Review Letters 119 (2017), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.119.225301.
  * Cover and Thomas (2006) T. M. Cover and J. A. Thomas, _Elements Of Information Theory_ (Wiley Interscience, 2006).
  * Horodecki et al. (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Reviews of Modern Physics 81, 865–942 (2009), ISSN 1539-0756, URL http://dx.doi.org/10.1103/RevModPhys.81.865.
  * Serrano et al. (2009) M. Á. Serrano, M. Boguñá, and A. Vespignani, Proceedings of the National Academy of Sciences 106, 6483 (2009), ISSN 0027-8424, eprint https://www.pnas.org/content/106/16/6483.full.pdf, URL https://www.pnas.org/content/106/16/6483.
  * Newman (2003) M. E. Newman, SIAM review 45, 167 (2003).
  * Silvi et al. (2019) P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, SciPost Physics Lecture Notes p. 8 (2019), ISSN 2590-1990, eprint 1710.03733.
  * Haegeman et al. (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.165116.
  * OpenMPS Documentation contributors (2020) OpenMPS Documentation contributors, _OpenMPS Documentation_ (2020), online; accessed 30-May-2020, URL https://openmps.sourceforge.io/pyautodoc/mpoterm.html#mpoterm.MBStringTerm.
  * Jaschke et al. (2018) D. Jaschke, M. L. Wall, and L. D. Carr, Computer Physics Communications 225, 59 (2018), ISSN 0010-4655, URL http://www.sciencedirect.com/science/article/pii/S0010465517304204.
  * Wall and Carr (2012) M. L. Wall and L. D. Carr, New Journal of Physics 14, 125015 (2012), URL https://doi.org/10.1088%2F1367-2630%2F14%2F12%2F125015.
  * Medenjak et al. (2020) M. Medenjak, B. Buca, and D. Jaksch, Phys. Rev. B 102, 041117 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.041117.
  * Buca et al. (2019) B. Buca, J. Tindall, and D. Jaksch, Nature Communications 10, 1730 (2019), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-019-09757-y.
  * Schachenmayer et al. (2013) J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. Daley, Physical Review X 3 (2013), ISSN 2160-3308, URL http://dx.doi.org/10.1103/PhysRevX.3.031015.
  * Doggen et al. (2021) E. V. H. Doggen, Y. Gefen, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, _Generalized quantum measurements with matrix product states: Entanglement phase transition and clusterization_ (2021), eprint 2104.10451.
  * Flach and Willis (1998) S. Flach and C. R. Willis, Phys. Reports 295, 181 (1998).
  * Farrelly (2020) T. Farrelly, Quantum 4, 368 (2020), ISSN 2521-327X, URL https://doi.org/10.22331/q-2020-11-30-368.
  * Arrighi (2019) P. Arrighi, _An overview of quantum cellular automata_ (2019), eprint 1904.12956.
  * Fujita et al. (2018) H. Fujita, Y. O. Nakagawa, S. Sugiura, and M. Watanabe, Journal of High Energy Physics 2018 (2018), ISSN 1029-8479, URL http://dx.doi.org/10.1007/JHEP12(2018)112.
  * Bianchi and Donà (2019) E. Bianchi and P. Donà, Physical Review D 100 (2019), ISSN 2470-0029, URL http://dx.doi.org/10.1103/PhysRevD.100.105010.
  * Gray (2009) R. M. Gray, _Toeplitz and Circulant Matrices: A review_ (2009), URL https://ee.stanford.edu/~gray/toeplitz.pdf.