# Real-time dynamics of a critical Resonating Valence Bond spin liquid

Ravi Teja Ponnaganti  Laboratoire de Physique Théorique, C.N.R.S. and
Université de Toulouse, 31062 Toulouse, France  Matthieu Mambrini  Laboratoire
de Physique Théorique, C.N.R.S. and Université de Toulouse, 31062 Toulouse,
France  Didier Poilblanc  Laboratoire de Physique Théorique, C.N.R.S. and
Université de Toulouse, 31062 Toulouse, France

###### Abstract

Implementation of the hardcore-dimer Hilbert space in cold Rydberg-atom
simulators opens a new route of investigating real-time dynamics of dimer
liquids under Hamiltonian quench. Here, we consider an initial Resonating
Valence Bond (RVB) state on the square lattice realizing a critical Coulomb
phase with algebraic and dipolar correlations. Using its representation as a
special point of a broad manifold of SU($2$)-symmetric, translationally
invariant, Projected Entangled Pair States (PEPS), we compute its non-
equilibrium dynamics upon turning on inter-site Heisenberg interactions. We
show that projecting the time-evolution onto the PEPS manifold remains
accurate at small time scales. We also find that the state evolves within a
PEPS sub-manifold characterized by a U($1$) gauge symmetry, suggesting that
the Coulomb phase is stable under such unitary evolution.

##  I Introduction

The search for spin liquids in condensed matter materials is a very rapidly
developing area of quantum magnetism [1]. Spin liquids, as the Resonating
Valence Bond (RVB) state proposed by Anderson [2], are prototypical states of
matter showing no symmetry breaking, even down to zero temperature, due to
enhanced zero-point quantum fluctuations. Their highly entangled nature leads
to unique physical aspects [3, 4], such as emerging non-local excitations,
topological properties, etc…

Interestingly, spin liquids have also become accessible to other experimental
set-ups based on ultracold atoms loaded into a two-dimensional optical lattice
[5] or in two-dimensional Rydberg arrays [6, 7] offering an alternative route
to emulate quantum simulators [8, 9], in addition to condensed matter
superconducting circuits [10].

The rapid progress in cold atom experiment setups call for new efficient
theoretical tools to investigate non-equilibrium dynamics of isolated pure
quantum systems. Apart from one-dimension for which efficient techniques exist
[11, 12], computing the non-equilibrium dynamics that follows a quantum quench
[13] is a tedious task in two-dimensional quantum spin systems [14, 15]. Here
we shall investigate the non-equilibrium dynamics of the RVB state, the most
paradigmatic example of a spin liquid. For simplicity, we shall consider the
case of the square lattice. A pictorial representation of the nearest-neighbor
RVB state is shown on Fig. 1(a), consisting only of resonating nearest-
neighbor (NN) valence bond configurations [2]. Interestingly, it has been
shown that the NN RVB state exhibits critical dimer-dimer correlations
connected to a local U(1) gauge symmetry and characteristic of a Coulomb phase
[16, 17]. One of the main goals of this work is to investigate the stability
of the Coulomb phase following a quantum quench as well as the dynamics of its
entanglement.

For such a purpose, through out this paper, we shall use the tensor network
formalism using a variational Projected Entangled Pair State (PEPS) ansatz of
the time-evolving many-body wave function [18]. This procedure is, in spirit,
similar to ref. [8] where the preparation dynamics of a Rydberg quantum
simulator is approximated by projecting it on a tensor network manifold. In
our case, the PEPS is defined by a unique time-dependent tensor ${\cal A}(t)$
placed on all the sites of a two-dimensional square lattice. Using Penrose
graphical representation [19], the on-site tensor and its corresponding tensor
network are shown in Fig. 2(a) and (b). Hence ${\cal A}(t)$ encodes locally,
at all times, the coefficients of the many-body wave function
$|\Psi(t)\big{>}$ in the exponentially-large $S_{z}$-basis
$\\{\sigma_{1},\sigma_{2},\cdots,\sigma_{N}\\}$, $\sigma_{i}=\pm 1/2$, as
shown in Fig. 2(b). Note that the entanglement (which grows with time) is
controlled by the (virtual) bond dimension $D$ of the tensor ${\cal A}(t)$.
For simplicity we shall also take the limit of an infinite system,
$N\rightarrow\infty$, using the infinite-PEPS (iPEPS) framework [20].

Figure 1:  (a) A typical valence bond configuration of the NN RVB state. (b) A
valence bond configuration including longer range singlets. (c) A typical
configuration under time-evolution of the NN RVB state. Blue (red) dimers are
singlet bonds build from two NN virtual spin-1/2 (spin-1) of the PEPS. The
characteristic U($1$) gauge symmetry of the Coulomb phase is broken in state
(b).  Figure 2:  (a) The (time-dependent) site tensor (blue dot) contains 4
virtual bonds (in black) of dimension $D$ and a physical leg (in red) of
dimension $d=2$ (spin-1/2 or ”qubit”). (b) The iPEPS is obtained from an
infinite square-lattice array of site tensors contracted over all virtual
indices.

The paper is organized as follow; First, one describes the set-up of the
global quench in Section II: A description of the physical nature of the
initial spin liquid state and a summary of the standard Trotter-Suzuki
procedure used to obtain the time-evolution are given in Subsections II.1 and
II.2, respectively. We detail the implementation of the simple update (SU)
framework keeping track of the space-group and spin-SU($2$) symmetries of the
initial state and of the Hamiltonian in Section III. This procedure involves
the determination of the optimal virtual space described in Subsection III.2
and a careful gauge fixing implementation described in Subsection III.3.
Details on the recovery of the exact point group and SU($2$) symmetries at
each Trotter step is left to Appendix A. Results are provided and discussed in
Section IV: First, the SU singular value spectrum is analyzed as a function of
propagating time in Subsection IV.1. Then, we propose the reversal probability
(Lodschmidt echo) as an interesting tool to measure truncation errors in
Subsection IV.2. The relevance of the finite-$D$ ansatz is controlled by
investigating the energy conservation in Subsection IV.3. Finally, the
critical property of the time-evolved state is investigated via the scalings
of the correlation length and the entanglement entropy of the boundary state
in Subsection IV.4. A summary is given and possible further developments are
discussed in Section V. Additional material is provided in two other
Appendices. The new algorithms specific to the factorisation of complex
symmetric tensors appearing in this work are described in Appendix B. Other
specificities of the tensor contraction algorithm are provided in Appendix C.

##  II Quench protocol

###  II.1 Resonating Valence Bond States and PEPS representations

Let us first start by describing the simple system setup we have considered.
Our initial quantum state $|\Psi_{0}\big{>}$ is a Resonating Valence Bond
(RVB) spin liquid on an infinite square lattice and a global Hamiltonian
quench is assumed, at time $t=0$, by turning on the antiferromagnetic nearest-
neighbor (NN) Heisenberg Hamiltonian,

${\cal H}(t)=\begin{cases}\;0,&\text{for }t\leq 0\\\ \;H=\sum_{\langle
x,y\rangle}H_{xy},&\text{for }t>0\end{cases}$ (1)

where

$H_{xy}=J{\bf S}_{x}\cdot{\bf S}_{y}.$ (2)

The NN RVB state, consisting of resonating NN singlets as shown on Fig. 1(a),
is in fact a special point of an extended one-dimensional RVB family [4]
including longer-range valence (singlet) bonds (see Fig. 1(b)). This RVB
family is conveniently represented by a simple PEPS manifold spanned by two
single-site tensors with full lattice ($C_{4v}$) and spin-rotation (SU($2$))
symmetries; the four virtual legs have virtual space ${\cal
V}_{0}=0\oplus\frac{1}{2}$ and are contracted, whereas the physical legs
correspond to the spins in the lattice [3, 21] (see Fig. 2). Hence, one can
tune the initial state by simply varying the ratio $\lambda_{2}/\lambda_{1}$
of the coefficients of the on-site tensor,

${\cal A}(0)=\lambda_{1}{\cal A}_{1}+\lambda_{2}{\cal A}_{2}\,,$ (3)

where the tensor ${\cal A}_{1}$ defines the NN RVB state and the tensor ${\cal
A}_{2}$ induces longer-range singlets by “teleportation”. These tensors simply
differ by the occupation of the spin-$0$ and spin-$1/2$ states on the four
virtual bonds, $n_{\rm occ}=\\{3,1\\}$ and $n_{\rm occ}=\\{1,3\\}$,
respectively. Recent work [4, 22] suggested that topological order appears
whenever longer-range bonds are present (pictorially shown in Fig. 1(b)), i.e.
$\lambda_{2}\neq 0$, breaking the U($1$) gauge symmetry to $\mathbb{Z}_{2}$.
Here, we shall take advantage of the small bond dimension $D_{0}=3$, and of
the full symmetries of our initial state and of the Hamiltonian to study the
time evolution over a small time interval.

###  II.2 Time evolution

During its time evolution the RVB state is expected to preserve its global
singlet character (SU(2) rotation symmetry) and the full lattice symmetry
(C4v). In order to compute, for $t>0$,

$|\Psi(t)\big{>}=\exp{(-iHt)}|\Psi_{0}\big{>}$ (4)

we have used a simple update (SU) method [23] which can be implemented in a
way that preserves all the symmetries under consideration. Hereafter time $t$
will be measured in unit of $1/J$. Here the time-evolved state is defined by a
unique on-site complex PEPS tensor of bond dimension $D>D_{0}$ (see Fig. 2(a))
expanded in a fully lattice-C4v/spin-SU(2) symmetric (real) tensor basis
$\\{T_{a}\\}$,

${\cal A}(t)=\sum_{a=1}^{M}\mu_{a}(t)\,T_{a},$ (5)

with $\mu_{a}(t)\in\mathbb{C}$. This is a simple generalization [18] of the
symmetric PEPS construction of the initial RVB state. Since entanglement grows
with time, it is necessary to include new virtual states, hence increasing the
bond dimension and the number of elementary tensors. Later, we will show that,
for the NN RVB $\lambda_{2}=0$ (or $\lambda_{2}$ small enough) it is
sufficient to consider a virtual space $V=0\oplus\frac{1}{2}\oplus 1$, i.e. a
bond dimension $D=6$, to describe the time-evolution for $t\lesssim 1$ (see
Table 1).

$D$ | ${\cal V}$ | Cs | Cs/U(1) | C4v | C4v/U(1)
---|---|---|---|---|---
3 | $0\oplus\frac{1}{2}$ | 7 | 4 | 2 | 1
6 | $0\oplus\frac{1}{2}\oplus 1$ | 41 | 30 | 11 | 8
Table 1: Number of SU($2$)-symmetric local tensors of bond dimension $D=3$ or
$D=6$, virtual space $\cal V$, Cs or C4v point-group symmetry, and
with/without U($1$) gauge symmetry (see text). The number of $C_{s}$ symmetric
tensors gives the number of degrees of freedom available at every time sub-
step. We explicitly find that the non-U($1$) symmetric tensors keep a
vanishing weight under time evolution. Figure 3:  Simple update scheme: (a)
the site tensor (blue sphere) is split by SVD to isolate the active bond
(black leg on the red sphere); (b) The complex symmetric gate ${\cal
G}_{xy}^{\alpha}(\tau)$ applied to the two-site bond (referred in the text as
“SU matrix”) is split using a Autonne-Takagi factorization [24, 25] (further
details are given in Appendix B); (c) the two sides are used to reconstruct
the new site tensors (blue sphere with updated blue leg). Note that for step
(a) the SVD spectrum $S$ (diagonal matrix) is absorbed entirely in the active
bond tensor (depicted in red), i.e. $USV^{\dagger}=U\left(SV^{\dagger}\right)$
while for step (b) the Autonne-Takagi spectrum $S$ is symmetrically absorbed
on both sides to fulfill the reflexion symmetry (depicted as a gray plane),
i.e. $USU^{T}=\left(U\sqrt{S}\right)\left(U\sqrt{S}\right)^{T}$.

Once the virtual space is defined, in order to obtain the time-evolution of
the coefficients $\\{\mu_{a}(t)\\}$ we used a standard Trotter-Suzuki (TS)
decomposition [26] of the unitary time-evolution operator

$\exp{(-iHt)}=\prod_{1}^{N_{\tau}}\exp{(-iH\tau)},$

where $\tau=t/N_{\tau}$ is a small time step (such that $\tau\ll 1$). The
Heisenberg Hamiltonian is then split into four parts,

$H=H^{A}+H^{B}+H^{C}+H^{D}\,,$ (6)

each acting on one of the four staggered configurations ${\cal C}^{\alpha}$ of
disconnected horizontal or vertical bonds labelled by $\alpha=A,B,C,D$. The
action of the elementary time-evolution operator $\exp{(-iH\tau)}$ can then be
approximated by the successive actions of four unitary gates,

$\displaystyle{\cal G}^{\alpha}(\tau)$ $\displaystyle=$
$\displaystyle\exp{(-iH^{\alpha}\tau)}$ (7) $\displaystyle=$
$\displaystyle\prod_{\rm\langle x,y\rangle\in
C_{\alpha}}\exp{(-iH_{xy}^{\alpha}\tau)},$

involving the standard systematic TS error vanishing in the limit
$\tau\rightarrow 0$. Depending on the method, the update of the coefficients
$\mu_{a}(t+\tau)$ under the action of all disconnected gates

${\cal G}_{xy}^{\alpha}(\tau)=\exp{(-iH_{xy}^{\alpha}\tau)}$ (8)

may be obtained locally (SU method) or take into account the environment
around each of the disconnected $(x,y)$ bonds. We shall here focus on the
simple update method, which we describe below, and will report on a time
dependent variational optimization method (involving the environment)
elsewhere. Note that, at every sub-step defined by the action of all
disconnected gates $\alpha$, the lattice $C_{4v}$ point-group symmetry is
broken down to $C_{s}$ involving only the reflection w.r.t. the direction of
the bonds $(x,y)$. Therefore, the updated one-site tensor has a basis
decomposition ( 5) involving a larger set of $C_{s}$-symmetric tensors (see
Table 1). Only after a full step of four sub-steps $\alpha=A,B,C,D$ is the
point-group $C_{4v}$ symmetry (approximately) restored.

##  III Simple update numerical algorithm

###  III.1 Algorithmic steps

The main lines of the Simple Update (SU) scheme we use can be summarized in
the following steps :

  1. 1.

The tensors at sites $x$ and $y$ are first split by singular value
decomposition (SVD) to isolate the active bond on which the complex symmetric
gate ${\cal G}_{xy}^{\alpha}(\tau)$ is applied (see Fig. 3(a)).

  2. 2.

The resulting two-site symmetric complex matrix (“SU matrix” of Fig. 3(b)) is
decomposed using an Autonne-Takagi factorization [24, 25] and Appendix B for
details.

  3. 3.

The two sides are used to reconstruct the first update of the tensors at site
$x$ and $y$ shown in Fig. 3(c).

  4. 4.

This procedure is then repeated for the other three bonds connected to site
$x$.

  5. 5.

At last, gauge fixing and projection of ${\cal A}(t+\tau)$ onto the symmetric
basis $\\{T_{a}\\}$ enables us to obtain the new set of coefficients
$\mu_{a}(t+\tau)$.

However, the actual implementation of this method (and especially point 5.)
deserves a particular attention in the context of SU(2) invariant tensors.
More specifically, we stress two important issues in the course of SU : (i)
Re-defining an identical gauge convention between the four bonds after the
sequence of independent factorizations and (ii) Enforcing a fully $C_{4v}$ and
SU(2) evolution at each time step. These points are discussed in detail in
section III.3.

###  III.2 Determination of the relevant virtual space

It is particularly interesting to first have a close look at the first
application of the ${\cal G}_{xy}^{A}(\tau)$ gate onto the initial (general)
RVB state. The result of this action will guide us to select the relevant
choice of the virtual space needed to approximate the time evolution. Fig. 4
shows the 12 singular values of the Autonne-takagi factorization of the matrix
in Fig. 3(b), grouped in SU(2) multiplets, as a function of
$\theta=\tan^{-1}{(|\lambda_{2}|/\lambda_{1})}$. Remarkably, at
$\lambda_{2}=0$ only 6 singular values are non-zero corresponding to spin-$0$,
spin-$1/2$ and spin-$1$ multiplets. This provides support for using the
virtual space ${\cal V}=0\oplus\frac{1}{2}\oplus 1$ to describe time evolution
at finite time $t$. Such a PEPS can be pictorially represented in the RVB
language as in Fig. 1(c). The range of validity in time of this approximation
will be discussed later on. Note that, when turning on $\lambda_{2}$ the
additional spin-$\frac{1}{2}$ and spin-$\frac{3}{2}$ multiplets acquire some
weights which increase with increasing $\lambda_{2}$ and, hence, could not be
neglected anymore.

Figure 4:  Singular values of the Autonne-Takagi factorization of the gate of
Fig. 3(b) computed for $\tau=0.025$ as a function of $\theta$ with
$\lambda_{1}=\cos(\theta)$ and $\lambda_{2}=\sin(\theta)$.  Figure 5: Gauge
fixing and symmetries in SU : Starting from a fully SU(2) and $C_{4v}$
symmetric site tensor (a) the four legs are updated using the SU method (see
Fig. 3) using ${\cal G}^{\alpha}(\tau)$ gate. (b) Due to gauge freedom the
basis used is generically different on the four virtual legs (depicted as
distinct colors). (c-e) A uniform fixed gauge is obtained in a 3-step process
in which the norm of the 2 reflection dissymetries (c-d) and rotation
dissymetry (e) are minimized by a complete enumeration of all possible basis
(see Eq.(9)). The resulting tensor (f) is expressed in a uniform gauge and is
approximately $C_{4v}$ symmetric (up to $\tau^{2}$ corrections). (h) The exact
$C_{4v}$ symmetry is explicitly restored. (i) Projection into the SU(2)
symmetric tensor basis leads to a SU(2) and $C_{4v}$ symmetric tensor suitable
for the next $\tau$-step evolution.

###  III.3 Gauge fixing

In this Subsection we describe the implementation of the tensor symmetrisation
at every time step (algorithmic step 5. mentioned in III.1). Readers
interested mostly by the physical problem may jump directly to Section IV.

According to Trotter-Suzuki decomposition (see Eqs. (6),(7),(8) and Fig. 3),
starting from a fully symmetric tensor (Fig. 5(a)), four Autonne-Takagi
transformations are required at every time step. Since tensors are SU(2)
invariant, the singular values spectra splits into degenerate sectors
corresponding to the various multiplets (see e.g. Fig. 4). As a direct
consequence, the unitary transformation involved in the factorization is not
unique. More precisely, any block SU(2) rotation acting on each multiplet
subspace ($S=0,S=1/2,\ldots$) leads to an equally valid decomposition. This
results in a continuous gauge freedom and generically a mismatch of the basis
used on the four virtual legs (see Fig. 5(b)), which prevents any further
computation.

If the gauge is uniformly fixed on the four virtual legs, the updated tensor
remains symmetric under $C_{4v}$ transformations up to $\tau^{2}$ corrections,
which corresponds to the lowest order where non-commutativity effects occurs
in the Trotter-Suzuki decomposition. This fact leads to a natural criterion to
fix a uniform gauge on all four virtual legs, optimizing the tensor point
group symmetry. However minimizing the tensor dissymetry under a continuous
set of unitary transformations reveals intractable in practice.

A way to circumvent this problem is to apply to the tensor a multiplicative
random noise (controlled by its amplitude $\varepsilon$) before performing the
4-step update. This kind of noise breaks SU(2) symmetry but preserve the U(1)
symmetry related to charge conservation ($S_{z}$ in the context of SU(2)).
Hence the singular values spectra appearing in the various decompositions get
the structure of slightly splited SU(2) multiplets. Interestingly, in this
scheme, $\varepsilon$ determines the relative splitting of the noise-free
multiplet singular value. This makes it possible to clearly identify the
multiplets in the spectra, regardless of the amplitude of the considered
singular value.

But, even more importantly, disorder reduces drastically the gauge freedom to
a discrete set of transformations. Keeping in mind the properties of Autonne-
Takagi transformation detailed in Appendix B.1 (in particular eq. 14), the
gauge transformation relating two equivalent Autonne-Takagi factorizations in
the case of a non-degenerate spectrum is just a diagonal matrix with $\pm 1$
entries. On top of this, the multiplicative noise can cause a reordering of
states inside each multiplet in a way that does not match the canonical order
(e.g $S_{z}=+1,0,-1$ for $S=1$). Hence the most general gauge transformation
summarizes as :

${\cal U}=\begin{pmatrix}\pm 1&&\\\ &\ddots&\\\ &&\pm
1\end{pmatrix}\begin{pmatrix}{\cal P}_{v_{1}}&&\\\ &\ddots&\\\ &&{\cal
P}_{v_{n}}\end{pmatrix},$ (9)

where ${\cal V}=\oplus_{i=1}^{n}v_{i}$ and ${\cal P}_{v_{i}}$ is a permutation
matrix in the $v_{i}$ subspace ($n$ stands here for the total number of
species). Hence the maximal number of distinct transformation is
$2^{D}\prod_{i=1}^{n}{\text{Dim}(v_{i})}!$ which is small enough in practical
applications to allow direct enumeration in the course of tensor dissymetry
minimization. The dissymetry over the two reflections and $\pi/2$-rotation is
described if Fig. 5 (c-e).

Once a uniform consistent gauge is fixed for the four virtual legs, the next
task is to project the evolved tensor into the symmetric basis $\\{T_{a}\\}$.
In this process the resulting tensor becomes fully symmetric under $C_{4v}$
and SU(2) (see Fig. 5 (h-i)). We demonstrate in Appendix A that, for small
enough $\tau$, the modification caused to the tensor by the symmetrization is
small compared to the increment of the tensor due to the time evolution
itself. Fig. 6 summarizes the time evolution of the complex components
$\mu_{a}(t)$ of ${\cal A}(t)$ expressed in the symmetric basis, for
$t\in[0,2]$, starting from the NN RVB state.

Figure 6:  Time evolution of ${\cal A}(t)$ complex components $\mu_{i}$ (see
Eq. 5) using $\tau=0.025$ in the time range $t\in[0,2]$, starting from the
initial tensor ${\cal A}(0)={\cal A}_{1}$ (NN RVB state). In each panel $\mu$
is displayed in the complex plane (black dots) as well as its real part (blue
spheres) and imaginary part (red spheres). The components $\mu_{2},\mu_{3}$
and $\mu_{4}$ of the non-U($1$) symmetric tensors are found to be identically
zero and, hence, are not displayed.

##  IV Results

In this section we shall focus on the NN RVB ($\lambda_{2}=0$) as an initial
state for which short-time evolution with the $D=6$ PEPS ansatz has been
justified above.

###  IV.1 Singular value spectrum

Since entanglement quickly grows with time, we expect that more virtual states
(i.e. larger $D$) may become necessary as time goes on. To control the
validity of our fixed $D=6$ approximation, we have examined, at each step
after applying the first gate ${\cal G}_{xy}^{A}(\tau)$ on $|\Psi(t)\big{>}$,
the singular value spectrum of the Autonne-Takagi factorization of the SU
matrix (defined in Fig. 3(b)), as a function of time $t$. We see in Fig. 7(a)
that the 3 multiplets of largest weights always stay well separated from the
rest of the spectrum. As expected, their weights tend to become equal, in
order to saturate the maximum available entanglement entropy $\ln{D}=\ln{6}$
per site, as shown in Fig. 7(b). This is obviously an artefact of the
truncation into the $D=6$ virtual space which, although involves only a small
error at every step, leads to a significant cumulative error when $t\sim 1$.
Hence further tests are needed to establish the maximum range of validity of
our approximation.

Figure 7:  (a) Singular values $\lambda_{i}$ of the Autonne-Takagi
factorization of the SU gate matrix represented in Fig. 3(b) and computed for
$\tau=0.025$, as a function of time. (b) Entropy per bond as a function of
time captured by the representations $0\oplus 1/2\oplus 1$ used in the SU
procedure.

###  IV.2 Fidelities and Loschmidt echo

In order to establish the ultimate maximum time above which our procedure
breaks down, we have performed the following “time reversal” procedure; (i)
time evolution is performed with $\exp{(-iHt)}$ from $t=0^{+}$ to $t=t_{R}$;
(ii) then, at $t=t_{R}$ time is “reversed”, $t\rightarrow-t$, i.e. time
evolution is performed with $\exp{(iHt)}$. In the case of an exact unitary
evolution, one should recover exactly the initial NN RVB state at time
$t=t_{R}$. In other words, (the modulus of the overlap) $|\big{<}{\cal
A}(t)|{\cal A}(0)\big{>}|$ of the time-dependent tensor with the initial
tensor (defining the $\lambda_{2}=0$ NN RVB state) should come back to 1 at
$t=2t_{R}$, see Fig. 8(a). Upon increasing $t_{R}$, due to various cumulative
errors, the evolution will stop being perfectly unitary after some time
$t<t_{R}$. If this happens then the Loschmidt echo $|\big{<}{\cal
A}(2t_{R})|{\cal A}(0)\big{>}|$ will start deviating from 1. This is exactly
what is observed on Fig. 8(b) suggesting that our procedure breaks down above
some “upper bound” time scale $t_{\rm max}\simeq 1$. Several sources of error,
at different degrees, are responsible for this breakdown, including the
Trotter-Susuki decomposition, the SU scheme instead of a full update scheme,
the symmetrisation procedure and the restriction of a limited virtual space
($D=6$ here).

Figure 8:  (a) Overlap $|\big{<}{\cal A}(t)|{\cal A}(0)\big{>}|$ of the time-
dependent tensor with the initial tensor defining the NN RVB state
($\lambda_{2}=0$) as a function of time (blue dots). Red dots are used after
time reversal at $t=t_{R}$, going backwards in time. (b) Loschmidt echo ${\cal
L}(t_{R})=|\big{<}{\cal A}(t=2t_{R}-\tau)|{\cal A}(\tau)\big{>}|$ vs $t_{R}$.
Note that for convenience ${\cal L}(t_{R})$ is not evaluated between times $0$
and $2t_{R}$ but between times $\tau$ and $2t_{R}-\tau$ to avoid dimensional
jump from $D=3$ at time $0$ to $D=6$ at time $t>0$. In (a) and (b) the SU
method is used with $\tau=0.025$.

###  IV.3 Energy conservation

Under unitary time evolution, the energy – defined as the expectation value of
$H$ in the time-evolving state – should be conserved. However, even if the
simple update (symmetric) scheme succeeds (up to $t_{\rm max}$) to optimize
properly the local tensor, the $D=6$ PEPS may stop accounting properly for the
increase of entanglement above some intermediate time $t_{\rm max}^{\rm
D=6}<t_{\rm max}$. We believe a good criterion to estimate this intermediate
time scale is to examine the possible deviation of the energy w.r.t. its $t=0$
value. The latter is then computed using the SU tensors via a Corner Transfer
Matrix Renormalization Algorithm (CTMRG) iPEPS algorithm [27, 28, 29, 30].
Note that the energy of the NN RVB converges slowly with the CTMRG environment
dimension $\chi$ (due to its critical nature) so that several values of $\chi$
will be considered here. Fig. 9 shows the energy per site as a function of
time $t$ for $\chi$ ranging from $D^{2}$ to $3D^{2}$. At the smallest
$\chi=D^{2}$ value, one clearly sees a plateau at small time $t\lesssim 0.25$
and, then, a clear deviation from the initial $t=0$ energy. Upon increasing
$\chi$ the deviation seems to occur a bit sooner. We believe this is a sign
that the restriction to $D=6$ is no longer accurate when $t>t_{\rm max}^{\rm
D=6}\simeq 0.25$ and/or a full update scheme would become necessary. Notice
however the energy scale used in Fig. 9(b), showing that the energy deviation
still remains relatively small up to intermediate time, e.g. $0.6\%$ at
$t=0.2$.

###  IV.4 Critical behavior, maximum correlation length and central charge

U(1) gauge symmetry – The basis tensors (see Table 1) involved in the
construction of our local tensor in Eq. (5) can be grouped according to their
occupation number of each of the virtual states on the four bonds $n_{\rm
occ}=\\{n_{0},n_{1/2},n_{1}\\}$ (with $\sum_{\alpha}n_{\alpha}=4)$. For
example, the initial NN RVB state has $n_{\rm occ}=\\{3,1,0\\}$. The integer
$n_{1/2}$ defines a U($1$) quantum number characterizing a U($1$) gauge
symmetry. Because of the SU($2$) fusion rules, $n_{1/2}$ can be either $1$ or
$3$ defining two separate classes of U($1$)-symmetric PEPS. We argue here that
time evolution preserves the $n_{1/2}=1$ U(1) quantum number because the
action of a two-site gate can only change the spin quantum number of the
virtual state by 1 on the bond it is applied to. We have checked numerically
that this conjecture is correct, more precisely the coefficients of all the
tensors with $n_{1/2}=3$ remains identically zero, at all times. The preserved
U(1) symmetry suggests that the time evolved state remains critical, as a
mapping to a height representation [16, 17] would imply.

Figure 9:  Energy per site – ideally a constant of motion (horizontal dashed
lines) – versus time $t$, using a $D=6$ SU($2$)-symmetric iPEPS ansatz for
$\chi$ varying from $D^{2}$ to $3D^{2}$. The local tensor ${\cal A}$(t) is
obtained via a SU procedure and $\tau=0.005$. (a) Full time range. (b) Zoom of
the small time region.

We believe the U(1)-symmetry of the time-evolving spin liquid is not fine-
tuned in the sense that it is not restricted to the $D=6$ ansatz we are using.
Indeed, for any family of (singlet) SU(2)-symmetric PEPS with bond dimension
$D$ one can always isolate a large subset of PEPS with U(1)-gauge symmetry
(which should exhibit algebraic correlations). The (simple) argument is as
follows: each local tensor is a projection from ${\cal V}^{\otimes 4}$ into
the spin-1/2 subspace. From the SU(2) conservation rules, one can only have an
odd number, one or three, of virtual states carrying half-integer spins (which
could be different like $1/2$ and $3/2$). Hence the SU($2$) tensor basis can
be split into two sets, and from each of them one can build two separate
families of U($1$)-symmetric tensors. However, combining basis tensors of the
two sets, i.e. with different numbers of half-integer virtual spins, will
break the U(1)-gauge symmetry into $\mathbb{Z}_{2}$.

There is however some sort of fine-tuning (or ”protection”) in the set-up
itself: this is in fact the form of the quench Hamiltonian which i) is
SU(2)-symmetric and ii) involves only NN bonds. In the special case ii) one
can use the simple TS decomposition in terms of 2-site gates. From i) each
gate acting on a bond can only change the virtual spin $S$ into $S+1$ or $S-1$
and, hence, will not change the U(1)-gauge symmetry of the state. Such a
property will no longer be true if either i) or ii) is broken. In particular,
the gauge U($1$) symmetry will be spoiled by a small amount of disorder
introducing (small) violations of i) or ii).

Transfer matrix – To investigate further the expected critical nature of the
state we have computed the spectrum of the transfer matrix (TM) build from two
$T$ environment tensors [4]. The (modulus of the) leading eigenvalues are
shown in Fig. 10 for $t=0.15$ (a) and $t=0.5$ (b) (but similar results are
also found for other values of $t$), normalizing the spectrum such that the
largest eigenvalue is $\lambda_{1}=1$. The behavior of the spectra with
increasing $\chi$ suggests a vanishing gap $\lambda_{1}-\lambda_{2}$ between
the largest eigenvalue and the subleading one $\lambda_{2}$. The maximum
correlation length $\xi_{\rm max}=-1/\ln{(\lambda_{2}/\lambda_{1})}$ is then
expected to diverge with $\chi$, as confirmed in Fig. 11 showing a linear
behavior of $\xi_{\rm max}$ with $\chi$ and no sign of saturation. In contrast
to the gapless singlet ($g=1$) spectrum a gap is seen in the doublet
(spin-1/2) and triplet ($g=3$, spin-1) sectors. Note that the doublet spectrum
exhibits an extra two-fold degeneracy ($g=4$). Such results are consistent
with short-range spin-spin correlations, as discussed later.

Figure 10:  TM spectrum at a fixed value of $t=0.15$ (a) and $t=0.5$ (b),
computed for several values of $\chi$. Different symbols are used to highlight
particular levels, the subleading $g=1$ eigenvalue (red dot) and the largest
eigenvalues $\lambda_{\rm spinon}$ and $\lambda_{T}$ with degeneracy $g=4$
(green dot) and $g=3$ (black dot), respectively. A gapless continuous singlet
spectrum is expected in the $\chi\rightarrow\infty$ limit.

Central charge – From the CTMRG environment one can build the boundary MPS of
bond dimension $\chi$ and “physical” dimension $D^{2}$. The critical nature of
the bulk PEPS is reflected in the critical nature of the boundary chain
characterized by a central charge $c=1$ for all time $t$. We have confirmed
this feature by computing the Von Neumann entanglement entropy $S_{\rm vN}$ of
the MPS. The later is shown in Fig. 12 as a function of $\ln{(\xi_{\rm max})}$
for small $t$ values. Fitting the data as $S_{\rm
vN}\sim\frac{c}{6}\ln{(\xi_{\rm max})}$[31, 32] one obtains $c\simeq 1$
consistent with the expected result.

Figure 11:  Maximum correlation length $\xi_{\rm max}$ for several values of
time $t$ showing absence of saturation as a function of $\chi$.  Figure 12:
Entanglement entropy of the boundary MPS versus $\ln{(\xi_{\rm max})}$ for
different time $t=0,0.05,0.1,0.15,0.3$. The dotted line correspond to the
behavior expected for $c=1$.

###  IV.5 Spin-spin correlations

The previous results suggest that the time-evolving state bears similar
properties as the initial NN RVB state, although with rapidly growing
entanglement. First, the diverging correlation length corresponds to power-law
decaying dimer-dimer correlations. Secondly, spin-spin correlations are
expected to be short-range. This is indeed seen in Fig. 13(a) for $t=0.15$.
The corresponding spin-spin correlation length $\xi_{S}$ is shown in Fig.
13(b) as a function of time $t$. Interestingly, the values extracted from fits
of the spin-spin correlations match the values $\xi_{T}$ obtained from the TM
spectra (see e.g. Fig. 10 (a,b) for $t=0.15$ and $t=0.5$): considering the
leading spin-1 ($g=3$) eigenvalue $\lambda_{T}$, one gets
$\xi_{T}=-1/\ln{(\lambda_{T}/\lambda_{1})}\simeq\xi_{S}$. Fig. 13(b) also
shows the behavior of the spinon correlation length $\xi_{\rm
spinon}=-1/\ln{(\lambda_{\rm spinon}/\lambda_{1})}$ where $\lambda_{\rm
spinon}$ is the leading eigenvalue with $g=4$ corresponding to two degenerate
spin-doublets ($S=1/2$). Remarkably, we observe that all these correlation
lengths vary moderately as a function of time.

Figure 13:  (a) Spin-spin correlation versus distance for $t=0.15$ in semi-log
scale for several values of $\chi$. The linear fit corresponds to an
exponential decay. (b) Correlation lengths at $\chi=108$ versus time $t$:
$\xi_{S}$ extracted from fits of the spin-spin correlations; $\xi_{T}$ and
$\xi_{\rm spinon}$ extracted from the TM spectra shown in Fig. 10 (using same
colors of dots).

##  V Conclusions

In this work we have considered a simple quench setup to investigate the non-
equilibrium time dynamics of a genuine critical spin liquid i.e. realizing a
Coulomb phase. Our description uses the iPEPS algorithm to approximate the
time evolution under the application of a NN Heisenberg interaction, allowing
for the preservation of all symmetries, both the lattice symmetry (by
considering a unique $C_{4v}$-symmetric site tensor) and the spin-rotation
SU(2) symmetry directly encoded at the level of the site tensor. At every time
step, the inerrant breaking of the point group $C_{4v}$ symmetry due to the
successive application of four (non-commuting) two-site gates followed by SVD
truncations is repaired by subtle gauge transformations on the four tensor
virtual legs.

Although our procedure is accurate only at small time due to the limited
tractable bond dimension $D$, we argue that the observed stability of the
critical nature of the evolving state is valid at all times in the case of
nearest neighbor interactions. The robustness of the U(1) gauge symmetry, the
intrinsic origin of the criticality, then suggests absence of
“thermalization”. Extension of this work to an initial $\mathbb{Z}_{2}$
topological spin liquid and/or a longer-range (spin symmetric) interaction is
left for future studies.

_Acknowledgments—_ We acknowledge inspiring discussions with Ji-Yao Chen,
Olivier Gauthé, Norbert Schuch, Luca Tagliacozzo, Laurens Vanderstraeten,
Frank Verstraete and support from the TNTOP ANR-18-CE30-0026-01 grant awarded
by the French Research Council. This work was granted access to the HPC
resources of CALMIP center under the allocation 2017-P1231.

##  Appendix A Implementation of $C_{4v}$ and SU(2) symmetries

Here we provide details on the steps represented in Fig. 5 (h-i). Once a
uniform consistent gauge is fixed for the four virtual legs (see text), the
evolved tensor is projected onto the symmetric basis $\\{T_{a}\\}$. In this
process the resulting tensor becomes fully symmetric under $C_{4v}$ and SU(2),
in particular the noise $\varepsilon$ introduced during the gauge fixing
process is removed.

Figure 14: (a) Tensor increment $T_{i}$ and (b) Tensor dissymetry $T_{d}$
using the same definitions as in Fig. 5.

In order to test the accuracy of this procedure, let us define two $d\times d$
matrices according to Fig. 14:

  * •

The tensor increment $T_{i}$ evaluates the variation of site tensor under time
evolution. From its definition as a bilayer, it is obviously gauge
independent. Since ${\cal A}(t+\tau)-{\cal A}(t)=\alpha\tau+{\cal
O}(\tau^{2})$, the norm $||T_{i}||$ is expected to scale like $\tau$,

  * •

Tensor dissymetry $T_{d}$ measures the effect of explicit $C_{4v}$
symmetrization. It is evaluated using the gauge fixed tensor, as $C_{4v}$
symmetrization would be meaningless otherwise. Corrections induced by
symmetrization on ${\cal A}(t+\tau)$ are due to the non-commutativity of the 4
substeps and thus expected to occur at most at order $\tau^{2}$. Hence
$||T_{d}||$ should scale as $\tau^{2}$.

Figure 15: Small $\tau$ dependence of $||T_{i}||\sim 1.6287\tau$ (upper panel)
and $||T_{d}||\sim 0.817375\tau^{2}$ (lower panel).

We checked these scalings using ${\cal V}=0\oplus\frac{1}{2}\oplus 1$ (the
justification for using this ansatz was given in section III.2). Figure 15
shows that indeed $||T_{i}||\sim\tau$ and $||T_{d}||\sim\tau^{2}$. Hence it is
always possible to choose $\tau$ small enough such that
$||T_{d}||\ll||T_{i}||$. It is typically the case for $\tau=0.025$ used in the
following.

In a final stage (see Fig 5 (i)), the tensor is projected in the SU(2)
symmetric basis $\\{T_{a}\\}$. This basis being orthogonal and normalized, the
error of this projection is simply evaluated as $1-\sum_{a}|\langle{\cal
A}(t)|T_{a}\rangle|^{2}$ and we checked that it never exceeds $10^{-10}$ for
$\tau=0.025$.

##  Appendix B Factorisations of symmetric PEPS

In this section, we briefly describe the matrix factorizations that are
frequently used in the current work. In particular, we are interested in
factorizing a complex symmetric matrix into a diagonal form i.e.,

$A=VDV^{T}$ (10)

where $D$ is a diagonal matrix and $V$ is an invertible matrix. We also
explain how these matrices can be used to truncate the matrix $A$.

###  B.1 Autonne-Takagi Decomposition

The Autonne-Takagi factorization of a complex symmetric matrix $A$ is defined
as,

$A=USU^{T}$ (11)

Where $S$ is a diagonal matrix with positive real entries called the singular
values of the matrix A and U is a complex unitary. In order to find the
unitary $U$, we follow the procedure described in Ref. [33].

For any complex matrix, we define the singular value decomposition as,

$A=USV^{\dagger}$ (12)

where $S$ is the diagonal with singular values and $U$ and $V$ are unitary
matrices with corresponding singular vectors and $\dagger$ denotes Hermitian
conjugate.

We define the unitary matrix $Z=U^{\dagger}V^{*}$ where $V^{*}$ denotes the
complex conjugate of $V$. Since $A$ is complex symmetric, we can rewrite the
SVD as,

$\displaystyle USV^{\dagger}=$ $\displaystyle A=A^{T}=V^{*}SU^{T}$
$\displaystyle\Rightarrow$ $\displaystyle SV^{\dagger}U^{*}=U^{\dagger}V^{*}S$
$\displaystyle\Rightarrow$ $\displaystyle ZS=SZ^{T}$
$\displaystyle\Rightarrow$ $\displaystyle ZSZ^{*}=S$

Since the above equation preserves the spectrum of S, it is a similarity
transformation. Hence, $Z$ and $Z^{*}$ should be multiplicative inverse of
each other. i.e. $Z^{\dagger}=Z^{*}\Rightarrow\quad Z=Z^{T}$. This implies
that $Z$ and $S$ commute with each other.

$\displaystyle(ZS)_{ij}=$ $\displaystyle
z_{ij}S_{jj}\;\text{and}\;(SZ)_{ij}=S_{ii}z_{ij}$ $\displaystyle\Rightarrow$
$\displaystyle(z_{ij}S_{jj}-S_{ii}z_{ij})=0$ $\displaystyle\Rightarrow$
$\displaystyle z_{ij}(s_{i}-s_{j})=0$ $\displaystyle\Rightarrow$
$\displaystyle z_{ij}=0\quad\forall s_{i}\neq s_{j}$

Hence, we note that any such matrix $Z$ that commutes with a diagonal matrix
$S$ should be diagonal if all the entries in $S$ are distinct. If $S$ has
repeated singular values, $Z$ has a block diagonal structure with the size of
the blocks equal to the degeneracy of the singular values i.e.,
$Z=\bigoplus_{k}B_{n_{k}\times n_{k}}$ where $n_{k}$ is the multiplicity of
the $k$-th distinct singular value and $B_{n_{k}\times n_{k}}$ is the
corresponding block diagonal part of $Z$. Similarly, we can rewrite
$S=\bigoplus_{k}s_{k}\cdot{\mathbbm{1}}_{n_{k}\times n_{k}}$. Since every
block $B_{n_{k}\times n_{k}}$ commutes with the identity matrix, the overall
matrix too commutes with the diagonal matrix $S$. Since any power of $Z$ can
also be written in similar block diagonal form, they too commute with the
diagonal matrix $S$. Hence we can write, $Z^{\frac{1}{2}}S=SZ^{\frac{1}{2}}$.

$\displaystyle A$ $\displaystyle=$ $\displaystyle USV^{\dagger}$ (13)
$\displaystyle=$ $\displaystyle USV^{\dagger}U^{*}U^{T}$ $\displaystyle=$
$\displaystyle USZ^{T}U^{T}$ $\displaystyle=$ $\displaystyle
U(Z^{\frac{1}{2}})^{T}S(Z^{\frac{1}{2}})^{T}U^{T}$ $\displaystyle=$
$\displaystyle U_{Z}SU_{Z}^{T}$

We therefore get $A=U_{Z}SU_{Z}^{T}$ with $U_{Z}=U(Z^{\frac{1}{2}})^{T}$. Note
that the resultant unitary $U_{Z}$ is not unique, since the singular matrices
$U$ and $V$ themselves are not unique. Let two such unitaries $U_{Z_{1}}$ and
$U_{Z_{2}}$ be related by $U_{Z_{1}}=U_{Z_{2}}\phi$, where $\phi$ is a unitary
matrix, then,

$\displaystyle A$ $\displaystyle=$ $\displaystyle U_{Z_{1}}S(U_{Z_{1}})^{T}$
(14) $\displaystyle=$ $\displaystyle U_{Z_{2}}(\phi S\phi^{T})(U_{Z_{2}})^{T}$
$\displaystyle=$ $\displaystyle U_{Z_{2}}S(U_{Z_{2}})^{T}$

i.e. $(\phi S\phi^{T})=S$. Since this transformation preserves the spectrum of
S, it has to be a similarity transformation, i.e. $\phi\phi^{T}=I$. This also
implies that the matrix $\phi$ has to commute with $S$. Hence, if the singular
values $S$ are all distinct, the matrix $\phi$ is a orthogonal diagonal matrix
i.e. it’s a diagonal matrix with only $\pm 1$ as its diagonal entries. If $S$
has multiplicities, the matrix $\phi$ can take a block diagonal structure with
$n_{k}\times n_{k}$ orthogonal blocks where $n_{k}$ is the multiplicity of the
corresponding singular value $S_{k}$.

###  B.2 Orthogonal Decomposition

If the complex symmetric matrix $A$ is diagonalizable, one can diagonalize it
by using a set of eigenvectors i.e.,

$\displaystyle AE$ $\displaystyle=$ $\displaystyle ED$
$\displaystyle\Rightarrow A$ $\displaystyle=$ $\displaystyle EDE^{-1}$ (15)

where $D$ is a diagonal matrix whose entries are the eigenvalues of $A$ and E
is the matrix with the corresponding eigenvectors [34]. Now, we use the
symmetry argument $A=A^{T}$,

$\displaystyle A=$ $\displaystyle EDE^{-1}=A^{T}=(E^{-1})^{T}DE^{T}$
$\displaystyle\Rightarrow$ $\displaystyle ED=(E^{-1})^{T}DE^{T}E$
$\displaystyle\Rightarrow$ $\displaystyle ZD=DZ,$

where we have defined $Z=E^{T}E$. Note that eigenvectors of different
eigenvalues of any complex symmetric matrix are always orthogonal (not
orthonormal). Indeed, let $u$ and $v$ be two eigenvectors corresponding to two
different eigenvalues $\lambda$ and $\mu$, then,

$u^{T}\cdot Av-(u)^{T}A^{T}\cdot v=(\mu-\lambda)(u^{T}\cdot v)$ (16)

If the matrix $A$ is symmetric, the LHS of the equation should be zero, and
since $\lambda$ and $\mu$ are distinct, the dot product $u^{T}\cdot v$ should
be zero which can only happen when $u$ and $v$ are orthogonal. Hence the
matrix $Z=E^{T}E$ is a diagonal matrix if the A has nonrepeating distinct
eigenvalues. However, in the presence of multiplicities, the matrix exhibits a
block diagonal form. Unlike Takagi decomposition, this block diagonal form has
to do with the absolute values of the diagonal entries D. Eigenvalues which
have the same absolute value but are a complex phase away(like the complex
conjugate) can still have the same eigenvectors because if $u$ is an
eigenvector, the vector $e^{i\phi}u$ is also an eigenvector. Following the
same arguments as in the previous section, we get
$Z^{\frac{1}{2}}D=DZ^{\frac{1}{2}}$. Notice that the square root of Z should
also have the same block diagonal form as Z. Then, we get

$\displaystyle A$ $\displaystyle=$ $\displaystyle EDE^{-1}$ $\displaystyle=$
$\displaystyle EDE^{-1}(E^{T})^{-1}E^{T}$ $\displaystyle=$ $\displaystyle
EDZ^{-1}E^{T}$ $\displaystyle=$ $\displaystyle EZ^{-1/2}DZ^{-1/2}E^{T}$
$\displaystyle=$ $\displaystyle(EZ^{-1/2})D(EZ^{-1/2})^{T}$ $\displaystyle=$
$\displaystyle ODO^{T}\,.$

Thus the matrix $O=EZ^{-1/2}=E(E^{T}E)^{-1/2}$ is a set of orthogonal
(complex) eigenvectors of A. Note that in practice, eigenvectors of
eigenvalues which are distinct but whose absolute values are quite close can
generate very similar eigenvectors.

The major difference between the Autonne-Takagi decomposition and the
orthogonal decomposition is that, since the orthogonal decomposition results
in orthogonal projectors, the operation $O^{T}BO$ is a similarity
transformation on the matrix B whereas $U^{T}BU$ is not. Both methods can be
used to a factorize any symmetric matrix into a product of a matrix and its
transpose, i.e. $A=MM^{T}$ where $M=OD^{1/2}$ or $M=US^{1/2}$.

###  B.3 Relation between Singular values and Eigenvalues

Since both methods would give us a decomposition of the matrix into a product
of a matrix and it’s transpose, we can try to establish a relation between
them by comparing the resultant matrices. Since the number of non-zero
elements of $S$ and $D$ is a equal to the rank of the matrix, the size of non-
zero elements of $S$ and $D$ should be the same. Let us write
$A=USU^{T}=U\sqrt{S}\sqrt{S}U^{T}=(U\sqrt{S})(U\sqrt{S})^{T}$ or, similarly,
$A=ODO^{T}=O\sqrt{D}\sqrt{D}O^{T}=(O\sqrt{D})(O\sqrt{D})^{T}$. The
identification of the two forms,

$(U\sqrt{S})(U\sqrt{S})^{T}=(O\sqrt{D})(O\sqrt{D})^{T}\,,$

implies that

$(U\sqrt{S})\phi=(O\sqrt{D})$

where $\phi$ is a complex orthogonal matrix. Hence,

$O=U\sqrt{S}\phi\sqrt{D}^{-1}$

Making use of $O^{T}O=I$, we then obtain

$\displaystyle(U\sqrt{S}\phi\sqrt{D}^{-1})^{T}(U\sqrt{S}\phi\sqrt{D}^{-1})=I$
$\displaystyle\Rightarrow$ $\displaystyle\sqrt{S}U^{T}U\sqrt{S}=\phi
D\phi^{T}$ (18)

Comparing the diagonal elements in LHS and RHS of Eq. B.3, we get,

$\displaystyle C_{ii}$
$\displaystyle=\sum_{k}\sqrt{s}_{i}U_{ik}^{T}U_{ki}\sqrt{s_{i}}=\sum_{k}d_{k}\phi_{ik}\cdot\phi_{ki}^{T}$
(19)
$\displaystyle\Rightarrow\sum_{i=1}^{k}d_{i}=\sum_{i=1}^{k}s_{i}(U_{i}^{T}\cdot
U_{i})$

Hence, it is clear that the eigenvalues and singular values are equal if and
only if the unitary $U$ is purely real. Additionally, if the singular values
are degenerate, the pre-factor term $U^{T}_{i}U_{i}$ need not be the same.
Hence, $D$ might not always have the exactly same multiplet structure of $S$.
This is particularly important while truncating the matrix. Since we do not
want to cut through the multiplets, we should be more lenient with cutting
through multiplets of $D$, often grouping values that are slightly off as
degenerate. Additionally, we have,

$\sum_{i}|d_{i}|=\sum_{i}|s_{i}||(U_{i}^{T}U_{i})|\leq\sum_{i}s_{i}||U_{i}||\leq\sum_{i}s_{i}$

Thus, since sum total of magnitude of the first $n$ eigenvalues is lower than
the sum of first $n$ singular values, we can infer that the eigenvalue
spectrum decays quickly when compared to the singular values.

###  B.4 Truncation

Matrix truncation or dimensionality reduction is a transformation in which a
matrix is projected from a higher dimensional space to a lower dimensional
space such that the resultant matrix retains the maximum information and
properties. In other words, for a given square matrix $M_{p\times p}$ we find
a matrix that is lower order $N_{q\times q}$ for a given $q<p$ such that the
matrix $N_{q\times q}$ can replace $M_{p\times p}$ with as little change as
possible.
In order to truncate a matrix, we have to find an appropriate projector
$U_{p\times q}$ which projects the linear vector space of dimension $p$ to the
vector space of dimension $q$, i.e.,

$U^{T}M_{p\times p}V=N_{q\times q}$ (20)

Since we are dealing with complex symmetric matrices, we can assume without
the loss of generality that if the resultant matrix is a complex symmetric,
the projectors U and V are identical. As we project the matrix to a lower
dimensional space, there is a loss of information due to the truncation. We
quantify this by bringing the matrix $N_{q\times q}$ back to the vectors space
of $M_{p\times p}$ and take the norm of the difference of the resultant matrix
from the original. This reconstructed matrix $\tilde{U}N_{q\times q}\tilde{U}$
is called a Low Rank Approximation of the original matrix $M$ [35].

${\text{Loss}}=||M-\tilde{U}N\tilde{U^{T}}||$ (21)

where $\tilde{U}$ is the pseudo-inverse of matrix $U$. Note that from here on,
the word norm shall be used to refer the Frobenius norm of a matrix. Ideally,
in a tensor contraction, we apply $\tilde{U}U$ on every bond and by absorbing
them to the corresponding tensors, we reduce the dimension of the bond.

###  B.5 Loss in Autonne-Takagi truncation

In Autone-Takagi factorization, we use the complex conjugate of the unitary
$U^{\dagger}$ obtained by the factorization to diagonalize the matrix i.e.
$M=USU^{T}$. We then truncate the resultant matrix by keeping the $q$ of the
diagonal matrices and setting the rest to zero. Let us denote the truncated
diagonal matrix by $\tilde{S}(q)$.

In this case, the loss due to truncation is given by,

Loss $\displaystyle=$
$\displaystyle||M-U\tilde{S}(q)U^{T}||=||U(S-\tilde{S}(q))U^{T}||$ (22)
$\displaystyle=$
$\displaystyle||(S-\tilde{S}(q))||=\sqrt{\sum_{k=q+1}^{p}s_{k}^{2}}$

Hence, to minimize this loss, we retain the singular values that are largest
in magnitude. Once the appropriate values of $S$ which have to be deleted are
determined, we remove the corresponding columns in the matrix $U$ to get the
required projector.

###  B.6 Loss in Orthogonal truncation

In case of orthogonal decomposition, we use the orthogonal matrices obtained
by Eigenvalue decomposition in Eq. B.2. Since the resultant diagonal entries
can be complex, we truncate by retaining the $q$ values that have the largest
magnitude.

Estimation of loss due to truncation is tricky in case of complex matrices
that are not Hermitian. Even in complex symmetric case, we can never get a
exact rule estimating the loss. Hence, we try to find an upper bound for the
error. By proving that the upper bound reduces with increase in the final
dimension $q$, we comment that the quality of truncation, in general,
increases with the dimension.

Figure 16: Normalized eigenvalues (in log scale) of $720\times 720$ corner
matrices obtained from random PEPS ansätze given by Eq. 5 with bond dimension
$D=6$. Figure 17: Relative truncation error (in log scale) of $720\times 720$
corner matrices obtained from random PEPS ansätze given by Eq. 5 with bond
dimension $D=6$. The dotted lines indicate the upper bound calculated by Eq.
25 while the corresponding solid lines indicate the actual errors. Notice that
the general trend of error decreasing with decrease in truncation size.

In order to get an upper bound, we make use of the inequality that the
Frobenius norm of the product of two matrices is lesser than or equal to the
product of the norms of the individual matrices i.e., $||AB||\leq||A||||B||$
and the fact that a matrix and its transpose have the same entries and by
extension, the same Frobenius norm. From Eq. B.2, we Let us denote the
truncated diagonal matrix by $\tilde{D}(q)$ and the matrix with discarded
values by $\Delta(D)$. Thus the loss due to truncation is given by,

Loss $\displaystyle=$
$\displaystyle||M-O\tilde{D}(q)O^{T}||=||O\Delta(D)O^{T}||$ $\displaystyle=$
$\displaystyle||O\sqrt{\Delta(D)}(O\sqrt{\Delta(D)})^{T}||=||PP^{T}||\leq||P||^{2}$
$\displaystyle\Rightarrow{\text{Loss}}$ $\displaystyle\leq$
$\displaystyle||O\sqrt{\Delta(D)}||^{2}\leq||O||^{2}||{\Delta(D)}||.$ (23)

Now, using the fact that $\Delta(D)$ has only the $(p-q)$ smallest
eigenvalues, and the rest set to zero, the $O$ in the above equation can be
replaced by only the corresponding $(p-q)$ orthogonal eigenvectors i.e.,

$O\Delta\sqrt{D(q)}=O_{p\times(p-q)}\sqrt{\Delta{D(q)}_{(p-q)\times(p-q)}}$

Let us call this truncated $O$ with $(p-q)$ eigenvectors as $O_{reduced}$. Now
since, $O$ comprises of orthogonal vectors,

$\displaystyle p-q$ $\displaystyle=$ $\displaystyle||I_{(p-q)\times(p-q)}||$
$\displaystyle=$
$\displaystyle||O_{\text{reduced}}O^{T}_{\text{reduced}}||\leq||O_{\text{reduced}}||^{2}$
$\displaystyle\Rightarrow||O_{\text{reduced}}||$ $\displaystyle\geq$
$\displaystyle\sqrt{(p-q)}$ (24)

Hence, we get,

${\text{Loss}}\leq{(p-q)}||\Delta(D)||=(p-q)\sqrt{\sum_{k=q+1}^{p}|d_{k}|^{2}}$
(25)

For a given matrix of size $p$, the upper limit of the loss due to truncation
will go down with increase in the size of the truncated matrix $q$.

Notice that this upper limit doesn’t neccesarily imply that the error will
always keep decreasing with increase in $q$. as shown in Fig. 16 and Fig. 17
where we plot the eigenvalues of 25 random $729\times 729$ matrices obtained
from initial corners of CTMRG process. We can observe that the general trend
of truncation error goes down with the decrease in the magnitude of
truncation, though it is not always strictly decreasing.

###  B.7 Preferences

We shall now comment the usability of the above mentioned decompositions in a
given problem. It is clear from the truncation error that for any given case,
the singular value decomposition will have a much accurate low-rank
approximation. Hence, in any dimension reduction problems that involve a
single isolated matrix (or tensor), one has to always use singular value
decomposition. Examples include dimension reduction in the simple-update
procedure III.

Orthogonal decomposition is preferred in cases where one has to reduce the
bond of a tensor that is a part of closed loops or in cases where the overall
trace has to be minimized. In cases where the tensor contraction can be
reduced to the trace of product of two matrices, the trunction has to be done
by inserting an isometry on the legs of the bonds that are to be truncated.
Since we are dealing with complex symmetric matrices, for the resultant
isometries to preserve the complex symmetric nature of the matrices, the
projectors forming the isometry have to be complex symmetric, for which one
has to resort to the orthogonal decomposition. In this work, we use orthogonal
decomposition to renormalize the corner matric and the edge tensor in the
CTMRG procedure.

Figure 18: Single site symmetric environment obtained from the CTMRG
procedure.

##  Appendix C CTMRG

In this section, we provide a description of the CTMRG procedure used to
contract the iPEPS states used in this work (in order to compute observables
like energy, correlation functions, etc…). The details of the CTMRG are
similar to the one described in reference [30]. More specifically its
translationally invariant single-site $C_{4v}$-symmetric version [18]. In this
scheme, when we contract all the active bonds of the bi-layer tensors in the
thermodynamic limit, it leads to a $C_{4v}$-symmetric SU(2) environment with
adjustable bond dimension $\chi$ as shown in Fig. 18. The environment consists
of a corner matrix $C$ of dimension $\chi\times\chi$ and a rank-3 transfer
tensor $T$ of dimension $\chi\times\chi\times D^{2}$. However, the major
difference is that the corner in this work is no longer Hermitian but complex
symmetric. The edge tensor $T$ too is complex and exhibits reflection symmetry
along the axis through its $D^{2}$ leg. To reach the fixed point tensors,
first $C$ and $T$ are initialized by contracting the corresponding legs of ket
and bra tensor ${\cal A}^{\dagger}{\cal A}$. We then proceed to apply the
renormalization procedure which is continued till the environment converges to
a fixed point.

Figure 19: Single iteration of CTMRG. (a-b) It involves absorption of the
single site tensor $E$ into the $C$ (corner) and $T$ (edge) tensors increasing
the bond dimension from $\chi$ to $\chi D^{2}$ and leading to $C^{\prime}$ and
$E^{\prime}$ tensors. (c) Isometries (in red) obtained by orthogonal
decomposition of the enlarged corner are placed on the enlarged bonds. Their
absorption into the $T^{\prime}$ tensor defines the renormalized $\tilde{T}$
tensor, diagonalizes the corner $C^{\prime}$ to $\tilde{C}$ and reduces the
dimension back to $\chi$.

Each renormalization step consists of three parts.

  1. 1.

Absorption. At each step, the C and T tensors absorb a single site bi-layer
tensor as shown in Fig. 19. This raises the bond dimension of $C$ and $T$ to
$\chi D^{2}\times\chi D^{2}$ and $\chi D^{2}\times\chi D^{2}\times D^{2}$.

  2. 2.

Corner renormalization. Once absorbed, we can represent the enlarged
contraction of the bi-layer network of tensors by just contracting the four
identical enlarged corners. By leaving out one bond untraced, we can equate
the resultant matrix to the reduced density matrix (RDM) of the system. Since
this matrix is complex symmetric by construction, we can use the Orthogonal
decomposition as described in Appendix A to obtain the orthogonal projectors.
For stability, we just perform the orthogonal decomposition on the corner
since it has the same eigenvectors as that of the RDM. We diagonalize the
corner and retain the largest $\chi$ values. In order to preserve the SU(2)
symmetry of the environment, we make sure that the truncation is done between
values that are not part of a multiplet.

  3. 3.

Edge renormalization. We use the isometries obtained from the orthogonal
decomposition of the corner to renormalize the edge tensor $T$.
The procedure is repeated until convergence, which is identified by the
measuring the change in absolute values of the corner matrix.

## References

  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Reports on Progress in Physics 80, 016502 (2016).
  * Anderson [1973] P. W. Anderson, Resonating valence bonds: A new kind of insulator?, Materials Research Bulletin 8, 153 (1973).
  * Poilblanc  _et al._ [2012] D. Poilblanc, N. Schuch, D. Pérez-García, and J. I. Cirac, Topological and entanglement properties of resonating valence bond wave functions, Phys. Rev. B 86, 014404 (2012).
  * Chen and Poilblanc [2018] J.-Y. Chen and D. Poilblanc, Topological $\mathbb{Z}_{2}$ resonating-valence-bond spin liquid on the square lattice, Phys. Rev. B 97, 161107 (2018).
  * Rutkowski and Lawler [2016] T. C. Rutkowski and M. J. Lawler, Spin liquid phases of large-spin Mott insulating ultracold bosons, Phys. Rev. B 93, 094405 (2016).
  * Slagle  _et al._ [2022] K. Slagle, Y. Liu, D. Aasen, H. Pichler, R. S. K. Mong, X. Chen, M. Endres, and J. Alicea, Quantum spin liquids bootstrapped from Ising criticality in Rydberg arrays (2022).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021), https://www.science.org/doi/pdf/10.1126/science.abi8794 .
  * Giudici  _et al._ [2022] G. Giudici, M. D. Lukin, and H. Pichler, Dynamical preparation of quantum spin liquids in Rydberg atom arrays, Phys. Rev. Lett. 129, 090401 (2022).
  * Cheng  _et al._ [2021] Y. Cheng, C. Li, and H. Zhai, Variational approach to quantum spin liquid in a Rydberg atom simulator 10.48550/ARXIV.2112.13688 (2021).
  * Satzinger  _et al._ [2021] K. J. Satzinger, Y.-J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, C. Gidney, I. Aleiner, F. Arute, K. Arya, J. Atalaya, R. Babbush, J. C. Bardin, R. Barends, J. Basso, A. Bengtsson, A. Bilmes, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. R. Derk, D. Eppens, C. Erickson, L. Faoro, E. Farhi, A. G. Fowler, B. Foxen, M. Giustina, A. Greene, J. A. Gross, M. P. Harrigan, S. D. Harrington, J. Hilton, S. Hong, T. Huang, W. J. Huggins, L. B. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, T. Khattar, S. Kim, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, A. Locharla, E. Lucero, O. Martin, J. R. McClean, M. McEwen, K. C. Miao, M. Mohseni, S. Montazeri, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, T. E. O’Brien, A. Opremcak, B. Pató, A. Petukhov, N. C. Rubin, D. Sank, V. Shvarts, D. Strain, M. Szalay, B. Villalonga, T. C. White, Z. Yao, P. Yeh, J. Yoo, A. Zalcman, H. Neven, S. Boixo, A. Megrant, Y. Chen, J. Kelly, V. Smelyanskiy, A. Kitaev, M. Knap, F. Pollmann, and P. Roushan, Realizing topologically ordered states on a quantum processor, Science 374, 1237 (2021), https://www.science.org/doi/pdf/10.1126/science.abi8378 .
  * Alba and Calabrese [2017] V. Alba and P. Calabrese, Entanglement and thermodynamics after a quantum quench in integrable systems, Proceedings of the National Academy of Sciences 114, 7947 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1703516114 .
  * Robinson  _et al._ [2021] N. J. Robinson, A. J. J. M. de Klerk, and J.-S. Caux, On computing non-equilibrium dynamics following a quench, SciPost Phys. 11, 104 (2021).
  * Das [2020] S. R. Das,  _Quantum Quench and Universal Scaling_ (Oxford University Press, 2020).
  * Mitra [2018] A. Mitra, Quantum quench dynamics, Annual Review of Condensed Matter Physics 9, 245 (2018), https://doi.org/10.1146/annurev-conmatphys-031016-025451 .
  * Guardado-Sanchez  _et al._ [2018] E. Guardado-Sanchez, P. T. Brown, D. Mitra, T. Devakul, D. A. Huse, P. Schauß, and W. S. Bakr, Probing the quench dynamics of antiferromagnetic correlations in a 2d quantum Ising spin system, Phys. Rev. X 8, 021069 (2018).
  * Moessner and Sondhi [2003] R. Moessner and S. L. Sondhi, Ising and dimer models in two and three dimensions, Phys. Rev. B 68, 054405 (2003).
  * Moessner and Raman [2011] R. Moessner and K. S. Raman, Quantum dimer models, in  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 437–479.
  * Mambrini  _et al._ [2016] M. Mambrini, R. Orús, and D. Poilblanc, Systematic construction of spin liquids on the square lattice from tensor networks with SU(2) symmetry, Phys. Rev. B 94, 205124 (2016).
  * Penrose [1971] R. Penrose, Applications of negative dimensional tensors, in  _Combinatorial Mathematics and its Applications,_ (Academic Press, 1971).
  * Jordan  _et al._ [2008] J. Jordan, R. Orús, G. Vidal, V. F., and J. Cirac, Classical simulation of infinite-size quantum lattice systems in two spatial dimensions, Phys. Rev. Letters 101, 250602 (2008).
  * Schuch  _et al._ [2012] N. Schuch, D. Poilblanc, J. I. Cirac, and D. Pérez-García, Resonating valence bond states in the PEPS formalism, Physical Review B 86, 115108 (2012).
  * Dreyer  _et al._ [2020] H. Dreyer, L. Vanderstraeten, J.-Y. Chen, R. Verresen, and N. Schuch, Robustness of critical U(1) spin liquids and emergent symmetries in tensor networks, arXiv preprint arXiv:2008.04833 10.48550/ARXIV.2008.04833 (2020).
  * Jiang  _et al._ [2008] H. C. Jiang, Z. Y. Weng, and T. Xiang, Accurate determination of tensor network state of quantum lattice models in two dimensions, Phys. Rev. Lett. 101, 090603 (2008).
  * Autonne [1915] L. Autonne, Sur les matrices hypohermitiennes et sur les matrices unitaires, Ann. Univ. Lyon 38, 1 (1915).
  * Takagi [1924] T. Takagi, On an algebraic problem related to an analytic theorem of Carathéodory and Fejér and on an allied theorem of Landau, Jpn. J. Math. 1, 83 (1924).
  * Suzuki [1990] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations, Physics Letters A 146, 319 (1990).
  * Nishino and Okunishi [1996] T. Nishino and K. Okunishi, Corner transfer matrix renormalization group method, Journal of the Physical Society of Japan 65, 891 (1996).
  * Nishino  _et al._ [1996] T. Nishino, K. Okunishi, and M. Kikuchi, Numerical renormalization group at criticality, Physics Letters A 213, 69 (1996).
  * Orús and Vidal [2009] R. Orús and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Physical Review B 80, 094403 (2009).
  * Orús [2012] R. Orús, Exploring corner transfer matrices and corner tensors for the classical simulation of quantum lattice systems, Physical Review B 85, 205117 (2012).
  * Calabrese and Cardy [2004] P. Calabrese and J. Cardy, Entanglement entropy and quantum field theory, Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
  * Calabrese  _et al._ [2009] P. Calabrese, J. Cardy, and B. Doyon, Entanglement entropy in extended quantum systems, Journal of Physics A: Mathematical and Theoretical 42, 500301 (2009).
  * Chebotarev and Teretenkov [2014] A. M. Chebotarev and A. E. Teretenkov, Singular value decomposition for the Takagi factorization of symmetric matrices, Applied Mathematics and Computation 234, 380 (2014).
  * Strang [2009] G. Strang,  _Introduction to Linear Algebra_ , 4th ed. (Wellesley-Cambridge Press, Wellesley, MA, 2009).
  * Eckart and Young [1936] C. Eckart and G. Young, The approximation of one matrix by another of lower rank, Psychometrika 1, 211 (1936).