# Fisher zeros and persistent temporal oscillations in non-unitary quantum
circuits

Sankhya Basu  Physics program and Initiative for the Theoretical Sciences, The
Graduate Center, CUNY, New York, NY 10016, USA  Department of Physics and
Astronomy, College of Staten Island, CUNY, Staten Island, NY 10314, USA
Daniel P. Arovas  Department of Physics, University of California at San
Diego, La Jolla, California 92093, USA  Sarang Gopalakrishnan  Department of
Physics, The Pennsylvania State University, University Park, PA 16802, USA
Physics program and Initiative for the Theoretical Sciences, The Graduate
Center, CUNY, New York, NY 10016, USA  Department of Physics and Astronomy,
College of Staten Island, CUNY, Staten Island, NY 10314, USA  Chris A. Hooley
SUPA, School of Physics and Astronomy, University of St Andrews, North Haugh,
St Andrews, Fife KY16 9SS, United Kingdom  Vadim Oganesyan  Physics program
and Initiative for the Theoretical Sciences, The Graduate Center, CUNY, New
York, NY 10016, USA  Department of Physics and Astronomy, College of Staten
Island, CUNY, Staten Island, NY 10314, USA

###### Abstract

We present a quantum circuit with measurements and post-selection that
exhibits a panoply of space- and/or time-ordered phases, from ferromagnetic
order to spin-density waves to time crystals. Unlike the time crystals that
have been found in unitary models, those that occur here are _incommensurate_
with the drive frequency. The period of the incommensurate time-crystal phase
may be tuned by adjusting the circuit parameters. We demonstrate that the
phases of our quantum circuit, including the inherently non-equilibrium
dynamical ones, correspond to complex-temperature equilibrium phases of the
exactly solvable square-lattice anisotropic Ising model.

##  I Introduction

For a many-body quantum system with Hamiltonian operator ${\hat{H}}$, there is
an evident formal similarity between the unitary time-evolution operator,
$e^{-i{\hat{H}}t/\hbar}$, and the density operator for a thermal equilibrium
state, $e^{-\beta{\hat{H}}}$. Since the 1950s this has led to very fruitful
cross-fertilization between the theory of quantum dynamics and the equilibrium
statistical mechanics of quantum systems. Perhaps the most influential of
these is the Matsubara formalismMatsubara (1955), where the thermal density
operator is regarded as an evolution operator in imaginary time: this allows
many of the tools of diagrammatic perturbation theory to be copied more or
less directly from the dynamical to the statistical case.

In recent decades, the development of the theory of open quantum systems has
led to a broadening of interest on the dynamical side of the dynamical-
statistical correspondence, since interactions between the quantum system of
interest and its environment generically induce (effectively) non-unitary
evolution. The quantum circuits in which we shall be mainly interested here
exhibit many-body mixed dynamics, with unitary evolution interrupted by
projection operations meant to model measurements by a classical environment.
Crucially, the many-body system is allowed to continue evolving after such
measurements, and displays a host of novel phenomena due to the tunable
interplay of non-unitary measurements and the intrinsic unitary
dynamicsSkinner  _et al._ (2019); Li  _et al._ (2018, 2019); Chan  _et al._
(2019); Cao  _et al._ (2019); Choi  _et al._ (2020); Gullans and Huse (2019,
2020); Zabalo  _et al._ (2020); Jian  _et al._ (2019); Bao  _et al._ (2020);
Tang and Zhu (2020); Fan  _et al._ (2020); Li and Fisher (2020); Lavasani  _et
al._ (2020); Sang and Hsieh (2020); Ippoliti  _et al._ (2020); Alberton  _et
al._ (2020); Fidkowski  _et al._ (2020).

There have been parallel broadenings of interest on the statistical side.
Starting from the early 1950s, Lee, Yang, Fisher and othersSuzuki and Fisher
(1971); Yang and Lee (1952); Lee and Yang (1952); Griffiths (1969) pioneered
the extension of conventional statistical mechanics to the case where the
coupling constants in the Hamiltonian, or even the inverse temperature $\beta$
itself, are considered to be complex quantities. This opens up the possibility
of points in the complex $\beta$-plane where the partition function vanishes,
something which is not possible for real temperature. For simple models, such
as the isotropic zero-field Ising model on the square lattice, these ‘Fisher
zeros’ occur on contours in the complex $\beta$-plane which cut the real
$\beta$-axis at positions corresponding to the critical temperatures of phase
transitions in the model. If the density of zeros vanishes as the real
$\beta$-axis is approached, the transition is continuous; if the density
remains finite, the transition is first-order.

Figure 1: Phase diagram showing the various spatially and/or temporally
ordered phases of our quantum circuit. PM (paramagnet): short-range order in
both space and time. FM (ferromagnet): uniform long-range order in both space
and time. NFM1 (non-ferromagnet 1): modulated quasi-long-range order in time
(‘algebraic time crystal’). NFM2 (non-ferromagnet 2): modulated quasi-long-
range order in space. The color scheme is used to indicate the evolution of
period of oscillations. In particular, it vanishes exactly on the NFM1-FM
boundary and reaches maximum at the unitary point $\tanh\beta\to i$ where all
damping processes cease and we reach a simple (Floquet) transverse field Ising
model spectrum.

Complex-coupling approaches to statistical mechanics have tended to be seen as
an essentially formal tool. In light of recent progress in understanding the
rich phase structure of _non-unitary_ circuits such as those under continuous
measurement, it is timely to revisit canonical statistical mechanics models at
complex temperature _considered as descriptions of non-unitary evolution_.
This is the perspective we adopt here. Classical measurements generically
introduce randomness in discrete space-time. By contrast, the canonical Ising
model is disorder-free. The correspondence that we shall demonstrate therefore
requires post-selection of measurement outcomes Lee and Chan (2014); Biella
and Schiró (2020); Gopalakrishnan and Gullans (2020); Nahum  _et al._ (2020);
Jian  _et al._ (2021), i.e. only certain outcomes are allowed to continue
evolving and contribute to the eventual disorder-free ensemble of
trajectories. Some of the features of the statistical side of our
correspondence were anticipated in prior work by some of usBeichert  _et al._
, demonstrating the existence of long-range incommensurately modulated
correlations underlying Fisher zeros in the thermodynamics of Ising ladders.

We shall show that the dynamics of the corresponding $M$-qubit circuits
exhibit long-range correlations that coalesce, in the $M\to\infty$ limit, into
extended ordered regions (see Fig. 1). These can be interpreted as phases of
the anisotropic 2D Ising model at complex temperature. They include relatively
conventional short-range ordered ‘paramagnetic’ phases and long-range-ordered
ferromagnetic and antiferromagnetic phases, but also somewhat peculiar
incommensurate critical phases. These latter phases exhibit spatially and/or
temporally modulated correlators with a dynamically determined modulation
period untethered from the underlying lattice. At least one of these latter
phases bears a phenomenological resemblance to the time crystals recently
discussed in the context of unitary dynamics of isolated many-body localized
systemsKhemani  _et al._ (2019). However, it does not fit into the
classification presented in that work, since the circuits we consider are non-
unitary, and the no-go theorems Watanabe and Oshikawa (2015) forbidding time-
crystalline order consequently do not apply.

Before turning to our results, we briefly discuss some connections to
superficially similar questions discussed in previous literature. Temporally
modulated phases in open quantum systems (i.e. limit cycles) have been shown
to exist in more than two spatial dimensions Chan  _et al._ (2015); Lazarides
_et al._ (2020). The non-unitary quantum circuits we consider can be regarded
as Trotterized non-Hermitian Hamiltonians, which have been extensively
explored Ashida  _et al._ (2020); Lee and Chan (2014); Biella and Schiró
(2020); Gopalakrishnan and Gullans (2020). Unlike these works, we keep the
Trotter “time-step” finite, so the models we consider are two-dimensional
statistical mechanics models with a transfer matrix that may be contracted
either sideways or from top to bottom. In addition, many-body entanglement
properties have been computed by contracting transfer matrices sideways in a
series of recent works Bertini  _et al._ (2019); Gopalakrishnan and Lamacraft
(2019); Foss-Feig  _et al._ (2020); Garratt and Chalker (2020a); Lerose  _et
al._ (2020); Garratt and Chalker (2020b); Sonner  _et al._ (2020); Ippoliti
and Khemani (2020); Lu and Grover (2021); Ippoliti  _et al._ (2021), but
primarily in contexts where the dynamics is unitary along one or both
directions. In particular, Refs. Foss-Feig  _et al._ , 2020; Ippoliti and
Khemani, 2020; Brandner  _et al._ , 2017 have proposed experimental protocols
to study the non-Hermitian dynamics of large systems using spacetime duality.
On the statistical mechanics side, Ref. García-Saez and Wei, 2015 used tensor-
network methods similar to those we use here Levin and Nave (2007) to
characterize the thermodynamics of the Yang-Lee model. The present work
applies tools from the complex-temperature statistical mechanics literature
Levin and Nave (2007); García-Saez and Wei (2015) to discuss the unexplored
physics of spatio-temporal correlations in non-unitary quantum circuits. So
far these circuits have primarily been studied for their entanglement
properties; we demonstrate here that even their conventional correlation
functions can exhibit striking phenomena that would be forbidden by unitarity
(in closed systems) or by dimensionality (in open systems Chan  _et al._
(2015); Lazarides  _et al._ (2020)).

The remainder of the paper is organized as follows. In Section II we define
our quantum circuit by explicitly constructing the local gates necessary to
reproduce the complex-temperature statistical mechanics of the Ising model. We
also discuss observables of interest, both ‘thermodynamic’ quantities and two-
point correlation functions, and construct the transfer matrices that govern
the complex-temperature statistical mechanics of interest. In Section III we
study the anisotropic complex-temperature Ising model analytically using
fermionization and also numerically using a tensor-network coarse-graining
scheme. The fermionization treatment is restricted to the case without an
externally applied magnetic field, but the tensor-network coarse-graining
approach allows us to identify additional transitions as a function of field
strength. We conclude with a summary and a discussion of open problems.

##  II Quantum circuits, observables, and transfer matrices

Figure 2: Sketch of the quantum circuit defined in eqn. (1).

###  II.1 Quantum circuit formalism

A quantum circuit consists of a set of qubits, which we may label
$\\{1,\ldots,M\\}$, successively subjected to operations in the form of
quantum gates. In the scenario investigated here, the qubits are Ising spins
and the gates are of two types: two-qubit gates
$W^{\vphantom{\dagger}}_{m,m+1}=\exp(JZ^{\vphantom{\dagger}}_{m}Z^{\vphantom{\dagger}}_{m+1})$
and single qubit gates $V^{\vphantom{\dagger}}_{m}=A\exp(\gamma
X^{\vphantom{\dagger}}_{m})$, where
$\\{X^{\vphantom{\dagger}}_{m},Y^{\vphantom{\dagger}}_{m},Z^{\vphantom{\dagger}}_{m}\\}$
are the Pauli spin matrices for qubit $m$. Defining

$\begin{split}{\cal W}&=W^{\vphantom{\dagger}}_{1,2}\cdots
W^{\vphantom{\dagger}}_{M-1,M}=\exp\\!\bigg{(}J\sum_{m=1}^{M-1}Z^{\vphantom{\dagger}}_{m}\otimes
Z^{\vphantom{\dagger}}_{m+1}\bigg{)}\\\ {\cal
V}&=V^{\vphantom{\dagger}}_{1}V^{\vphantom{\dagger}}_{2}\cdots
V^{\vphantom{\dagger}}_{M}=A^{M}\exp\\!\bigg{(}\gamma\sum_{m=1}^{M}X^{\vphantom{\dagger}}_{m}\bigg{)}\quad,\end{split}$
(1)

the single-step evolution operator ${\cal T}$ for each temporal slice of the
circuit is ${\cal T}={\cal V}{\cal W}$, and the evolution operator for the
full circuit is ${\cal T}^{L}$, where $L$ is the total number of temporal
slices. Note that our qubits are arranged in a finite chain vis-à-vis the
operator ${\cal W}$, rather than a ring with periodic boundary conditions.
This is a common context for models of quantum circuits; what is new here is
that the parameters $J$, $\gamma$, and $A$ are all allowed to be complex,
hence the individual quantum gates are not in general unitary. Acting on a
density matrix $\rho^{\vphantom{*}}_{j}$, the result of each temporal slice of
the circuit is the update $\rho^{\vphantom{*}}_{j+1}={\cal
T}\rho^{\vphantom{*}}_{j}{\cal T}^{\dagger}$. Requiring that
$\textrm{Tr\,}\,\rho^{\vphantom{*}}_{j}$ be preserved then imposes a
relationship among the complex parameters $\\{J,\gamma,A\\}$.

###  II.2 Single-qubit transfer matrix

At each site, our single-qubit transfer matrix $V$ is a product of a unitary
$U=\exp(i\omega\hat{\bm{d}}\cdot{\bm{\sigma}})$ and a POVM (positive operator-
valued measure),

$P(\hat{\bm{n}},\phi\,|\,\alpha)={\textstyle{1\over\sqrt{2}}}\big{(}\cos\phi\
\mathbb{1}+\alpha\sin\phi\,\hat{\bm{n}}\cdot{\bm{\sigma}}\big{)}\quad,$ (2)

where both $\hat{\bm{d}}$ and $\hat{\bm{n}}$ are unit vectors on
$\textsf{S}^{2}$, and where $\alpha=\pm 1$ is the specified measurement
outcome. The fact that
$\sum_{\alpha}P^{\dagger}(\hat{\bm{n}},\phi\,|\,\alpha)\,P(\hat{\bm{n}},\phi\,|\,\alpha)=\mathbb{1}$
is what makes $P(\hat{\bm{n}},\phi\,|\,\alpha)$ a POVM; the fact that there is
one such operator for each measurement outcome $\alpha$ means that the
measurement is ‘efficient’. We choose
$\hat{\bm{d}}=\hat{\bm{n}}=\hat{\bm{x}}$, and we write

$V(\alpha)=e^{i\omega X}\cdot{\textstyle{1\over\sqrt{2}}}\big{(}\cos\phi\
\mathbb{1}+\alpha\sin\phi\,X\big{)}\equiv A\,e^{\alpha\gamma X}$ (3)

where $A=\sqrt{\cos(2\phi)/2}$ and

${\rm Re}\,\gamma=\tanh^{-1}\tan\phi\quad,\quad{\rm Im}\,\gamma=\omega\quad.$
(4)

We define the operator

${\cal V}(\alpha^{\vphantom{*}}_{1},\ldots,\alpha^{\vphantom{*}}_{M})\equiv
V^{\vphantom{*}}_{1}(\alpha^{\vphantom{*}}_{1})\cdots
V^{\vphantom{*}}_{M}(\alpha^{\vphantom{*}}_{M})\quad.$ (5)

Below we shall post-select $\alpha^{\vphantom{*}}_{m}=1$ for all sites, hence

${\cal V}\equiv{\cal
V}(1,1,\ldots,1)=A^{M}\exp\\!\bigg{(}\gamma\sum_{m=1}^{M}X^{\vphantom{*}}_{m}\bigg{)}\quad.$
(6)

The single site transfer matrix $V=A\exp(\gamma X)$ is identical to that of
the one-dimensional classical Ising model
$H=-J_{x}\sum_{j=1}^{L}\sigma^{\vphantom{*}}_{j}\sigma\,^{\vphantom{*}}_{j+1}$
:

$V^{\vphantom{*}}_{\sigma\sigma^{\prime}}=e^{\beta
J_{x}\sigma\sigma^{\prime}}=e^{\beta J_{x}}\big{(}\mathbb{1}+e^{-2\beta
J_{x}}X\big{)}^{\vphantom{*}}_{\sigma\sigma^{\prime}}$ (7)

with $\tanh\gamma=e^{-2\beta J_{x}}$ and $A=\sqrt{\coth\gamma-\tanh\gamma}$ .
Throughout the remainder of this paper we shall set $J_{x}\equiv 1$. Assuming
periodicity in this (temporal) direction, the classical partition function is
${\cal
Z}^{\vphantom{*}}_{L}(\beta)=\textrm{Tr\,}V^{L}=\cosh^{L}\\!\\!\beta+\sinh^{L}\\!\\!\beta$
and the condition ${\cal Z}^{\vphantom{*}}_{L}(\beta)=0$ requires
$\tanh\beta=e^{(2\ell+1)\pi i/L}$, occurring at $L$ equally spaced values
$\ell\in\\{0,\ldots,L-1\\}$ around the circle $|\tanh\beta\,|=1$ in the
complex $\tanh\beta$ plane, as noted by Beichert et al.Beichert  _et al._
Interleaved with these Fisher zeros are $L$ points $\tanh\beta=e^{2\pi
i\ell/L}$ where the correlation function

$\begin{split}C(r;\beta)=\langle\sigma^{\vphantom{*}}_{j}\sigma^{\vphantom{*}}_{j+r}\rangle&=\textrm{Tr\,}(Z\,V^{r}Z\,V^{L-r})\\\
&={\tanh^{r}\\!\beta+\tanh^{L-r}\\!\beta\over
1+\tanh^{L}\\!\\!\beta}\end{split}$ (8)

is long ranged, with $C(r,\beta)=\cos(2\pi\ell r/L)$ (restricting $0\leqslant
r\leqslant L$), corresponding to a wavevector $q={\rm
arg}\tanh\beta=2\pi\ell/L$. In the thermodynamic limit $L\to\infty$, the
Fisher zeros coalesce into a branch cut along the unit circle, with the free
energy exhibiting a simple first-order-like cusp nonanalyticity across the
cut.

Along the contour of Fisher zeros $\tanh\beta=e^{i\theta}$, with $\theta$
real. It follows that $e^{-2\beta}=\tanh\gamma=-i\tan(\theta/2)$. Thus the
entire unit circle in the complex $\tanh\beta$ plane corresponds to simple
unitary stroboscopic precession, i.e. coherent spin-flipping. The value
$\theta=0$ corresponds to a perfectly static spin with no flipping and perfect
persistence, analogous to the ferromagnetic ground state of the corresponding
statistical mechanics problem. The value $\theta=\pi$ corresponds to a
complete spin flip with no identity component and no persistence, analogous to
the antiferromagnetic ground state.

Values of $\tanh\beta$ that do not lie on the unit circle are associated with
circuits that include measurement. Such circuits in the single-qubit case
generically exhibit exponentially decaying temporal correlations, just as the
analogous complex-temperature Ising models exhibit exponentially decaying
spatial correlations. However, as we shall see, the multi-qubit case is
richer: when $M>1$ we shall find ‘decoherence-free subspaces’ along lines in
the complex-temperature plane, though fully unitary evolution occurs only at
isolated points.

###  II.3 Two-qubit transfer matrix

As noted above, our single step evolution operator is a product of single
qubit and two qubit terms. For $M>1$ we introduce the two-qubit transfer
matrix $W^{\vphantom{*}}_{m,m+1}$, which can be expressed as conditioning a
symplectic operation $\exp(\eta JZ^{\vphantom{*}}_{m})$ on the POVM
$P^{\vphantom{*}}_{m+1}(\hat{\bm{z}},{\textstyle{\pi\over 4}}\,|\,\eta)$,
where $\eta=\pm 1$, for each qubit $m\in\\{1,\ldots,L-1\\}$. Explicitly, we
have

$\begin{split}W^{\vphantom{*}}_{m,m+1}&\equiv e^{JZ^{\vphantom{*}}_{m}}\otimes
P^{\vphantom{*}}_{m+1}(\hat{\bm{z}},{\textstyle{\pi\over 4}}\,|\,1)\\\ &\hskip
43.36243pt+e^{-JZ^{\vphantom{*}}_{m}}\otimes
P^{\vphantom{*}}_{m+1}(\hat{\bm{z}},{\textstyle{\pi\over 4}}\,|\,-1)\\\
&=\exp(JZ^{\vphantom{*}}_{m}\otimes Z^{\vphantom{*}}_{m+1})\quad,\end{split}$
(9)

and ${\cal W}\equiv\prod_{m=1}^{M-1}W^{\vphantom{*}}_{m,m+1}$ as given in eqn.
(1).

###  II.4 Mapping to complex-temperature statistical mechanics

The correspondence between statistical mechanics in $d+1$ dimensions and
quantum mechanics in $d$ dimensions is well known Hertz (1976). Traditionally
it requires fine-tuning to a critical point of some sort (the so-called
$\tau$-continuum limit) to enable passing to continuum time where the
correspondence is most powerful. More generally, however, the correspondence
is to a discrete-time ‘kicked’ quantum evolution, of the type exhibited by our
quantum circuit.

By construction, the one-step evolution operator of our quantum circuit
resembles the transfer matrix of a statistical mechanical system:
specifically, an anisotropic 2D Ising model. That model is characterized by
its couplings in the $x$\- and $y$-directions, $J_{x}$ and $J_{y}$
respectively, and by its inverse temperature, $\beta$. These are functions of
the quantum circuit parameters $\\{J,\gamma,A\\}$. In what follows, we shall
consider the subspace of circuit parameters for which the inverse temperature
$\beta$ is complex, while the coupling constants $J_{x}$ and $J_{y}$ are real.

One important question is which direction in our two-dimensional statistical
model will be identified as the time direction in the quantum circuit. In
isotropic models one typically chooses a diagonal direction; here, by
contrast, we choose the $x$-axis of our 2D anisotropic Ising model – the one
with the stronger coupling ($J_{x}>J_{y}$) – to correspond to the time
direction of our quantum circuit.

###  II.5 The $M$-qubit transfer matrix

For complex $\\{J,\gamma,A\\}$, the transfer matrix ${\cal T}={\cal V}{\cal
W}$, which is of dimension $2^{M}\times 2^{M}$, is in general not normal, i.e.
it does not in general commute with its Hermitian conjugate. Nevertheless, any
non-normal complex matrix ${\cal V}$ can be brought to Jordan canonical form
by a similarity transformation ${\cal T}^{\prime}={\cal R}^{-1}{\cal T}{\cal
R}$, where ${\cal R}$ is invertible. If we assume there are no Jordan blocks,
then ${\cal T}$ may be decomposed in terms of its eigenvalues and its left and
right eigenvectors, viz.

${\cal
T}=\sum_{a=0}^{2^{M}-1}\lambda^{\vphantom{*}}_{a}\,{|\hskip-1.2pt|\,R^{\vphantom{*}}_{a}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,L^{\vphantom{*}}_{a}\,|\hskip-1.2pt|}\quad,$
(10)

where
${\langle\hskip-2.2pt\langle\,L^{\vphantom{*}}_{a}\,|\hskip-1.2pt|\,R^{\vphantom{*}}_{b}\,\rangle\hskip-2.2pt\rangle}=\delta^{\vphantom{*}}_{ab}$
, and where there is no complex conjugation implied in the bra vector
${\langle\hskip-2.2pt\langle\,L^{\vphantom{*}}_{a}\,|\hskip-1.2pt|}$ with a
doubled bracket. The eigenvalues $\\{\lambda^{\vphantom{*}}_{a}\\}$ are in
general complex. If we order the eigenvalues such that
$|\lambda^{\vphantom{*}}_{a}|>|\lambda^{\vphantom{*}}_{a+1}|$ for all $a$,
then assuming the largest eigenvalue $\lambda^{\vphantom{*}}_{0}$ is
nondegenerate, after a sufficiently large number of iterations $s$ we have

${\cal
T}^{s}=\lambda_{0}^{s}\,{|\hskip-1.2pt|\,R^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,L^{\vphantom{*}}_{0}\,|\hskip-1.2pt|}+{\cal
O}\big{(}|\lambda^{\vphantom{*}}_{1}/\lambda^{\vphantom{*}}_{0}|^{s}\big{)}\quad.$
(11)

It is convenient to here and henceforth implement a similarity transformation
and redefine ${\cal T}\equiv{\cal V}^{1/2}\,{\cal W}\,{\cal V}^{1/2}$, which
is manifestly symmetric: ${\cal T}={\cal T}^{\top}$. The corresponding right
and left eigenvectors of ${\cal T}$ are then mutual transposes, with no
complex conjugation, which we write as
${|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{a}\,\rangle\hskip-2.2pt\rangle}$ and
${\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{a}\,|\hskip-1.2pt|}$ ,
respectively.

We consider two natural correlation functions which may be used to
characterize the properties of the circuit. The first is the quantum two time
correlator,

$C(s;i,j)={\textrm{Tr\,}\big{[}Z^{\vphantom{*}}_{i}{\cal
T}^{s}Z^{\vphantom{*}}_{j}\,\rho^{\vphantom{*}}_{0}\,({\cal
T}^{\dagger})^{s}\big{]}\over\textrm{Tr\,}\big{[}{\cal
T}^{s}\,\rho^{\vphantom{*}}_{0}\,({\cal T}^{\dagger})^{s}\big{]}}\quad.$ (12)

With $\rho^{\vphantom{*}}_{0}=\mathbb{1}$, we have

$\displaystyle C(s;i,j)$
$\displaystyle={{\langle\hskip-2.2pt\langle\,\Psi^{*}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{i}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{*}_{0}\,\rangle\hskip-2.2pt\rangle}\over\big{|}{\langle\hskip-2.2pt\langle\,\Psi^{*}_{0}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}\big{|}^{2}}\
+$ (13) $\displaystyle\hskip-14.45377pt2\,{\rm
Re}\,\Bigg{\\{}\\!\bigg{(}{\lambda^{\vphantom{*}}_{1}\over\lambda^{\vphantom{*}}_{0}}\bigg{)}^{\\!\\!s}\>{{\langle\hskip-2.2pt\langle\,\Psi^{*}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{i}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{1}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{1}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{*}_{0}\,\rangle\hskip-2.2pt\rangle}\over\big{|}{\langle\hskip-2.2pt\langle\,\Psi^{*}_{0}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}\big{|}^{2}}\Bigg{\\}}+\ldots$

The second is the statistical correlator,

$\displaystyle C^{\vphantom{*}}_{L}(s;i,j)$
$\displaystyle={\textrm{Tr\,}\big{[}Z^{\vphantom{*}}_{i}{\cal
T}^{s}Z^{\vphantom{*}}_{j}{\cal T}^{L-s}\big{]}\over\textrm{Tr\,}\big{[}{\cal
T}^{L}\big{]}}$ (14)
$\displaystyle={\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{i}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}\
+$
$\displaystyle\quad\bigg{(}{\lambda^{\vphantom{*}}_{1}\over\lambda^{\vphantom{*}}_{0}}\bigg{)}^{\\!\\!s}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{i}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{1}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{1}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}\
+$
$\displaystyle\quad\bigg{(}{\lambda^{\vphantom{*}}_{1}\over\lambda^{\vphantom{*}}_{0}}\bigg{)}^{\\!\\!L-s}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{1}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{i}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{1}\,\rangle\hskip-2.2pt\rangle}+\ldots$

Let ${\cal X}=\prod_{j=1}^{M}X^{\vphantom{*}}_{j}$. Since $[{\cal T},{\cal
X}]=0$, assuming
${|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}$ is
nondegenerate, ${\cal
X}\,{|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}=\pm\,{|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}$.
Then

${\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}={\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,{\cal
X}Z^{\vphantom{*}}_{j}{\cal
X}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}=-{\langle\hskip-2.2pt\langle\,\Psi^{\vphantom{*}}_{0}\,|\hskip-1.2pt|\,Z^{\vphantom{*}}_{j}\,|\hskip-1.2pt|\,\Psi^{\vphantom{*}}_{0}\,\rangle\hskip-2.2pt\rangle}$
(15)

and thus the $s$-independent terms in the above two correlators both vanish.
We then have that both $C(s;i,j)$ and $C^{\vphantom{*}}_{L}(s;i,j)$ decay
exponentially in the time direction with a correlation time
$\tau=1/\ln|\lambda^{\vphantom{*}}_{0}/\lambda^{\vphantom{*}}_{1}|$ and a
frequency
$\omega=\textrm{arg}(\lambda^{\vphantom{*}}_{1}/\lambda^{\vphantom{*}}_{0})$
which is generally incommensurate (i.e. irrational). When the spectral gap
collapses, both correlation functions become long-ranged.

At short times we do not expect them to agree, e.g., for unitary circuits,
quantum correlators obey rigid Lieb-Robinson bounds with strictly vanishing
correlators outside the light cone, while statistical correlators are small
but finite for spacelike separations at short times.

It will also be useful to extend some of our expressions from real-temperature
thermodynamics to complex inverse temperature, $\beta$. We examine the modulus
the partition function and define the free energy density accordingly
$f\equiv\log|Z|/ML\beta$, where $M\times L$ is the total number of spins,
followed by the internal energy density and the specific heat capacity

$u\equiv\left|{\partial f\over\partial\beta}\right|\qquad,\qquad
c\equiv\left|{\partial u\over\partial\beta}\right|\quad.$ (16)

We typically plot our results not in the complex $\beta$-plane, but rather in
the complex $\tanh\beta$ plane, which we shall refer to simply as the ‘complex
temperature plane’.

##  III Large-$M$ limit: the anisotropic 2D Ising model

In this Section we consider our quantum circuit in the limit of a large number
of qubits, $M\gg 1$. In the $M\to\infty$ limit, the dynamical correlation
functions of the circuit may be written in terms of the statistical
correlations of an anisotropic 2D Ising model at complex temperature. We
analyze this model using three complementary methods: analytic continuation of
the Onsager solution; numerical evaluation of a tensor-network representation
of the partition and correlation functions; and exact fermionization of the
zero-field problem using the Jordan-Wigner transformation.

###  III.1 Thermodynamics from the Onsager solution

In the isotropic Ising model, the zeros of the partition function lie on
linear contours in the complex-temperature plane. In such a case, provided
that the linear density of zeros reaches a finite value in the thermodynamic
limit, we expect a simple slope discontinuity in the free energy, as already
noted by Fisher Brittin  _et al._ (1964). The anisotropic Ising model was
examined similarly; however, in that case we observe a far more complicated
situation with patterns of zeros that appear to occupy extended regions in the
complex-temperature plane. This makes the expected behavior of the free energy
less clear.

We note at this point that Fisher’s observation and the majority of others
that have followed it are in fact based on a portion of Onsager’s result; as
demonstrated in appendix A, it is manifestly incorrect for the case of finite-
width Ising ladders (corresponding to circuits with a finite number of
qubits). Nevertheless, one might anticipate that the approximation remains
asymptotically exact for 2D bulk (intensive) quantities. The Onsager
expression for the real part of the (dimensionless) free energy per spin of
the anisotropic model is

$\begin{split}\beta f&=\ln 2+{1\over 2}\int\limits_{-\pi}^{\pi}\\!{dk_{x}\over
2\pi}\int\limits_{-\pi}^{\pi}\\!{dk_{y}\over 2\pi}\,\ln\big{|}\cosh j_{x}\cosh
j_{y}\\\ &\hskip 43.36243pt-\,\sinh j_{x}\cos k_{x}-\sinh j_{y}\cos
k_{y}\big{|}\quad,\end{split}$ (17)

where $j_{x,y}\equiv 2\beta J_{x,y}$, where $J_{x}\equiv 1$ and
$J^{\vphantom{*}}_{y}\equiv J$. For the case $J_{y}=0.1$ and complex $\beta$
we have evaluated this for the infinite system numerically. Fig. 1 was
obtained by taking numerous cuts through the complex temperature plane. One
particular cut that is especially revealing is a radial cut away from the
real-temperature axis (Fig. 3) that clearly displays the continuous nature of
the PM-NFM1 transition.

Figure 3: Top left: The internal energy as a function of distance along a
radial line at angle $2\pi/9$ to the positive real axis in the complex
$\tanh\beta$ plane. The two transitions are those into and out of the ‘NFM1’
phase (see Fig. 1). Top right: The specific heat capacity along the same
contour (see text for the precise definition of ‘specific heat capacity’ at
complex temperature). Bottom panels: The same specific heat graph, but
multiplied by the specified factor, demonstrating that the two singularities
are one-sided square-root singularities.

###  III.2 Correlations from tensor-network renormalization

We would like to characterize the different phases that appear in Fig. 1,
especially the NFM1 and NFM2 phases that are not simple continuations of real-
temperature phases. For this we need to know the spin-spin correlation
functions in both the $x$ (time) and $y$ (qubit array) directions. Our most
general method for determining these, which has the additional advantage of
allowing the inclusion of a longitudinal magnetic field, is via a
renormalization group algorithm based on tensor networks - in particular, the
Tensor Renormalization Group (TRG) Levin and Nave (2007), a method that
involves representing classical partition functions as tensor networks and
coarse-graining these tensor networks numerically.

Most of our results are obtained with bond dimensions up to $50$, and we use
relatively modest convergence goals which we check throughout, e.g. that the
free energy density is converged to $\sim 0.001$. In the $J_{y}\to 0$ limit,
our system is a set of uncoupled Ising chains; we know that, in this limit,
the entire complex-temperature plane is paramagnetic with the exception of the
unit circle $|\tanh\beta|=1$. We therefore expect that, for $J_{y}\ll 1$,
correlated non-paramagnetic phases will be concentrated near the unit circle.

We have used finite-field TRG to establish the disappearance of uniform
ferromagnetic order as we traverse the unit circle $|\tanh\beta|=1$
counterclockwise from the real-temperature line, for several values of the
coupling $J_{y}$ — see Fig. 4. It is clear that the FM phase is progressively
suppressed as the interchain coupling is reduced. We interpret this as the
gradual reversion to the incommensurately modulated order seen on the unit
circle in the case of decoupled chains. This perspective already suggests that
the NFM1 phase exhibits some form of long-range incommensurate order; we show
below that that is essentially true.

Figure 4: The spontaneous magnetization of the anisotropic 2D Ising model as a
function of angle from the real-temperature axis on the contour
$|{\tanh\beta}|=1$, for four different values of the interchain coupling
$J_{y}$. The jump corresponds to the FM-NFM1 transition, shown for the
$J_{y}=0.1$ case in Fig. 1.

Additional work is required to use TRG to compute correlation functions. As
the method is multiplicative, i.e. distances are reduced by a factor of 2 per
iteration, evaluating the correlation function at separations $2^{n}$ is
relatively easy. These are useful in cases for which we expect simple power-
law decays. Here, however, we are interested in modulated correlators, which
requires careful renormalization at short distances.

Our TRG-computed correlation functions in the NFM1 phase are shown in the
left-hand panels of Fig. 5. We observe modulated correlations along the
direction of strong coupling ($x$ in the statistical mechanics setting; time
in the quantum circuit picture). The correlations along the weak direction
($y$ in the statmech picture; inter-qubit in the quantum circuit picture) are
non-oscillatory and apparently power-law decaying. In the next section we
shall provide an interpretation of these results in terms of a fermionized
version of the model.

Figure 5: Spin-spin correlations in the NFM1 and NFM2 phases, obtained via
tensor renormalization group (TRG) algorithm. Left-hand panels: The spin-spin
correlation function in the $x$ (time) direction (upper panel) and in the $y$
(intra-register) direction (lower panel) in the NFM1 phase, at the point
$\tanh\beta=e^{i\pi/6}$. Middle panels: The same, but in the NFM2 phase, at
the point $\tanh\beta=1.3e^{i\pi/90}$. Notice that the oscillations in the
NFM2 phase run in the $y$ (intra-register) direction, while those in the NFM1
phase run in the $x$ (time) direction. Right-hand upper panel: The spin-spin
correlation function in the $x$ (time) direction in the NFM1 phase, at the
point $\tanh\beta=e^{i\pi/3}$, for three different values of the applied
longitudinal magnetic field. Right-hand lower panel: The magnetization
$m=-\partial\,{\rm Re}(\beta f)/\partial h$ as a function of applied external
magnetic field $h$, at various points $\tanh\beta=e^{i\theta}$. The point
$\theta=10^{\circ}$ is in the FM phase; the others are in the NFM1 phase.

Our TRG results for the correlation functions in the NFM2 phase are shown in
the middle panels of Fig. 5. Surprisingly, we discover that the directions of
the modulated and simply-decaying correlations are swapped in NFM2 relative to
NFM1. It is worth noting that for the NFM1 phase the $J_{y}\to 0$ limit is
solvable and contains modulated correlations already, whereas the NFM2 phase
does not exist in the decoupled limit, owing its existence to interchain
interactions.

The upper right-hand panel of Fig. 5 shows the correlation function in the $x$
(time) direction at a different point in the NFM1 phase, where the period of
the temporal oscillations is shorter. It also shows the effect on these
correlations of the application of a longitudinal magnetic field. We see that
they survive almost entirely unaltered, i.e. that this phenomenon is
insensitive to the breaking of integrability.

The lower right-hand panel of Fig. 5 shows the magnetization as a function of
longitudinal field in the modulated phases. Up to some critical field (which
depends on the extent of anisotropy), these modulated phases are stable. At
this critical field a metamagnetic transition takes place and the system exits
the modulated phase.

The TRG method makes it straightforward to compute correlations even for large
systems, for particular exponentially spaced sets of distances. Evaluating the
spatial decay of the correlator in the NFM1 phase, we see clear signs of
algebraic decay (Fig. 6), exactly as one would expect in a one-dimensional
critical phase.

###  III.3 Fermionization and the origin of oscillations

Finally, we present an approach to the zero-field problem that uses the
Jordan-Wigner transformation to represent the spins/qubits in terms of
fermionic degrees of freedom. We shall show that the occurrence of
correlations is due to a certain type of resonance between two eigenvalues of
these fermionic operators. We may use this picture to predict both the
temporal period of the oscillations in the NFM1 phase and the spatial period
of the oscillations in the NFM2 phase; in both cases, we find good agreement
with our TRG results presented above.

In the Jordan-Wigner representation, we write the Pauli matrices on site $j$
as follows:

$\begin{split}X^{\vphantom{*}}_{j}&=2c^{\dagger}_{j}c^{\vphantom{\dagger}}_{j}-1\\\
Z^{\vphantom{*}}_{j}&=(c^{\dagger}_{j}+c^{\vphantom{\dagger}}_{j})\,\prod_{l=1}^{j-1}e^{i\pi
c^{\dagger}_{l}c^{\vphantom{\dagger}}_{l}}\quad,\end{split}$ (18)

where the operator $c_{j}$ annihilates a (spinless) fermion on site $j$ of the
qubit register. It follows that

$Z^{\vphantom{\dagger}}_{j}Z^{\vphantom{*}}_{j+1}=c^{\dagger}_{j}c^{\vphantom{\dagger}}_{j+1}+c^{\dagger}_{j+1}c^{\vphantom{\dagger}}_{j}-c^{\vphantom{\dagger}}_{j}c^{\vphantom{\dagger}}_{j+1}-c^{\dagger}_{j+1}c^{\dagger}_{j}\quad.$
(19)

We now Fourier transform the fermion operators with respect to the index $j$,
i.e. we move to a plane-wave basis in our qubit register:

$c_{j}={\textstyle{1\over\sqrt{M}}}{\sum_{k}}^{\prime}e^{-ikj}c_{k}\quad,\quad
c^{\dagger}_{j}={\textstyle{1\over\sqrt{M}}}{\sum_{k}}^{\prime}e^{+ikj}c^{\dagger}_{k}\quad,$
(20)

where the prime on the sum restricts $k$ to the first Brillouin zone, i.e.
$k\in[-\pi,\pi)$. In terms of these plane-wave operators, the components of
our single-step time-evolution operator become

$\displaystyle{\cal W}$
$\displaystyle=\exp\\!\bigg{(}2J\,{\sum_{k>0}}^{\prime}\Big{[}\cos
k\,\big{(}c^{\dagger}_{k}c^{\vphantom{\dagger}}_{k}+c^{\dagger}_{-k}c^{\vphantom{\dagger}}_{-k}\big{)}$
(21) $\displaystyle\hskip 108.405pt+i\sin
k\,\big{(}c^{\vphantom{\dagger}}_{-k}c^{\vphantom{\dagger}}_{k}-c^{\dagger}_{k}c^{\dagger}_{-k}\big{)}\Big{]}\bigg{)}$
$\displaystyle{\cal V}$
$\displaystyle=\exp\\!\bigg{(}2\gamma\,{\sum_{k>0}}^{\prime}\Big{[}c^{\dagger}_{k}c^{\vphantom{\dagger}}_{k}+c^{\dagger}_{-k}c^{\vphantom{\dagger}}_{-k}-1+\gamma^{-1}\ln
A\Big{]}\bigg{)}\quad,$

with ${\cal T}={\cal V}^{1/2}\,{\cal W}\,{\cal V}^{1/2}$.

We may streamline our notation using Anderson pseudospin operators
$\tau^{\alpha}_{k}$, defined as follows:

$\tau_{k}^{\alpha}\equiv\begin{pmatrix}c^{\dagger}_{k}&c^{\vphantom{\dagger}}_{-k}\end{pmatrix}\sigma^{\alpha}\begin{pmatrix}c^{\vphantom{\dagger}}_{k}\\\
c^{\dagger}_{-k}\end{pmatrix}\quad,$ (22)

where $\alpha\in\\{0,1,2,3\\}$. In terms of these operators, the gates can now
be re-written as follows:

$\displaystyle{\cal W}$
$\displaystyle={\prod_{k>0}}^{\prime}\exp\\!\Big{(}2J\,\big{(}\tau^{z}_{k}\cos
k+\tau^{y}_{k}\sin k\big{)}\Big{)}$ (23) $\displaystyle{\cal V}$
$\displaystyle={\prod_{k>0}}^{\prime}A^{2}\exp\\!\big{(}2\gamma\tau^{z}_{k}\big{)}$
(24) $\displaystyle{\cal V}^{1/2}$
$\displaystyle={\prod_{k>0}}^{\prime}A\exp\\!\big{(}\gamma\tau^{z}_{k}\big{)}\quad.$
(25)

The partition function ${\cal Z}=\textrm{Tr\,}\big{(}{\cal T}^{L}\big{)}$ may
be expressed as the product

$\displaystyle{\cal Z}$
$\displaystyle={\prod_{k>0}}^{\prime}\textrm{Tr\,}\big{[}\Theta_{k}^{L}\big{]}$
(26) $\displaystyle\Theta^{\vphantom{*}}_{k}$ $\displaystyle\equiv
A^{2}\exp(\gamma\tau^{z}_{k})\exp\\!\big{[}2J(\tau^{z}_{k}\cos
k+\tau^{y}_{k}\sin k)\big{]}\exp(\gamma\tau^{z}_{k})\quad.$

For each wavevector $k>0$, the operator $\Theta^{\vphantom{*}}_{k}$ has two
eigenvalues, $\lambda_{k,\pm}$ (see appendix C),

$\lambda^{\vphantom{*}}_{\pm}(k)=2\big{(}h^{\vphantom{*}}_{k}\pm\delta^{\vphantom{*}}_{k}\big{)}\quad,$
(27)

where $h^{\vphantom{*}}_{k}$ and $\delta^{\vphantom{*}}_{k}$ are given by

$h^{\vphantom{*}}_{k}=\cosh(2\beta)\cosh(2J)+\sinh(2J)\cos(k)$ (28)

and

$\displaystyle\delta^{\vphantom{*}}_{k}$
$\displaystyle=\Big{[}\sinh^{2}(2\beta)\sinh^{2}(2J)\sin^{2}\\!k$ (29)
$\displaystyle\hskip 14.45377pt+\big{[}\cosh(2J)+\cosh(2\beta)\sinh(2J)\cos
k\big{]}^{2}\,\Big{]}^{1/2}\quad,$

where we have used a mixed notation, trading coupling constants $\gamma$ and
$A$ for $\beta$ – recall the relations $\tanh\gamma=e^{-2\beta}$ and
$A=\sqrt{\coth\gamma-\tanh\gamma}$. To analyze the late-time properties of the
evolution, we find the largest-amplitude eigenvalue of
$\Theta^{\vphantom{*}}_{k}$ (eq.26) for each wavenumber $k$. We denote this as
$\lambda^{\vphantom{*}}_{0}(k)$ and the corresponding right eigenvector as
${|\hskip-1.2pt|\,\psi^{\vphantom{*}}_{0}(k)\,\rangle\hskip-2.2pt\rangle}$.
These eigenvalues determine the decay rate and precession of a typical initial
condition at late times ($t\to\infty$):

${|\,0\,\rangle}\to{\prod_{k>0}}^{\prime}\big{[}\lambda^{\vphantom{*}}_{0}(k)\big{]}^{t}\,{|\hskip-1.2pt|\,\psi^{\vphantom{*}}_{0}(k)\,\rangle\hskip-2.2pt\rangle}\quad.$
(30)

Figure 6: Correlations along the non-oscillatory $y$ direction on a log-log
scale, extracted from TRG for $\tanh\beta=e^{i\pi/6}$ and $J_{y}=0.1J_{x}$
(NFM1 phase), showing clear evidence of power-law decay with exponent $\sim
0.33$.

We now discuss the behavior of this late-time state in terms of properties of
the fermionic spectrum. Outside the NFM1 phase, either the $+$ or $-$ branch
is consistently larger-amplitude throughout the Brillouin zone, and there is
generically a unique steady state
${|\hskip-1.2pt|\,0\,\rangle\hskip-2.2pt\rangle}$. In the NFM1 phase, however,
the branches “invert” as a function of $k$, which is to say
$|\lambda^{\vphantom{*}}_{+}(k)|>|\lambda^{\vphantom{*}}_{-}(k)|$ for
$k\approx 0$ but the opposite inequality holds for $k\approx\pm\pi$. At
special momenta $\pm k^{*}$, the two eigenvalues are degenerate. Therefore, in
the subspace of $\pm k^{*}$, the system never reaches a unique steady state,
and instead one has persistent oscillations at the frequency

$\omega^{*}=\big{|}\arg\lambda^{\vphantom{*}}_{+}(k^{*})-\arg\lambda^{\vphantom{*}}_{-}(k^{*})\,\big{|}\quad.$
(31)

The band inversion point $k^{*}$ sweeps across the Brillouin zone as one
progresses through the NFM1 phase, leading to incommensurate temporal
modulations of varying frequency.

By comparing with the numerical solutions we see that this phase exhibits
_temporal_ oscillations but no apparent spatial oscillations. While the
momentum $2k^{*}$ appears to be special in some sense, from the above
argument, there is no simple relation between spectral degeneracies of the
sort described above and spatial oscillations. To capture modulated
correlations in the NFM2 phase, it is convenient instead to fermionize the
model _sideways_ , by performing the Jordan-Wigner transformation along the
$x$ axis (which hitherto we took to be the temporal direction). The “band
inversion” described above now happens in the NFM2 phase, leading to
oscillations in the _spatial_ direction (i.e., along $y$).

Table 1: Comparison of numerically estimated periods of order parameter
oscillations ($T_{\rm TRG}$) against exact Jordan-Wigner fermion periods,
$T_{\rm JW}\equiv 4\pi/\omega^{*}$. Note the addtional factor of 2 due to a
two-site unit-cell in the time-direction implicit in the definition of the one
step evolution operator. First three rows corresponds to points inside NFM1,
where modulations are along the $x$-axis, while the last row is in NFM2.

##  IV Discussion

In the present work we have exploited the correspondence between non-unitary
quantum circuits and complex temperature statistical mechanics to construct a
simple quantum circuit that has a surprisingly rich phase diagram, including a
phase with incommensurate temporal order. Such incommensurate time crystals do
not seem to occur in closed systems; nor do they occur in one-dimensional open
quantum systems with short-range interactions, for entropic reasons. Our
results suggest that an important class of quantum circuits that exhibit
incommensurate time-crystalline order are spacetime duals of circuits that
realize incommensurate density-wave phases. In this simple one-dimensional
case, such phases occur only for complex couplings, but in more general
settings it might be possible to write down quantum circuits that cool the
system into a ground state with density-wave order Stoudenmire and White
(2010). These would also be spacetime dual to temporally modulated phases.

In practice, post-selection is an expensive operation requiring effort that
scales exponentially in the area of the quantum circuit. Thus, practical
realizations of the physics discussed here will be restricted to circuits that
are either very shallow or involve only a small number of qubits evolved for a
long time. These map onto Ising ladders at complex temperature, which can be
solved using the methods discussed above (App. A). We find that signatures of
the modulated phases are present even for systems with modest numbers of spins
($M=5$), which should be realistic to explore in a variety of present-day
experiments.

##  V Acknowledgments

We are grateful to F. Beichert, R. Moessner, F. Pollmann, D. Huse, E. Fradkin,
S. Kivelson, W. Bialek, Alexander Abanov, M. Stoudenmire, Tzu-Chieh Wei, N.
Pomata, John Cardy and Michael Fisher for stimulating discussions over the
years. CAH gratefully acknowledges financial support from UKRI under grant
number EP/R031924/1; he is also grateful to Rice University for a visiting
appointment in spring 2019, where part of this work was completed. SB and VO
acknowledge support from the NSF DMR Grant No. 1508538 and US-Israel BSF Grant
No. 2014265.

##  Appendix A Complex-temperature statistical mechanics of Ising ladders

In this appendix, we review previously obtained results on loci of Fisher
zeros as they coalesce into branch cuts for $M>1$. We also present some new
results where we identify the regions of the complex-temperature plane in
which the correlation lengths exceed ten lattice spacings. These latter
results foreshadow the form of the phase diagram in the $M\to\infty$ limit,
Fig. 1.

Quite generally, the only place where long-range order can occur in an
infinite Ising system is on the same contours where the partition function
zeros for the finite system concentrate. This is because the condition for the
two is the same, viz. that the largest two eigenvalues of the transfer matrix
become equimodular. As explicitly demonstrated analytically for $M=1$ and
numerically for $M>1$, Fisher zeros merge into $\approx M$ branch cuts in the
complex temperature plane in the thermodynamic (long-time) limit, thus
allowing for a smooth evolution of correlation functions, with the correlation
time becoming infinite on branch cuts and the correlations themselves
retaining an oscillatory character inherited from the relative phase of the
two dominant eigenvaluesBeichert  _et al._ .

What happens to the spin-spin correlation length (decay time) in the regions
between the contours? The answer is that, between the $M$ contours that are
‘coalescing’ into the ordered region in the 2D Ising model, the correlation
length stays very high. The existence of ferromagnetism over a finite region
of the complex-temperature plane may be anticipated by noticing that phases of
the largest and second-largest eigenvalues ‘lock’ to each other (as
$M\to\infty$). This appears to be the only type of correlated ‘phase’ that
occurs in the isotropic case. Anisotropic lattices, however, appear to support
another type of ‘gapless’ correlated phase, which exhibits multiple long
length-scales, and which shows precursor signatures in the behavior of the
correlation lengths in the finite-$M$ case.

We begin by reviewing our prior resultsBeichert  _et al._ , in which the
limiting behavior of Fisher zeros in ladders was computed — see Fig. 7. It can
be shown, via a low-temperature expansion, that the number of contours (branch
cuts) emerging from the two zero-temperature points $\tanh\beta=\pm 1$ is
equal to $M$. It is less clear how to compute the total number of contours,
although contours that do not emerge from $\tanh\beta=\pm 1$ do not appear to
proliferate, and may be strongly dependent on boundary conditions.

Fisher’s original proposal overlooked this behavior entirely. The exact
solutions of the Ising model by Onsager and several others in the years that
followedOnsager (1944) usually consist of several contributions, only one of
which dominates in the thermodynamic limit. Fisher’s original argument for
generalizing Yang-Lee results was based on a seemingly incorrect procedure
whereby he analytically continued only the portion of the result that was
important at real temperature. As explicitly demonstrated by Beichert et al.,
this produces entirely wrong patterns of zeros in ladders. Remarkably,
however, Fisher’s approximate solution is accurately reproduced by the
unbiased TRG computational scheme applied to the 2D lattice.

Figure 7: (Reproduced from F. Beichert et al.Beichert  _et al._ ) The contours
containing the zeros of the partition function in the complex temperature
plane for Ising spin ladders with 4 legs and isotropic spin-spin interactions
($J_{x}=J_{y}$). Dashed lines represent the expected location of Fisher zeros
of the isotropic 2D Ising modelBrittin  _et al._ (1964). The total number of
contours is equal to the number of legs; their location was found to be only
weakly sensitive to boundary conditions in the short (rung) direction and
insensitive to the boundary conditions in the long (infinite) direction.

Next we examine the growth of correlation lengths in the $x$ (time) direction
as we increase $M$. Each correlation length is controlled by the ratio between
one of the subdominant eigenvalues of the transfer matrix, $\lambda_{j}$
($j>0$), and the dominant one, $\lambda_{0}$. In the left-hand panels of Fig.
8, the shaded areas mark the regions of the complex-temperature plane in which
the longest correlation length, i.e. the one controlled by
$\lambda_{1}/\lambda_{0}$, is greater than ten lattice spacings. This is shown
for the two-qubit case (top left) and for the five-qubit case (bottom left).
It is clear that our arbitrarily determined threshold of 10 lattice sites is
already exceeded for $M=5$ in the entire crescent region outlined by Fisher’s
original 2D proposal. The right-hand panels show the same information but for
the second-longest correlation length, i.e. the one controlled by
$\lambda_{2}/\lambda_{0}$. We note that there is no other long correlation
length present; as we can see from these plots, the eigenvalue $\lambda_{2}$
does not approach $\lambda_{0}$ except near the unitary point $\tanh\beta=i$,
i.e. the ferromagnetic phase is ‘gapped’.

$M=2$, longest $M=2$, second-longest

$M=5$, longest $M=5$, second-longest

Figure 8: The top-right quadrant of the complex $\tanh\beta$ plane of Fig. 7
for an $M$-leg Ising spin ladder with isotropic interactions, $J_{x}=J_{y}=1$.
Upper panels: for a 2-leg ladder, the longest correlation length (left) and
the second-longest (right), both in the ‘long’ direction, i.e. along the legs
of the ladder. Lower panels: the same for a 5-leg ladder. The regions in which
the correlation length in question exceeds 10 lattice sites are shaded black.
Note that, except near the unitary point $\tanh\beta=i$, only the longest
correlation length shows significant structure.

We next turn to the case of anisotropic couplings, $J_{y}<1$. It was observed
by van Saarloos and Kurtze van Saarloos and Kurtze (1984) that, in this case,
Fisher’s approximation to the partition function produces highly complex
patterns of zeros. For simple integer fractions $J_{y}=1/n$ it is possible to
compute and plot the contoursS. Basu et al. and observe an erratic pattern
that does not show simple convergence to the limit of decoupled qubits.
Numerically exact computations for finite-$M$ transfer matrices, however,
produce a nicely regular and convergent pattern, from which Fig. 9 was
obtained. The ordered phase is reduced in its extent as we have reduced
$J_{y}$ tenfold. As with the isotropic model we observe a growth of the
longest correlation length as we increase $M$. However, in this anisotropic
case we also see additional regions away from the unitary point in which the
second-longest correlation length also becomes long: one tracking the unit
circle, and another small patch on the real $\tanh\beta$ axes at values
exceeding 1. These locations are suggestively similar in shape and location to
the NFM1 and NFM2 phases in Fig. 1.

$M=2$, longest $M=2$, second-longest

$M=5$, longest $M=5$, second-longest

Figure 9: The same as Fig. 8, but for the case of anisotropic interactions,
$J_{y}=J_{x}/10=0.1$. Note that, in a region roughly corresponding to the
‘NFM1’ and ‘NFM2’ regions of Fig. 1, the second-longest correlation length
also shows non-trivial structure.

##  Appendix B TRG results for the magnetization : isotropic case

Figure 10: Magnetization, $m=-\partial\,{\rm Re}(\beta f)/\partial h$ as a
function of applied uniform magnetic field in the two phases of the isotropic
2D Ising model in the complex $\sinh(2\beta)$ plane. Upper panel:
ferromagnetic phase, $|\sinh(2\beta)|=1.3$. Lower panel: metamagnetic response
in the paramagnetic phase, $|\sinh(2\beta)|=0.9$.

In this appendix we present our TRG results for the behavior of the isotropic
Ising model in a longitudinal field.

Based on the patterns of Fisher zeros, it seems that the isotropic Ising model
in zero applied magnetic field exhibits, for general complex temperature, a
first-order transition between a paramagnetic (PM) and a ferromagnetic (FM)
phase. We explicitly verify this through a direct computation of the
magnetization process in the vicinity of that phase boundary.

The phase boundary between the PM and FM phases forms the famous double
crescent in the $\tanh\beta$ plane, but becomes a yet simpler unit circle in
the complex $\sinh 2\beta$ plane: this appears to be a more natural variable
for the isotropic case. A trajectory with a fixed modulus of $\sinh 2\beta$
slightly above/below 1 allows us to study the evolution of the character of
the transition as we move away from the real-temperature critical point.

In Fig. 10, we show our results for the magnetization as a function of field
for a variety of angles along two such circular contours in the complex $\sinh
2\beta$ plane. In the upper panel, we observe that the FM phase indeed
exhibits spontaneous magnetization, the magnitude of which drops to zero as
the point $\arg\sinh 2\beta=\pi/2$ is approached. The PM phase shows initially
linear response, followed by a metamagnetic jump, as expected near
conventional first-order transitionsChaikin  _et al._ (1995). Interestingly,
there is a large domain of _diamagnetic_ response, which would be a rather
peculiar state of affairs for a more conventional statistical problem.

##  Appendix C Fermionization

In this appendix, we provide the details of the various steps of our
fermionization approach to the zero-field case.

In §III.3 it was shown that the partition function ${\cal
Z}=\textrm{Tr\,}\big{(}{\cal V}^{1/2}\,{\cal W}\,{\cal V}^{1/2}\big{)}^{\\!L}$
is given by ${\cal
Z}=\prod_{k>0}^{\prime}\textrm{Tr\,}\big{[}\Theta_{k}^{L}\big{]}$, where

$\Theta^{\vphantom{*}}_{k}\equiv
A^{2}\,e^{\gamma\tau^{z}_{k}}\,e^{2J(\tau^{z}_{k}\cos k+\tau^{y}_{k}\sin
k)}\,e^{\gamma\tau^{z}_{k}}\quad.$ (32)

For our purposes here, considering one $k$ mode at a time, we write
${\bm{\tau}}=\\{X,Y,Z\\}$ for the Pauli matrices
$\big{\\{}\tau^{x}_{k},\tau^{y}_{k},\tau^{z}_{k}\big{\\}}$. We then have

$e^{\gamma Z}=\cosh\gamma+Z\sinh\gamma$ (33)

and

$\displaystyle e^{2J(Y\\!\sin k+Z\cos k)}$ $\displaystyle=\cosh(2J)$ (34)
$\displaystyle\hskip 28.90755pt+\sinh(2J)\,(Y\\!\sin k+Z\cos k)\quad.$

We find
$\Theta^{\vphantom{*}}_{k}=2\big{[}d^{\vphantom{*}}_{0}(k)+{\bm{d}}(k)\cdot{\bm{\tau}}\big{]}$
with

$\begin{split}d^{\vphantom{*}}_{0}(k)&=\cosh(2\beta)\cosh(2J)+\sinh(2J)\cos
k\\\ d^{\vphantom{*}}_{x}(k)&=0\\\
d^{\vphantom{*}}_{y}(k)&=\sinh(2\beta)\sinh(2J)\sin k\\\
d^{\vphantom{*}}_{z}(k)&=\cosh(2J)+\cosh(2\beta)\sinh(2J)\cos
k\quad.\end{split}$ (35)

Thus the eigenvalues of $\Theta^{\vphantom{*}}_{k}$ are given by

$\begin{split}\lambda^{\vphantom{*}}_{\pm}(k)&=2d^{\vphantom{*}}_{0}(k)\pm
2\sqrt{d_{y}^{2}(k)+d_{z}^{2}(k)}\\\
&=2d^{\vphantom{*}}_{0}(k)\,\Big{\\{}1\pm\sqrt{\Omega_{+}(k)\,\Omega_{-}(k)}\Big{\\}}\end{split}$
(36)

where

$\Omega^{\vphantom{*}}_{\pm}(k)\equiv{d^{\vphantom{*}}_{z}(k)\pm
id^{\vphantom{*}}_{y}(k)\over d^{\vphantom{*}}_{0}(k)}\quad.$ (37)

Resonances : The resonance condition
$\big{|}\lambda^{\vphantom{*}}_{+}(k)\big{|}=\big{|}\lambda^{\vphantom{*}}_{-}(k)\big{|}$
thus pertains when

$\Omega_{+}(k)\,\Omega_{-}(k)=-\alpha^{2}\in\mathbb{R}_{-}\quad.$ (38)

This entails

$\begin{split}0&=d_{+}(k)\,d_{-}(k)+\alpha^{2}d_{0}^{2}(k)\\\
&=(1+\alpha^{2})\cosh^{2}(2\beta)\cosh^{2}(2J)-\sinh^{2}(2\beta)\\\ &\qquad\
+(1+\alpha^{2})\cosh(2\beta)\sinh(4J)\cos k\\\
&\qquad\qquad+(1+\alpha^{2})\sinh^{2}(2J)\cos^{2}\\!k\quad,\end{split}$ (39)

which is a quadratic equation in $\cos k$, with the solution

$\cos
k={1\over\sinh(2J)}\Bigg{\\{}\\!\\!-\cosh(2\beta)\cosh(2J)\pm{\sinh(2\beta)\over\sqrt{1+\alpha^{2}}}\Bigg{\\}}$
(40)

Now we know that $\alpha\in\mathbb{R}$ and hence the closed form expression
for the resonance condition can be derived from equating ${\rm
Im}(\alpha^{2})=0$ and ${\rm Re}(\alpha^{2})\geq 0$.

From (40), we can derive the expression for $\alpha^{2}$ as

$\alpha^{2}={\sinh^{2}(2\beta)\over[\cosh(2\beta)\cosh(2J)+\cos
k\sinh(2J)]^{2}}-1$ (41)

and the subsequent closed form resonance conditions for persistent
oscillations immediately follow

$\displaystyle{\rm
Im}\\!\left({\sinh^{2}(2\beta)\over[\cosh(2\beta)\cosh(2J)+\cos
k\sinh(2J)]^{2}}\right)$ $\displaystyle=0$ (42) $\displaystyle{\rm
Re}\\!\left({\sinh^{2}(2\beta)\over[\cosh(2\beta)\cosh(2J)+\cos
k\sinh(2J)]^{2}}\right)$ $\displaystyle\geq 1\ .$ (43)

Note that both conditions in (42) and (43) need to be satisfied simultaneously
in order to get the correct phase diagram, which matches the phase diagram
obtained from TRG and analytically continued Onsager solution (see Fig. 1).

Steady state – We now obtain an expression for the steady state, where each
$(k,-k)$ mode pair is in an eigenstate of $\Theta^{\vphantom{*}}_{k}$. With
$\Omega^{\vphantom{*}}_{\mu}(k)\equiv
d^{\vphantom{*}}_{\mu}(k)/d^{\vphantom{*}}_{0}(k)$ for $\mu\in\\{x,y,z\\}$, we
have

$\Theta^{\vphantom{*}}_{k}=2d^{\vphantom{*}}_{0}(k)\begin{pmatrix}1+\Omega^{\vphantom{*}}_{z}(k)&-i\Omega^{\vphantom{*}}_{y}(k)\\\
i\Omega^{\vphantom{*}}_{y}(k)&1-\Omega^{\vphantom{*}}_{z}(k)\end{pmatrix}\quad.$
(44)

As $\Theta^{\vphantom{*}}_{k}$ is in general non-Hermitian, its right and left
eigenvectors are not necessarily related by complex conjugation, and are given
by

$\begin{split}{|\hskip-1.2pt|\,R^{\vphantom{*}}_{\pm}(k)\,\rangle\hskip-2.2pt\rangle}&=\begin{pmatrix}\mu^{\vphantom{*}}_{\pm}(k)\\\
\mu^{\vphantom{*}}_{\mp}(k)\end{pmatrix}\\\
{\langle\hskip-2.2pt\langle\,L^{\vphantom{*}}_{\pm}(k)\,|\hskip-1.2pt|}&=\pm{\cal
N}^{\vphantom{*}}_{k}\begin{pmatrix}\mu^{\vphantom{*}}_{\pm}(k)&,&-\mu^{\vphantom{*}}_{\mp}(k)\end{pmatrix}\quad,\end{split}$
(45)

where

$\mu^{\vphantom{*}}_{\pm}(k)=\sqrt{\Omega^{\vphantom{*}}_{+}(k)}\pm\sqrt{\Omega^{\vphantom{*}}_{-}(k)}$
(46)

and ${\cal
N}^{\vphantom{*}}_{k}=1/4\sqrt{\Omega^{\vphantom{*}}_{+}(k)\,\Omega^{\vphantom{*}}_{-}(k)}$.
These states are normalized so that
${\langle\hskip-2.2pt\langle\,L^{\vphantom{*}}_{a}(k)\,|\hskip-1.2pt|\,R^{\vphantom{*}}_{b}(k)\,\rangle\hskip-2.2pt\rangle}=\delta^{\vphantom{*}}_{ab}$
(with no complex conjugation of the left eigenvector).

In the $t\to\infty$ limit, and at each wavevector $k\in(0,\pi)$, one of these
states is selected – the one corresponding to the greater value of
$\big{|}\lambda^{\vphantom{*}}_{\pm}(k)\big{|}$. The surviving state’s
wavefunction is given by the appropriately normalized right eigenvector, and
the asymptotic state is thus of the BCS form,

${|\,\Psi(t\to\infty)\,\rangle}={\prod_{k>0}}^{\prime}{\cal
C}^{\vphantom{*}}_{k}\,\Big{[}\mu^{\vphantom{*}}_{\mp}(k)+\mu^{\vphantom{*}}_{\pm}(k)\,c^{\dagger}_{k}c^{\dagger}_{-k}\Big{]}\,{|\,0\,\rangle}\quad,$
(47)

where ${\cal C}_{k}=1/\sqrt{2\big{(}|\Omega_{k,+}|+|\Omega_{k,-}|\big{)}}$,
and where ${|\,0\,\rangle}$ is the Fock space vacuum, equivalent to the state
${|\,\\!\\!\downarrow\downarrow\cdots\downarrow\,\rangle}$ for all the
$k$-space Anderson pseudospins.

## References

  * Matsubara (1955) T. Matsubara, Prog. Theor. Phys. 14, 351 (1955).
  * Skinner  _et al._ (2019) B. Skinner, J. Ruhman, and A. Nahum, Phys. Rev. X 9, 031009 (2019).
  * Li  _et al._ (2018) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 98, 205136 (2018).
  * Li  _et al._ (2019) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 100, 134306 (2019).
  * Chan  _et al._ (2019) A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Phys. Rev. B 99, 224307 (2019).
  * Cao  _et al._ (2019) X. Cao, A. Tilloy, and A. D. Luca, SciPost Phys. 7, 24 (2019).
  * Choi  _et al._ (2020) S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Phys. Rev. Lett. 125, 030505 (2020).
  * Gullans and Huse (2019) M. J. Gullans and D. A. Huse, arXiv e-prints , arXiv:1905.05195 (2019), arXiv:1905.05195 [quant-ph] .
  * Gullans and Huse (2020) M. J. Gullans and D. A. Huse, Phys. Rev. Lett. 125, 070606 (2020).
  * Zabalo  _et al._ (2020) A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, Phys. Rev. B 101, 060301 (2020).
  * Jian  _et al._ (2019) C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. W. Ludwig, arXiv e-prints , arXiv:1908.08051 (2019), arXiv:1908.08051 [cond-mat.stat-mech] .
  * Bao  _et al._ (2020) Y. Bao, S. Choi, and E. Altman, Phys. Rev. B 101, 104301 (2020).
  * Tang and Zhu (2020) Q. Tang and W. Zhu, Phys. Rev. Research 2, 013022 (2020).
  * Fan  _et al._ (2020) R. Fan, S. Vijay, A. Vishwanath, and Y.-Z. You, arXiv e-prints , arXiv:2002.12385 (2020), arXiv:2002.12385 [cond-mat.stat-mech] .
  * Li and Fisher (2020) Y. Li and M. P. A. Fisher, arXiv e-prints , arXiv:2007.03822 (2020), arXiv:2007.03822 [quant-ph] .
  * Lavasani  _et al._ (2020) A. Lavasani, Y. Alavirad, and M. Barkeshli, arXiv e-prints , arXiv:2004.07243 (2020), arXiv:2004.07243 [quant-ph] .
  * Sang and Hsieh (2020) S. Sang and T. H. Hsieh, arXiv e-prints , arXiv:2004.09509 (2020), arXiv:2004.09509 [cond-mat.stat-mech] .
  * Ippoliti  _et al._ (2020) M. Ippoliti, M. J. Gullans, S. Gopalakrishnan, D. A. Huse, and V. Khemani, arXiv preprint arXiv:2004.09560 (2020).
  * Alberton  _et al._ (2020) O. Alberton, M. Buchhold, and S. Diehl, arXiv e-prints , arXiv:2005.09722 (2020), arXiv:2005.09722 [cond-mat.stat-mech] .
  * Fidkowski  _et al._ (2020) L. Fidkowski, J. Haah, and M. B. Hastings, arXiv preprint arXiv:2008.10611 (2020).
  * Suzuki and Fisher (1971) M. Suzuki and M. E. Fisher, J. Math. Phys. 12, 235 (1971).
  * Yang and Lee (1952) C. N. Yang and T. D. Lee, Phys. Rev. 87, 404 (1952).
  * Lee and Yang (1952) T. D. Lee and C. N. Yang, Phys. Rev. 87, 410 (1952).
  * Griffiths (1969) R. B. Griffiths, Phys. Rev. Lett. 23, 17 (1969).
  * Lee and Chan (2014) T. E. Lee and C.-K. Chan, Phys. Rev. X 4, 041001 (2014).
  * Biella and Schiró (2020) A. Biella and M. Schiró, arXiv preprint arXiv:2011.11620 (2020).
  * Gopalakrishnan and Gullans (2020) S. Gopalakrishnan and M. J. Gullans, arXiv preprint arXiv:2012.01435 (2020).
  * Nahum  _et al._ (2020) A. Nahum, S. Roy, B. Skinner, and J. Ruhman, arXiv preprint arXiv:2009.11311 (2020).
  * Jian  _et al._ (2021) S.-K. Jian, Z.-C. Yang, Z. Bi, and X. Chen, arXiv preprint arXiv:2101.04115 (2021).
  * (30) F. Beichert, C. A. Hooley, R. Moessner, and V. Oganesyan, arXiv preprint 1304.6314 (2013).
  * Khemani  _et al._ (2019) V. Khemani, R. Moessner, and S. Sondhi, arXiv preprint arXiv:1910.10745 (2019).
  * Watanabe and Oshikawa (2015) H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015).
  * Chan  _et al._ (2015) C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys. Rev. A 91, 051601 (2015).
  * Lazarides  _et al._ (2020) A. Lazarides, S. Roy, F. Piazza, and R. Moessner, Phys. Rev. Research 2, 022002 (2020).
  * Ashida  _et al._ (2020) Y. Ashida, Z. Gong, and M. Ueda, arXiv preprint arXiv:2006.01837 (2020).
  * Bertini  _et al._ (2019) B. Bertini, P. Kos, and T. Prosen, Phys. Rev. Lett. 123, 210601 (2019).
  * Gopalakrishnan and Lamacraft (2019) S. Gopalakrishnan and A. Lamacraft, Phys. Rev. B 100, 064309 (2019).
  * Foss-Feig  _et al._ (2020) M. Foss-Feig, D. Hayes, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, J. M. Pino, and A. C. Potter, arXiv preprint arXiv:2005.03023 (2020).
  * Garratt and Chalker (2020a) S. J. Garratt and J. Chalker, arXiv preprint arXiv:2008.01697 (2020a).
  * Lerose  _et al._ (2020) A. Lerose, M. Sonner, and D. A. Abanin, arXiv preprint arXiv:2009.10105 (2020).
  * Garratt and Chalker (2020b) S. Garratt and J. Chalker, arXiv preprint arXiv:2012.11580 (2020b).
  * Sonner  _et al._ (2020) M. Sonner, A. Lerose, and D. A. Abanin, arXiv preprint arXiv:2012.00777 (2020).
  * Ippoliti and Khemani (2020) M. Ippoliti and V. Khemani, arXiv preprint arXiv:2010.15840 (2020).
  * Lu and Grover (2021) T.-C. Lu and T. Grover, arXiv preprint arXiv:2103.06356 (2021).
  * Ippoliti  _et al._ (2021) M. Ippoliti, T. Rakovszky, and V. Khemani, arXiv preprint arXiv:2103.06873 (2021).
  * Brandner  _et al._ (2017) K. Brandner, V. F. Maisi, J. P. Pekola, J. P. Garrahan, and C. Flindt, Phys. Rev. Lett. 118, 180601 (2017).
  * García-Saez and Wei (2015) A. García-Saez and T.-C. Wei, Phys. Rev. B 92, 125132 (2015).
  * Levin and Nave (2007) M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 (2007).
  * Hertz (1976) J. A. Hertz, Phys. Rev. B 14, 1165 (1976).
  * Brittin  _et al._ (1964) W. E. Brittin, A. O. Barut, and L. Marshall,  _Statistical Physics, Weak Interactions, Field Theory_, Lectures in Theoretical Physics, Vol. 7C (University of Colorado Press, 1964).
  * Stoudenmire and White (2010) E. Stoudenmire and S. R. White, New Journal of Physics 12, 055026 (2010).
  * Onsager (1944) L. Onsager, Phys. Rev. 65, 117 (1944).
  * van Saarloos and Kurtze (1984) W. van Saarloos and D. A. Kurtze, Journal of Physics A: Mathematical and General 17, 1301 (1984).
  * (54) S. Basu et al., in preparation.
  * Chaikin  _et al._ (1995) P. M. Chaikin, T. C. Lubensky, and T. A. Witten,  _Principles of Condensed Matter Physics_ , Vol. 10 (Cambridge university press Cambridge, 1995).