# Stochastic Sampling of Operator Growth Dynamics

Ayush De  Racah Institute of Physics, The Hebrew University of Jerusalem,
Jerusalem 91904, Israel  Umberto Borla  Racah Institute of Physics, The Hebrew
University of Jerusalem, Jerusalem 91904, Israel  Xiangyu Cao  Laboratoire de
Physique de l’École normale supérieure, ENS, Université PSL, CNRS, Sorbonne
Université, Université Paris Cité, F-75005 Paris, France  Snir Gazit  Racah
Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904,
Israel  The Fritz Haber Research Center for Molecular Dynamics, The Hebrew
University of Jerusalem, Jerusalem 91904, Israel

###### Abstract

We put forward a Monte Carlo algorithm which samples the Euclidean time
operator growth dynamics at infinite temperature. Crucially, our approach is
free of the numerical sign problem for a range of quantum many-body spin
systems, allowing for numerically exact and unbiased calculations. We apply
this methodological headway to study the high-frequency dynamics of the mixed-
field quantum Ising model (QIM) at one and two dimensions. Our results
corroborate the recently proposed operator growth hypothesis, predicting an
exponential fall-off of generic response functions at large frequencies.
Remarkably, our calculations are sufficiently sensitive to detect subtle
logarithmic corrections of the hypothesis in one dimension. In addition, in
two dimensions, we uncover a non-trivial crossover between two large frequency
decay rates. Lastly, we find that the growth rate of the operator support size
is dominated by boundary effects, with a logarithmic correction in one
dimension.

Introduction  — Within the Heisenberg representation of quantum mechanics, the
time evolution of operators is governed by their commutation relations with
the Hamiltonian. Under generic many-body dynamics, local operators are
expected to evolve into increasingly complex and non-local operator strings.
This process of “operator growth” is equivalent to the scrambling of quantum
information in complex systems, and closely related to quantum chaos. As such,
understanding the precise properties of operator growth is a fundamental line
of research that has driven a great amount of theoretical and experimental
work [1, 2, 3, 4, 5, 6, 7, 8, 9]. It has direct links to the question of how
the laws of quantum mechanics bring about thermalization or its failure in
isolated systems [10, 11, 12, 13, 14]. It is also a useful tool for
investigating the fate of information in a black hole and how to retrieve it
[15, 16, 17, 18, 19].

In this context, an interesting way to quantify the rate of operator growth,
besides the much studied out-of-time order correlator, is to examine the high-
frequency regime of the spectral density associated with generic correlators
[20]. Recently, Ref. [21] proposed an operator growth hypothesis that
translates to an exponentially decaying large frequency tail

$\Phi(\omega)\sim\exp(-|\omega|/\omega_{0})$ (1)

in generic many-body systems with finite local bandwidth. The energy scale
$\omega_{0}$ provides an upper bound to the Lyapunov exponent of out-of-time
order correlators [21, 22, 23], and is thus an alternative measure of operator
growth [24, 25, 26, 27, 28, 29, 30]. Crucially, this quantity is
experimentally accessible [31], e.g., from the pre-thermalization time in a
fast periodic drive [32].

Testing the operator growth hypothesis in concrete lattice models is
particularly challenging due to the inherent difficulty of simulating many-
body quantum dynamics with classical resources. Direct evidence for the
hypothesis is largely limited to large-N or semi-classical models [23, 30, 33,
20, 34] and brute-force numerics of small to intermediate system sizes that
suffer from sizable finite size effects [35]. Although analytical bounds
supporting the conjecture have been obtained for finite-dimensional systems
[36, 37], these results are not satisfactory in that they do not provide an
accurate estimate of $\omega_{0}$ or the precise functional form of the
spectral function. It is, therefore, desirable to devise efficient numerical
schemes that can reliably test the operator growth hypothesis.

(a)  (b)

Figure 1:  (a) A space-time snapshot of the Euclidean time operator growth
process extracted from the $400^{\text{th}}$ moment of the Pauli $Z$ operator
at site 0 of the 1D mixed field QIM. Transformations are generated by
Liouvillians between time slices (indexed by $k$). The colors correspond to
the different Pauli operators. (b) The average occupation probability on a
site by a non-identity Pauli operator resulting from operator growth arising
from the same Monte Carlo dynamics as (a).

In this Letter, we present a Monte Carlo method that samples the Euclidean
time operator growth dynamics at infinite temperature in an unbiased and
statistically exact manner. The Monte Carlo method tracks the Euclidean time
evolution of operators (Fig. 1a) instead of the standard wave function
approach that is common in conventional quantum Monte Carlo calculations.
Crucially, our scheme does not suffer from the notorious numerical sign
problem when applied to a broad family of lattice spin models of interest. We
employ our method to study the mixed-field QIM in both one and two dimensions.
Our results provide high-precision tests of the hypothesis. In particular, we
pinpoint the subtle log-correction to Eq. 1 in one dimension, and unveil an
intriguing crossover between two high-frequency tails in two dimensions.
Finally, we determine the dynamical scaling properties of the operator support
size.

Operator-space partition function — In statistical physics, Monte Carlo
schemes aim to sample from a partition function representing a weighted sum
over the configuration space. In our case, the relevant partition function is
the auto-correlation function of a local operator, $O$, analytically continued
to imaginary time:

$\mathcal{Z}(\tau):=(O|e^{\mathcal{L}\tau}O)\,.$ (2)

Here, $\mathcal{L}O=[H,O]$ is the Liouvillian super operator representing the
commutation with the Hamiltonian $H$, and
$(A|B)=\mathrm{Tr}[A^{\dagger}B]/\mathrm{Tr}[1]$ denotes the infinite-
temperature scalar product on operator space. We shall also assume that $O$
equals a product of one-site operators. To make contact with the real
frequency spectral density (and especially its large frequency tails), it is
useful to consider the Taylor expansion of $Z(\tau)$ at $\tau=0$:

$\displaystyle\mathcal{Z}(\tau)=\sum_{n=0}^{\infty}\mu_{2n}\frac{\tau^{2n}}{(2n)!}$
(3)
$\displaystyle\mu_{2n}=(O|\mathcal{L}^{2n}|O)=\int\Phi(\omega)\omega^{2n}\frac{\mathrm{d}\omega}{2\pi}.$
(4)

In the last equation, we identify the Taylor expansion coefficients
$\\{\mu_{2n}\\}$ with the moments of the spectral function $\Phi(\omega)$,
which follows from a Fourier transform of the analytic continuation to real-
time correlators, $\mathcal{C}(t)=\mathcal{Z}(\tau=it)$. Thus, high moments
encode information about the high-frequency tail of $\Phi(\omega)$. In
particular, an exponential tail (Eq. 1) yields:

$\mu_{2n}\simeq\left(\frac{2n\omega_{0}}{e}\right)^{2n}$ (5)

for large $n$. We recall that the moments grow rapidly, such that
$\mathcal{Z}(\tau)$ converges only when $\tau<\tau_{*}=1/\omega_{0}$, and has
a pole at $\tau_{*}$.

To evaluate the moments via stochastic sampling, we express $\mu_{2n}$ as a
sum over paths in operator space. Assuming local dynamics, we write
$\mathcal{L}=\sum_{i,a}\mathcal{L}_{i,a}$ , with
$\left\\{\mathcal{L}_{i,a}\right\\}$ denoting local Liouvillian of type $a$
acting on site $i$ , as dictated by the underlying Hamiltonian,
$H=\sum_{i,a}H_{i,a}$, with $\mathcal{L}_{i,a}O=[H_{i,a},O]$. Following the
above definitions, we expand the product $(O|\mathcal{L}^{2n}|O)$, by
introducing operator space resolution of identities between individual
Liouvillians:

$\mu_{2n}=\sum_{\mathcal{S}_{2n}}\underbrace{\prod_{k=0}^{2n-1}(O_{k+1}|\mathcal{L}_{i_{k},a_{k}}O_{k})}_{W_{\mathcal{S}_{2n}}}\,,$
(6)

where the sum is over all local Liouvillian sequences, specified by $2n$
position-type pairs, $\mathcal{S}_{2n}=\\{(i_{k},a_{k})\\}_{k=0}^{2n-1}$, with
nonzero weight $W_{\mathcal{S}_{2n}}\neq 0$. This construction implicitly
defines operator histories $(O_{0}=O,O_{1},\dots,O_{2n}=O_{0})$ and is
reminiscent of the Stochastic Series Expansion (SSE) [38, 39]. We emphasize
that, unlike standard quantum Monte Carlo approaches, here, the configuration
space is defined via the space of local operators, and transitions between
them are generated by the Liouvillian.

The weights $W_{\mathcal{S}}$ are generally non-positive and can lead to a
sign problem, even in cases where the associated Hamiltonian problem is sign
problem free. Remarkably, for a class of nontrivial models, the sum in Eq. 6
contains strictly non-negative real terms, $W_{\mathcal{S}}\geq 0$ [37]. This
class includes the mixed-field QIM, with the Hamiltonian,

$H=-J\sum_{\left\langle{i,j}\right\rangle}Z_{i}Z_{j}-\sum_{i}(h^{Z}Z_{i}+h^{X}X_{i})\,.$
(7)

Here, $\left\langle{i,j}\right\rangle$ define nearest neighbours on a
hypercubic lattice, and $X_{i},Z_{i}$ are the Pauli operators act on sites $i$
along the $x,z$ spin axes, respectively. In what follows, we will focus on the
numerical study of the mixed field QIM with positive and uniform $J$, $h^{Z}$,
$h^{X}$.

Figure 2:  A typical configuration contributing to the $4$-th moment with the
horizontal and vertical axes being temporal and spatial respectively. The
sites are populated by the $X$, $Y$, and $Z$ Pauli operators or the $I$
(identity) operator. The temporal direction is generated by the Liouvillians
of type $\mathcal{L}_{Z}$, $\mathcal{L}_{ZZ}$ and $\mathcal{L}_{X}$
corresponding to the three types of local Hamiltonian terms. The empty
Liouvillian slots correspond to null superoperators.

Numerical methods and observables — Following the expansion of
$\mathcal{Z}(\tau)$ in Eq. 3, we now turn to devise a Monte Carlo algorithm to
numerically estimate the moments $\mu_{2n}$. To that end, we set $\tau=1/J$ in
Eq. 2 and restrict the expansion order to obey $2n\leq 2n_{\max}$. Analogously
to SSE, we consider Liouvillian sequences, $\mathcal{S}_{2n_{\max}}$, of fixed
length $2n_{\max}$, which for an expansion order $n$ include additional
$(2n_{\max}-2n)$ factors of the null super-operator $\mathcal{L}_{0,0}\equiv
I$.

We note that the above truncation is natural in the context of imaginary time
Hamiltonian dynamics via SSE simulations. There, for given values of
microscopic parameters and temperature, the probability distribution of the
expansion order rapidly decays away from its mean value. By contrast, here, we
will be interested in targeting _all_ moments $\mu_{2n}$ for $2n\leq
2n_{\max}$. As will become clearer below, this entails visiting each expansion
order at roughly equal frequency.

The resulting partition function has the following form,

$\displaystyle G=$
$\displaystyle\sum_{\mathcal{S}_{2n_{\max}}}\frac{1}{J^{2n}(2n)!}$ (8)
$\displaystyle\times\frac{(2n_{\max}-2n)!(2n)!}{(2n_{\max})!}\prod_{k=0}^{2n_{\max}-1}(O_{k+1}|\mathcal{L}_{i_{k},a_{k}}|O_{k}),$

where the combinatorial factors arise from the redundancy in the assignment of
null operators. Eq. 8 has a natural diagrammatic interpretation as a sum over
$D=d+1$ dimensional space-time configurations consisting of on-site Pauli
operator, with the dynamics generated by insertions of local Liouvillian
super-operators, transitioning between Pauli strings, see Fig. 2.

In the following, we present a set of local Monte Carlo moves that enable
traversing between all allowed configurations: (1) add-drop – allows switching
between different $n$ sectors by adding or removing a pair of identical local
Liouvillians located on the same sites at two consecutive time slices. (2)
swap – swaps two commuting Liouvillians at different time slices, without
changing the expansion order $n$. This move is proposed only if there are no
non-commuting Liouvillians at intermediate time slices. (3) interact – here,
two temporally consecutive $\mathcal{L}_{ZZ}$ Liouvillians that share a single
spatial site “interact” to create a pair of temporally consecutive
$\mathcal{L}_{Z}$ Liouvillians acting on next to nearest neighbour sites. This
move allows accessing configurations with an odd number of $\mathcal{L}_{ZZ}$
and $\mathcal{L}_{Z}$ Liouvillians acting on a given site.

In practice, we found it essential to apply a reweighting approach, in the
spirit of the Wang-Landau algorithm [40, 41, 42], in order to ensure
ergodicity due to rapid fall off (or rise) of the expansion order distribution
towards large moments. This is achieved by iteratively tuning the relative
visiting rate of each expansion order until reaching an approximately uniform
distribution across expansion order sectors. In other words, we sample the
partition sum $\mathcal{Z}(\tau)$ near criticality, $\tau\nearrow\tau_{*}$.
Additional technical details on the Monte Carlo algorithm are relegated to
[43].

Our key physical observable is the moments’ ratio:

$r_{n}^{2}:=\frac{\mu_{2n+2}}{\mu_{2n}}=\frac{\left\langle\delta_{n^{\prime},n+1}\right\rangle_{n^{\prime}}}{\left\langle\delta_{n^{\prime},n}\right\rangle_{n^{\prime}}}\times
J^{2}(2n+2)(2n+1).$ (9)

In the above equation, averaging is carried out over the probability
distribution of the expansion order $P(n^{\prime})$, where the Kronecker delta
$\delta_{n^{\prime},n}$ counts Monte Carlo events satisfying $n^{\prime}=n$.
Effectively, $r^{2}_{n}$ is tied to the ratio between the visiting frequencies
of expansion orders $n+1$ and $n$. The large frequency behavior of Eq. 5,
predicts a linear growth of $r_{n}$ with $n$,

$r_{n}\simeq\frac{2\omega_{0}n}{e}.$ (10)

In addition, the prefactor allows a numerical estimation of the decay rate
$\omega_{0}$.

To estimate the spatial extent of operator growth, we define the support size
at time slice $k$ as:

$S_{k}=\langle\text{Volume}(B^{d}_{k})\rangle.$ (11)

Here, $B^{d}_{k}$ is geometrically defined as the largest region in $d$
dimensional space enclosed by non-trivial operators at temporal index $k$.

(a)  (b)  (c)  (d)

Figure 3:  (a) $r_{n}$ vs $n$ for the $Z$ operator with $J=h^{Z}=h^{X}=1$ in
1D. The inset zooms in on the first 35 moments, comparing our calculations
with exact brute force numerics. (b) $\Delta r_{n}$ vs $n$ for 1D. This
asymptotically falls to zero. Multiplying by $\log(n)$ recovers a constant
plateau. (c) Same as (a) for 2D. The inset has the first 16 moments compared
against exact computation. (d) $\Delta r_{n}$ vs $n$ for 2D. Unlike 1D here
the $\log$ multiplication isn’t needed. The inset shows $n\Delta r_{n}$ vs
$n$.

Numerical Results –  In the following, we present Monte Carlo data for the $Z$
operator dynamics in the mixed field QIM with the microscopic parameters
$h^{X}=h^{Z}=J=1$ in 1D and 2D. We begin our presentation by examining the 1D
case. We recall that, due to the restricted dynamics in 1D, the linear growth
in Eq. 10 receives a logarithmic correction such that $r_{n}\sim n/\log{n}$
[44, 36, 32, 45, 21]. In Fig. 3a, we depict $r_{n}$ as a function of $n$. We
first note the excellent agreement with brute force calculation (suffering
from an unfavourable exponential scaling with the expansion order) as shown in
the inset, serving as a non-trivial benchmark of our approach.

Furthermore, we observe that the large $n$ growth of $r_{n}$ deviates from a
linear trend. To pin down the source of this deviation and check whether it
can be attributed to the aforementioned logarithmic correction in Fig. 3b, we
plot $\Delta r_{n}\times\log(n)$, where $\Delta r_{n}$ is a numerical
derivative extracted using the five-point stencil method. This quantity is
expected to reach a plateau at large $n$ with the logarithmic correction.
Indeed, we find that only in the presence of the log multiplication, the
desired plateau is obtained. By contrast, $\Delta r_{n}$ asymptotically decays
to zero due to the $1/\log(n)$ factor. We note that this is the first
numerically exact confirmation for the predicted log correction associated
with the restricted 1D dynamics, as previous brute-force numerics cannot
resolve this subtle functional dependence convincingly due to the limited
number of accessible moments [21, 35, 46].

We now turn to study 2D systems. In Fig. 3c, we plot the evolution of $r_{n}$
with the expansion order $n$. As before, we note the precise agreement with
brute force computation, shown in the inset. Impressively, the stochastic
approach allows access to the first 180 moments, as compared to brute force
methods that are limited to the first 16 moments. We observe a clear linear
trend of $r_{n}$, which is the first extensive test of the operator growth
hypothesis in a concrete model beyond 1D systems.

In Fig. 3d, we plot $\Delta r_{n}$, and remarkably, we find an intriguing two-
plateau structure. We infer that the spectral density features a crossover
between two exponential tails. Indeed, if we suppose

$\Phi(\omega)\sim\begin{cases}e^{-|\omega|/\omega_{1}}&|\omega|\lesssim\omega_{c}\\\
e^{-|\omega|/\omega_{2}}&|\omega|\gtrsim\omega_{c}\end{cases}$ (12)

where $\omega_{1}$ and $\omega_{2}$ are the decay rates and $\omega_{c}$ is
the crossover scale. Using a saddle point approximation of Eq. 4, we may find
that

$\Delta r_{n}\approx\begin{cases}2\omega_{1}/e&n\lesssim
n_{1}=\omega_{c}/2\omega_{1}\\\ 2\omega_{2}/e&n\gtrsim
n_{2}=\omega_{c}/2\omega_{2}\end{cases}$ (13)

In particular, we have $n_{1}\omega_{1}=n_{2}\omega_{2}=\omega_{c}/2$. To show
our numerical results are consistent with Eq. 12, we plot $n\Delta r_{n}$ in
the inset of Fig. 3d, and find a remarkably large $\omega_{c}\approx 90J$.
This crossover would have been missed in a brute-force calculation.

Owing to the thermal nature of the dynamics, we expect the large $n$ behavior
of $r_{n}$ to be universal for typical local operators with a given set of
model parameters. In [43], we indeed see this anticipated behavior for
boundary $X$ operator, with $\log(n)\Delta r_{n}\simeq 2.2(1)J$ and $\Delta
r_{n}\simeq 1.26(4)J$ in 1D and 2D respectively, consistent with the plots
here. We also find that the structure of the intermediate plateau as seen in
2D is not universal. Nevertheless, the asymptotic plateau value is reached
only at a similar large value $n_{2}\approx 40$. Correspondingly, the high-
frequency tail of the spectral density is truly established only at a
frequency that was numerically inaccessible before the present work.

_Statistical mechanics perspective of operator growth_ – From the point of
view of statistical mechanics, it is natural to examine the geometric
properties of the operator string configurations Fig. 1b. In particular, the
configurations with large $n_{\max}$ describe the vicinity of the critical
point $\tau\nearrow\tau_{*}$ in the partition function (3) (in 1D
$\tau_{*}=\infty$ because of the log-correction). Therefore, these large
$n_{\max}$ configurations exhibit universal, scale-invariant behavior
ubiquitous to continuous phase transitions.

Figure 4:  Scaling collapse of the average operator support size (for the $Z$
operator), as function of $k$, for different values of $n_{\max}$, in 1D (a),
where $S_{\max}:=n_{\max}/\log n_{\max}$, and in 2D (b).

To illustrate this point, we measured the average support size of the operator
string as a function of temporal index $k$ in configurations of an increasing
range of maximal expansion order $n_{\max}$. The results, plotted in Fig. 4,
are in agreement with the following scaling ansatz in the large $n_{\max}$
limit:

$\displaystyle\left<S_{\text{2D}}\right>=n_{\max}f_{\text{2D}}(k/2n_{\max}),$
(14a)
$\displaystyle\left<S_{\text{1D}}\right>=\frac{n_{\max}}{\log(n_{\max})}f_{\text{1D}}(k/2n_{\max})$
(14b)

where $f_{\text{2D}}$ and $f_{\text{1D}}$ are scaling functions, which capture
the evolution of the operator size growth to an $n_{\max}$-dependent maximum
(and then decrease) in terms of rescaled time $k/2n_{\max}$ and size
$S/S_{\max}$.

The maximal size $S_{\max}\sim n_{\max}$ in 2D, and we expect the same in
higher dimensions as well, since this is the only asymptotics compatible with
the universal moment growth in Eq. 5. The slow convergence to the scaling
collapse echos the late establishment of the universal growth observed above.
In 1D, the maximal size $S_{\max}=n_{\max}/\log(n_{\max})$ has a log-
correction, which has the same geometric origin as the log-correction to the
moment growth: in 1D, the operator can only grow at the two extremities, which
is a severe entropy penalty. Therefore the dominant contributions must be a
compromise: the operator grows less in order to “scramble” more in the bulk.
Previous arguments [44, 36] indicated a three-stage process where the operator
grows to $\sim S_{\max}$ at $k\sim S_{\max}$, then stays as such until
$2n_{\max}-k\sim S_{\max}$ at which point the operator starts to shrink. In
units of rescaled time $k/2n_{\max}$, the growth and shrinking stages’
duration would be vanishing $\sim 1/\log(n_{\max})$, and the scaling function
$f_{\text{1D}}$ would be a constant. In Fig. 4 we see indeed that a plateau
establishes itself and extends to the whole interval as $n_{\max}$ increases,
however very slowly (due to the log dependence). In 2D, one may worry that a
similar “boundary-versus-bulk” entropy penalty applies as well. However,
delicate rigorous arguments confirm that this is not the case [36]. Put
simply, these arguments state that the operator support in 2D can have a
“rough” boundary (with length comparable to the bulk area) and thus enough
ways to grow in order to generate the moment asymptotics (5). Our results on
moment ratio asymptotics provide strong (yet indirect) evidence of this
unexpected geometric fact, whose direct observation will require pushing the
numerics to large $n$.

_Discussion_ – We conclude our presentation by flagging several research
directions motivated by our results. Our approach can be extended to models
whose associated Liouvillian dynamics are expressed as sign problem-free sums.
The quantum XY model is a notable example [37], which crucially admits a
global $\text{U}(1)$ symmetry and hence can provide access to transport
properties of conserved charges [47, 48]. The presence of disorder is also
compatible with our method, it will be interesting to investigate the
signature of Many-Body Localization in the evolution of operator strings.
Another important direction is extracting real frequency dynamical spectral
functions away from the high-frequency limit, as can, in principle, be via
Lanczos coefficients expressed in terms of moment ratios [49, 50, 51]. We note
that the latter might suffer from high sensitivity to statistical noise and
may require an improvement on current algorithms. Lastly, we expect the
support size scaling functions to be universal and valid for a generic class
of chaotic Hamiltonians, and predicting them theoretically is an interesting
open question. We leave the aforementioned possibilities to future efforts.

_Acknowledgements_ – We thank Ehud Altman and Assa Auerbach for fruitful
discussions. S.G. acknowledges support from the Israel Science Foundation
(ISF) Grant no. 586/22 and the US–Israel Binational Science Foundation (BSF)
Grant no. 2020. U.B. acknowledges support from the Israel Academy of Sciences
and Humanities through the Excellence Fellowship for International
Postdoctoral Researchers. Computational resources were provided by the Intel
Labs Academic Compute Environment and the Fritz Haber Center for Molecular
Dynamics, Hebrew University.

## References

  * Swingle  _et al._ [2016] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, Measuring the scrambling of quantum information, Phys. Rev. A 94, 040302 (2016).
  * Nahum  _et al._ [2018] A. Nahum, S. Vijay, and J. Haah, Operator spreading in random unitary circuits, Phys. Rev. X 8, 021014 (2018).
  * Khemani  _et al._ [2018a] V. Khemani, A. Vishwanath, and D. A. Huse, Operator spreading and the emergence of dissipative hydrodynamics under unitary evolution with conservation laws, Phys. Rev. X 8, 031057 (2018a).
  * Khemani  _et al._ [2018b] V. Khemani, D. A. Huse, and A. Nahum, Velocity-dependent lyapunov exponents in many-body quantum, semiclassical, and classical chaos, Phys. Rev. B 98, 144304 (2018b).
  * Gopalakrishnan  _et al._ [2018] S. Gopalakrishnan, D. A. Huse, V. Khemani, and R. Vasseur, Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems, Phys. Rev. B 98, 220303 (2018).
  * Qi and Streicher [2019] X.-L. Qi and A. Streicher, Quantum epidemiology: operator growth, thermal effects, and syk, Journal of High Energy Physics 2019, 12 (2019).
  * Roberts  _et al._ [2018] D. A. Roberts, D. Stanford, and A. Streicher, Operator growth in the syk model, Journal of High Energy Physics 2018, 122 (2018).
  * Xu and Swingle [2020] S. Xu and B. Swingle, Accessing scrambling using matrix product operators, Nature Physics 16, 199 (2020).
  * Swingle [2018] B. Swingle, Unscrambling the physics of out-of-time-order correlators, Nature Physics 14, 988 (2018).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Sekino and Susskind [2008] Y. Sekino and L. Susskind, Fast scramblers, Journal of High Energy Physics 2008, 065 (2008).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, Journal of High Energy Physics 2014, 67 (2014).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 2016, 106 (2016).
  * Brown  _et al._ [2023] A. R. Brown, H. Gharibyan, S. Leichenauer, H. W. Lin, S. Nezami, G. Salton, L. Susskind, B. Swingle, and M. Walter, Quantum gravity in the lab. i. teleportation by size and traversable wormholes, PRX Quantum 4, 010320 (2023).
  * Schuster  _et al._ [2022] T. Schuster, B. Kobrin, P. Gao, I. Cong, E. T. Khabiboulline, N. M. Linke, M. D. Lukin, C. Monroe, B. Yoshida, and N. Y. Yao, Many-body quantum teleportation via operator spreading in the traversable wormhole protocol, Phys. Rev. X 12, 031013 (2022).
  * Elsayed  _et al._ [2014] T. A. Elsayed, B. Hess, and B. V. Fine, Signatures of chaos in time series generated by many-spin systems at high temperatures, Phys. Rev. E 90, 022910 (2014).
  * Parker  _et al._ [2019] D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, A universal operator growth hypothesis, Phys. Rev. X 9, 041017 (2019).
  * Murthy and Srednicki [2019] C. Murthy and M. Srednicki, Bounds on chaos from the eigenstate thermalization hypothesis, Phys. Rev. Lett. 123, 230606 (2019).
  * Gu  _et al._ [2022] Y. Gu, A. Kitaev, and P. Zhang, A two-way approach to out-of-time-order correlators, Journal of High Energy Physics 2022, 133 (2022).
  * Barbón  _et al._ [2019] J. L. F. Barbón, E. Rabinovici, R. Shir, and R. Sinha, On the evolution of operator complexity beyond scrambling, Journal of High Energy Physics 2019, 264 (2019).
  * Rabinovici  _et al._ [2021] E. Rabinovici, A. Sánchez-Garrido, R. Shir, and J. Sonner, Operator complexity: a journey to the edge of krylov space, Journal of High Energy Physics 2021, 62 (2021).
  * Rabinovici  _et al._ [2022] E. Rabinovici, A. Sánchez-Garrido, R. Shir, and J. Sonner, Krylov complexity from integrability to chaos, Journal of High Energy Physics 2022, 151 (2022).
  * Avdoshkin and Dymarsky [2020] A. Avdoshkin and A. Dymarsky, Euclidean operator growth and quantum chaos, Phys. Rev. Res. 2, 043234 (2020).
  * Dymarsky and Smolkin [2021] A. Dymarsky and M. Smolkin, Krylov complexity in conformal field theory, Phys. Rev. D 104, L081702 (2021).
  * Avdoshkin  _et al._ [2022] A. Avdoshkin, A. Dymarsky, and M. Smolkin, Krylov complexity in quantum field theory, and beyond, arXiv:2212.14429 (2022).
  * Jian  _et al._ [2021] S.-K. Jian, B. Swingle, and Z.-Y. Xian, Complexity growth of operators in the syk model and in jt gravity, Journal of High Energy Physics 2021, 14 (2021).
  * Rubio-Abadal  _et al._ [2020] A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Floquet prethermalization in a bose-hubbard system, Phys. Rev. X 10, 021044 (2020).
  * Abanin  _et al._ [2017] D. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, A rigorous theory of many-body prethermalization for periodically driven and closed quantum systems, Communications in Mathematical Physics 354, 809 (2017).
  * Bhattacharjee  _et al._ [2023] B. Bhattacharjee, P. Nandy, and T. Pathak, Krylov complexity in large q and double-scaled syk model, Journal of High Energy Physics 2023, 99 (2023).
  * Bhattacharjee  _et al._ [2022] B. Bhattacharjee, X. Cao, P. Nandy, and T. Pathak, Krylov complexity in saddle-dominated scrambling, Journal of High Energy Physics 2022, 174 (2022).
  * Heveling  _et al._ [2022] R. Heveling, J. Wang, and J. Gemmer, Numerically probing the universal operator growth hypothesis, Phys. Rev. E 106, 014152 (2022).
  * Bouch [2015] G. Bouch, Complex-time singularity and locality estimates for quantum lattice systems, Journal of Mathematical Physics 56, 123303 (2015).
  * Cao [2021] X. Cao, A statistical mechanism for operator growth, Journal of Physics A: Mathematical and Theoretical 54, 144001 (2021).
  * Sandvik and Kurkijärvi [1991] A. W. Sandvik and J. Kurkijärvi, Quantum monte carlo simulation method for spin systems, Phys. Rev. B 43, 5950 (1991).
  * Sandvik [2010] A. W. Sandvik, Computational Studies of Quantum Spin Systems, AIP Conference Proceedings 1297, 135 (2010).
  * Wang and Landau [2001] F. Wang and D. P. Landau, Efficient, multiple-range random walk algorithm to calculate the density of states, Phys. Rev. Lett. 86, 2050 (2001).
  * Troyer  _et al._ [2003] M. Troyer, S. Wessel, and F. Alet, Flat histogram methods for quantum systems: Algorithms to overcome tunneling problems and calculate the free energy, Phys. Rev. Lett. 90, 120201 (2003).
  * Belardinelli  _et al._ [2008] R. E. Belardinelli, S. Manzi, and V. D. Pereyra, Analysis of the convergence of the $1/t$ and wang-landau algorithms in the calculation of multidimensional integrals, Phys. Rev. E 78, 067701 (2008).
  * [43] Attached supplemental material.
  * Araki [1969] H. Araki, Gibbs states of a one dimensional quantum lattice, Communications in Mathematical Physics 14, 120 (1969).
  * Ponte  _et al._ [2015] P. Ponte, Z. Papić, F. m. c. Huveneers, and D. A. Abanin, Many-body localization in periodically driven systems, Phys. Rev. Lett. 114, 140401 (2015).
  * Noh [2021] J. D. Noh, Operator growth in the transverse-field ising spin chain with integrability-breaking longitudinal field, Phys. Rev. E 104, 034112 (2021).
  * Lindner and Auerbach [2010] N. H. Lindner and A. Auerbach, Conductivity of hard core bosons: A paradigm of a bad metal, Phys. Rev. B 81, 054512 (2010).
  * Bhattacharyya  _et al._ [2024] S. Bhattacharyya, A. De, S. Gazit, and A. Auerbach, Metallic transport of hard-core bosons, Phys. Rev. B 109, 035117 (2024).
  * Viswanath and Müller [1990] V. S. Viswanath and G. Müller, Recursion method in quantum spin dynamics: The art of terminating a continued fraction, Journal of Applied Physics 67, 5486 (1990).
  * Viswanath and Müller [2008] V. Viswanath and G. Müller,  _The Recursion Method: Application to Many-Body Dynamics_, Lecture Notes in Physics Monographs (Springer Berlin Heidelberg, 2008).
  * Khait  _et al._ [2016] I. Khait, S. Gazit, N. Y. Yao, and A. Auerbach, Spin transport of weakly disordered heisenberg chain at infinite temperature, Phys. Rev. B 93, 224205 (2016).

Supplemental Materials:“Stochastic Sampling of Operator Growth Dynamics”

##  I The Monte Carlo algorithm

In this section, we outline our Marko chain Monte Carlo (MCMC) sampling
algorithm. We first derive the precise form of the partition function taking
into account the re-weighting procedure. Following this, we describe each of
the Monte Carlo moves and derive their transition probabilities based on the
detailed balance conditions.

We start by deriving Eq. 8 in the paper. Combining $\mathcal{Z}(\tau=1/J)$
from Eq. 3 with Eq. 6 in the main text, the partition function becomes:

$\mathcal{Z}(\tau=1/J)=\sum_{n=0}^{\infty}\sum_{\mathcal{S}_{2n}}\frac{1}{J^{2n}(2n)!}\prod_{k=0}^{2n-1}(O_{k+1}|\mathcal{L}_{i_{k},a_{k}}|O_{k})$
(S1)

where we recall that $\mathcal{S}_{2n}=\\{(i_{k},a_{k})\\}_{k=0}^{2n-1}$ is a
sequence of position-type pairs that specify the local Liouvillians
$\mathcal{L}_{i_{k},a_{k}}$. To proceed, we truncate this partition sum to
sequence $n\geq n_{\max}$. We also fix the Liouvillian sequences to be of
length $2n_{\max}$, by inserting $2n_{\max}-2n$ factors of the identity
Liouvillian $\mathcal{L}_{0,0}$ in all possible combinations. This procedure
however generates duplicate terms which are compensated for by dividing with
the combinatorial factor ${2n_{\max}}\choose{2n_{\max}-2n}$. This gives us the
desired expression of the truncated partition sum:

$G=\sum_{\mathcal{S}_{2n_{\max}}}\frac{1}{J^{2n}(2n)!}\times\frac{(2n_{\max}-2n)!(2n)!}{(2n_{\max})!}\prod_{k=0}^{2n_{\max}-1}(O_{k+1}|\mathcal{L}_{i_{k},a_{k}}|O_{k}).$
(S2)

As mentioned in the main text, the expansion order distribution associated
with $G$ falls typically exponentially on either side of its modal point. In
practice, this hinders the sampling of all the $n_{\max}+1$ sectors in an
ergodic way. To flatten out this distribution, we employ a variant of the
Wang-Landau algorithm [40, 42] that implements a re-weighting procedure.
Before explaining our scheme, it is illuminating to rewrite $G$ using a set of
generalized weights $\\{w_{n}\\}$:

$G^{\prime}=\sum_{\mathcal{S}_{2n_{\max}}}w_{n}\frac{(2n_{\max}-2n)!(2n)!}{(2n_{\max})!}\prod_{k=0}^{2n_{\max}-1}(O_{k+1}|\mathcal{L}_{i_{k},a_{k}}|O_{k})\,=\sum_{n}w_{n}\mu_{2n}.$
(S3)

The weights $\\{w_{n}\\}$ are initialized as $w_{n}=1/(J^{2n}(2n)!)$. In
practice one needs to only keep track of weight ratios between consecutive $n$
sectors $w_{n}/w_{n+1}$. The reweighting procedure works iteratively. At each
step, we take $N$ samples and determine the number of times $c_{n}$ that each
sector is visited by the MCMC. We define $c_{n}$ as:

$c_{n}=\sum_{i=1}^{N}\delta_{n,n_{i}}\,.$ (S4)

For a flat distribution, one expects $c_{n}\sim(n_{\max}+1)/N$. To bias our
distribution in this direction, we update the weights by setting
$w_{n}\rightarrow w_{n}N/(c_{n}(n_{\max}+1))$. The corresponding stored weight
ratios get updated as:
$w_{n}/w_{n+1}\rightarrow(w_{n}/w_{n+1})\times(c_{n+1}/c_{n})$. This
assignment penalizes sectors with $c_{n}>(n_{\max}+1)/N$ by reducing their
weight and reinforces weights for sectors with $c_{n}<(n_{\max}+1)/N$. This
procedure is carried out repeatedly until an approximately flat distribution
is obtained.

With generalized weights, we need to rescale the counts in Eq. 9 with
$\\{w_{n}\\}$ to obtain the moment ratios. Our observables now become:

$r_{n}^{2}=\frac{\mu_{2n+2}}{\mu_{2n}}=\frac{\left\langle\delta_{n^{\prime},n+1}\right\rangle_{n^{\prime}}}{\left\langle\delta_{n^{\prime},n}\right\rangle_{n^{\prime}}}\times\frac{w_{n}}{w_{n+1}}.$
(S5)

We now move our discussion towards the Monte Carlo moves. The moves are
designed to ensure ergodicity while only making local changes to the
configuration space. Below, we discuss the three classes of moves: add-drop,
swap and interact. For the rest of this section, we will work with a $d+1$
dimensional lattice with length $L$ and $2n_{\max}+1$ temporal slices. The
temporal axis is pictured in the vertical direction. While the schematic
diagrams show the moves on a $d=1$ lattice, they generally also apply to
higher dimensional cases. We will point out subtleties pertaining to $d>1$
dimensions when appropriate.

###  I.1 Add Drop

(a)  (b)  (c)

Figure S1:  Illustration of the add-drop moves with (a) $\mathcal{L}_{Z}$
Liouvillians (b) $\mathcal{L}_{X}$ Liouvillians (c) $\mathcal{L}_{ZZ}$
Liouvillians. The temporal direction is vertical.

The add-drop move is designed to allow transitioning between the different $n$
sectors by adding and removing a pair of identical local Liouvillians. As
illustrated in Fig. S1, the move is defined as follows:

  1. 1.

Randomly choose a time slice with uniform probability.

  2. 2.

If the local Liouvillians occupying both sides of the time slice are
identical, propose to remove both.

  3. 3.

Else, if the Liouvillian slots on both sides are unoccupied, randomly choose a
site with uniform probability and propose to add a pair of identical
Liouvillians onto the unoccupied slots at that site.

We now determine the transition probabilities. Let $i$ be a configuration in
the $n-1$ sector with the $(k-1)$-th and $k$-th Liouvillian slots unoccupied.
Let $f$ be a configuration in the $n$ sector that is identical to $i$ except
with these slots occupied at site $r$ by Liouvillians
$\mathcal{L}_{r_{k-1},a_{k-1}}=\mathcal{L}_{r_{k},a_{k}}=:\mathcal{L}_{r,a}$.
The two configurations are thus related by an add-drop move. Based on Eq. S3,
detailed balance imposes the constraint:

$\frac{T_{i\rightarrow f}A_{i\rightarrow f}}{T_{f\rightarrow i}A_{f\rightarrow
i}}=\frac{P_{f}}{P_{i}}=\frac{w_{n}(O_{k-1}|\mathcal{L}_{r,a}|O_{k})(O_{k}|\mathcal{L}_{r,a}|O_{k+1})(2n)(2n-1)}{w_{n-1}(2n_{\max}-2n+2)(2n_{\max}-2n+1)}\,.$
(S6)

Here, $P_{i},P_{f}$ are the respective weights, $T_{i\rightarrow
f},T_{f\rightarrow i}$ are probabilities to propose the add/drop transitions
and $A_{i\rightarrow f},A_{f\rightarrow i}$ are the acceptance probabilities
for the moves. Now, proposing an add operation involves choosing a time slice,
followed by choosing a site and then choosing a Liouvillian type; proposing a
drop amounts to simply choosing a time slice. Therefore the proposition
probabilities are

$\displaystyle T_{i\rightarrow
f}=\frac{1}{2n_{\max}+1}\times\frac{1}{L^{d}}\times\frac{(O_{k-1}|\mathcal{L}_{r,a}|O_{k})(O_{k}|\mathcal{L}_{r,a}|O_{k+1})}{\sigma}\,,\,T_{f\rightarrow
i}=\frac{1}{2n_{\max}+1}\,.$ (S7)

where we denoted $\sigma=4(dJ^{2}+(h^{X})^{2}+(h^{Z})^{2})$ as the sum of
weights of all possible Liouvillians at a site. Therefore, the following
“Metropolis” acceptance rates satisfy detailed balance:

$\displaystyle A_{i\rightarrow
f}=\text{min}\left(1,\frac{w_{n}L^{d}\sigma(2n)(2n-1)}{w_{n-1}(2n_{\max}-2n+2)(2n_{\max}-2n+1)}\right)\,,\,$
(S8) $\displaystyle A_{f\rightarrow
i}=\text{min}\left(1,\frac{w_{n-1}(2n_{\max}-2n+2)(2n_{\max}-2n+1)}{w_{n}L^{d}\sigma(2n)(2n-1)}\right)\,.$
(S9)

###  I.2 Swap

Figure S2:  Illustration of the swap move. We show two examples of
Liouvillians that commute with each other: (a) $\mathcal{L}_{Z}$ and
$\mathcal{L}_{ZZ}$ Liouvillians sharing a site. (b) Two $\mathcal{L}_{ZZ}$
Liouvillians sharing just one site. In addition to (a) and (b), identical
Liouvillians (those of the same kind, on the same site; and along the same
axis in case of $\mathcal{L}_{ZZ}$) and those on different sites with no
overlap, trivially commute.

The swap move permutes through different configurations in the same $n$ sector
using the commutation property of Liouvillian pairs, see Fig. Fig. S2 for an
illustration. Before discussing the steps, we first define the notion of
“ceiling” and “floor”. The ceiling (floor) of a Liouvillian is the temporally
closest Liouvillian above (below, respectively) it which doesn’t commute with
it. If a non-commuting Liouvillian is lacking, the ceiling/floor is then
defined to be the initial/final temporal boundary. Now, the swap steps are as
follows:

  1. 1.

Randomly, with uniform probability, choose a Liouvillian.

  2. 2.

Randomly, with uniform probability, choose a second Liouvillian slot amongst
those lying between the ceiling and floor of the first Liouvillian.

  3. 3.

If the newly chosen slot is empty, remove the Liouvillian from its existing
slot and place it in the new slot.

  4. 4.

Otherwise, that is, if the slot is occupied by another Liouvillian, check
whether the ceiling and floor of this second Liouvillian lie at an
intermediate temporal index between the indices of the two Liouvillians. If
that is not the case, swap the two Liouvillians.

Since the swap moves involve no weight change and everything is sampled
uniformly, it satisfies detailed balance automatically, that is, with $100\%$
acceptance rate.

###  I.3 Interact

Figure S3:  Illustration of the interact move. Two $\mathcal{L}_{ZZ}$
Liouvillians sharing only one site may “interact” to give a pair of
$\mathcal{L}_{Z}$ Liouvillians situated on next to nearest neighbour sites.

The add-drop move only allows Liouvillians to exist on a site in pairs. This
can however miss configurations with an odd number of bond Liouvillians on a
site. To enable this we introduce the interact move, as illustrated in Fig.
S3. It is to be noted that for $d>1$, the two interacting bonds can be along
different spatial axes. In such cases, two different $\mathcal{L}_{ZZ}$ pair
configurations can map to the same $\mathcal{L}_{Z}$ pair configuration. This
needs to be accounted for in detailed balance.

The interaction is implemented as:

  1. 1.

Randomly choose a time slice with uniform probability.

  2. 2.

If the Liouvillian on both sides of the time slice are of type
$\mathcal{L}_{ZZ}$ and share exactly one site, check if the shared site has a
$Z$ or $I$ operator on it. If so, replace the $\mathcal{L}_{ZZ}$ Liouvillians
with $\mathcal{L}_{Z}$ Liouvillians as in Fig. S3.

  3. 3.

If the Liouvillian slots on both sides host $\mathcal{L}_{Z}$ Liouvillians
placed at next to nearest neighbour sites, check if their common nearest
neighbouring site has a $Z$ or $I$ operator. If this is the case, replace the
$\mathcal{L}_{Z}$ Liouvillians with $\mathcal{L}_{ZZ}$ Liouvillians as in Fig.
S3.

Now, we derive the transition probabilities. Let $i$ be a configuration with
$\mathcal{L}_{ZZ}$ Liouvillians on the $k-1$-th and $k$-th time slots sharing
a single site $r$. Let $f$ be the configuration related to $i$ by the interact
move, hosting two $\mathcal{L}_{Z}$ Liouvillians. Now, detailed balance
implies:

$\frac{T_{i\rightarrow f}A_{i\rightarrow f}}{T_{f\rightarrow i}A_{f\rightarrow
i}}=\frac{P_{f}}{P_{i}}=\frac{(h^{Z})^{2}}{J^{2}}$ (S10)

where $P_{i},P_{f}$ are the respective weights, $T_{i\rightarrow
f},T_{f\rightarrow i}$ are proposition probabilities and $A_{i\rightarrow
f},A_{f\rightarrow i}$ the acceptance probabilities. We claim that the
proposition probabilities are as follows:

$T_{i\rightarrow f}=\frac{1}{2n_{\max}+1}\,,\,T_{f\rightarrow
i}=\frac{1}{(2n_{\max}+1)\alpha}\,,$ (S11)

where $\alpha=2$ if the $\mathcal{L}_{ZZ}$ Liouvillians in configuration $i$
are along different spatial directions (otherwise, $\alpha=1$). This is
because, proposing a $\mathcal{L}_{ZZ}\to\mathcal{L}_{Z}$ interact move
amounts to choosing a time slice. Meanwhile, to specify a
$\mathcal{L}_{Z}\to\mathcal{L}_{ZZ}$ interact move, we need choose a time
slice, and then pick a common neighboring site among the $\alpha$ choices.
Therefore, the following acceptance rates fulfill detailed balance:

$A_{i\rightarrow f}=\text{min}\left(1,\frac{(h^{Z})^{2}}{\alpha
J^{2}}\right)\,,\,A_{f\rightarrow i}=\text{min}\left(1,\frac{\alpha
J^{2}}{(h^{Z})^{2}}\right)\,.$ (S12)

##  II The $X$ Operator

(a)  (b)  (c)  (d)

Figure S4:  (a) $r_{n}$ vs $n$ for the $X$ operator with $J=h^{Z}=h^{X}=1$ in
1D. The inset focuses in on the first 30 moments, comparing our MCMC with
exact numerics. (b) $\Delta r_{n}$ and $\log(n)\Delta r_{n}$ vs $n$ for 1D.
(c) Same as (a) for 2D. The inset has the first 12 moments. (d) $\Delta r_{n}$
vs $n$ for 2D.

The statistical nature of infinite temperature operator growth dynamics in
non-integrable systems like the mixed field QIM implies that local operators
should typically give rise to the same generic behavior at a large enough
order $n$. Based on the Operator Growth Hypothesis [21], one expects local
operators to asymptotically grow linearly with the same slope. To check this,
we compute the same plots as Fig. 3 but for the $X$ operator in Fig. S4. We
see that for 1D $\log(n)\Delta r_{n}$ vs $n$ again gives us a constant curve
$\log(n)\Delta r_{n}\simeq 2.2(1)J$, the same as the $Z$ operator in the main
text. Similarly, in 2D we again get a plot consistent with linear growth with
slope $\simeq 1.26(4)J$. Another point to note is the non-universal nature of
the intermediate plateau in the 2D $\Delta r_{n}$ vs $n$ plot. We see that
here the plateau is much narrower than the one for the $Z$ operator as shown
in the main text.

##  III Additional move for the TFIM Limit

Figure S5: An example of a configuration in the $n=2$ sector that requires the
interact move to be generated. Temporal slices of the 2D lattice are indexed
by k. The first temporal slice is identical to the last because of the fixed
boundary conditions. The $\mathcal{L}_{ZZ}$ Liouvillians, shown as rectangular
boxes, transform local Pauli operators between slices.

In the $h^{Z}=0$ limit, the interact move breaks down. This is not a problem
in 1D as configurations with an odd number of identical $\mathcal{L}_{ZZ}$
Liouvillians also require the presence of $\mathcal{L}_{Z}$ Liouvillians and
thus do not exist in this limit. For $d>1$ it is possible for configurations
to have an odd number of identical $\mathcal{L}_{ZZ}$ Liouvillians without
$\mathcal{L}_{Z}$ Liouvillians. Fig. S5 shows a simple example with
$\mathcal{L}_{ZZ}$ Liouvillians _looping_ around a plaquette, satisfying the
boundary conditions. Note that each site hosts just one $\mathcal{L}_{ZZ}$
Liouvillian of a kind. One can generate larger _loops_ by removing
$\mathcal{L}_{ZZ}$ pairs on the common edges of overlapping smaller _loops_
using add-drop.

Figure S6:  Section of two configurations connected by the rotate move on a 2D
lattice. k corresponds to the temporal index. Note that both configurations
have the same temporal boundaries, implying the local nature of the move.

To access this part of the configuration space in absence of interact we
propose the rotate move. Fig. S6 shows that it is possible to rotate two
parallel $\mathcal{L}_{ZZ}$ Liouvillians sharing a plaquette to get another
pair of parallel $\mathcal{L}_{ZZ}$ Liouvillians on the same plaquette while
maintaining the same boundaries. Detailed balance for this move should be
trivial as the involved configurations have the same weight and the Metropolis
acceptance rate is $100\%$. It can be seen that a combination of add-drop and
rotate can generate the configuration in Fig. S5. We benchmark this move
against exact numerics in Fig. S7.

Figure S7:  Comparing our MCMC equipped with and without the rotate move
against exact computations for the 2D TFIM with the same boundary operator as
in Fig. S5.