††thanks: The first two authors contributed equally††thanks: These two authors
contributed equally

# Clustering of steady-state correlations in open systems with long-range
interactions

Andrew Y. Guo<EMAIL_ADDRESS>Joint Center for Quantum Information and Computer
Science, NIST/University of Maryland, College Park, MD 20742, USA  Joint
Quantum Institute, NIST/University of Maryland, College Park, MD 20742, USA
Simon Lieu  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, MD 20742, USA  Joint Quantum
Institute, NIST/University of Maryland, College Park, MD 20742, USA  Minh C.
Tran  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, MD 20742, USA  Joint Quantum
Institute, NIST/University of Maryland, College Park, MD 20742, USA  Center
for Theoretical Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA  Department of Physics, Harvard University,
Cambridge, Massachusetts 02138, USA  Alexey V. Gorshkov  Joint Center for
Quantum Information and Computer Science, NIST/University of Maryland, College
Park, MD 20742, USA  Joint Quantum Institute, NIST/University of Maryland,
College Park, MD 20742, USA

###### Abstract

Lieb-Robinson bounds are powerful analytical tools for constraining the
dynamic and static properties of non-relativistic quantum systems. Recently, a
complete picture for closed systems that evolve unitarily in time has been
achieved. In experimental systems, however, interactions with the environment
cannot generally be ignored, and the extension of Lieb-Robinson bounds to
dissipative systems which evolve non-unitarily in time remains an open
challenge. In this work, we prove two Lieb-Robinson bounds that constrain the
dynamics of open quantum systems with long-range interactions that decay as a
power-law in the distance between particles. Using a combination of these
Lieb-Robinson bounds and mixing bounds which arise from
“reversibility”—naturally satisfied for thermal environments—we prove the
clustering of correlations in the steady states of open quantum systems with
long-range interactions. Our work provides an initial step towards
constraining the steady-state entanglement structure for a broad class of
experimental platforms, and we highlight several open directions regarding the
application of Lieb-Robinson bounds to dissipative systems.

##  I Introduction

While the speed of information transfer is always bounded by the speed of
light, many quantum platforms operate in a non-relativistic regime where
typical velocities are far below this threshold. Nevertheless, the Schrödinger
equation admits fundamental limits to the rate at which correlations can
spread throughout the system. Such bounds are known as Lieb-Robinson bounds
and are connected to a diverse array of phenomena, including the decay of
correlations in the ground state [1], generation of topological order [2, 3],
efficiency of classical/quantum simulation [4, 5], hardness of bosonic
sampling tasks [6], heating rates in periodically driven Floquet systems [7,
8], and signatures of quantum chaos [9, 10].

To date, most formulations of Lieb-Robinson bounds apply to closed systems
that evolve via a unitary time-evolution operator. In such systems, recent
advances have proved tight information-transfer bounds for interaction ranges
that span the whole spectrum from local [11, 12] to highly non-local regimes
[5, 13, 14, 15], and have been saturated via explicit state-transfer protocols
[16, 17, 18]. While a complete picture for quantum information transfer has
emerged for closed quantum systems, the analogous question for systems that
evolve non-unitarily in time remains less well understood. For a broad range
of quantum platforms (including noisy quantum simulators), interactions with a
larger environment are unavoidable and must be taken into account to
accurately describe dynamics. While progress in this direction has been made
[19, 20, 21, 22, 23], the question of how the fundamental rate of information
transfer differs for systems that interact with some larger environment
remains unanswered.

Indeed, the notion of a Lieb-Robinson bound in an open system may seem a
priori surprising from the point of view of quantum trajectories [24]. In this
picture, in a time-step $dt$ the system’s wavefunction either evolves via a
non-unitary evolution operator, or a quantum jump discontinuously alters the
state. A specific trajectory belonging to a spatially-local Hamiltonian with
local dissipation can transfer information faster than the limit set by the
Hamiltonian’s Lieb-Robinson bound [25]. Intuitively, this is because
conditioning on measurements is an inherently nonlocal process. As an extreme
example, it is possible to create a highly-entangled (GHZ) state from a
product state in a time $t=\mathcal{O}\left\lparen 1\right\rparen$ using only
locally entangling gates and measurements, for a specific outcome of the
measurements [26]. This would violate the Lieb-Robinson bound for local
systems, which gives $t=\Omega\left\lparen r\right\rparen$ for distance $r$
[27]. After averaging over trajectories, the state of the system can be
represented via a density matrix $\rho$ which evolves via a master equation:
$d\rho/dt=\mathcal{L}(\rho)$. Subsequently, the notion of a Lieb-Robinson
bound is properly restored upon averaging over trajectories.

In this work, we make progress on the question of the fundamental rates of
information propagation in open systems by proving a broad class of Lieb-
Robinson bounds for systems with long-range interactions—specifically those
that decay as a power-law $1/r^{\alpha}$ in the distance $r$ between
particles, for some $\alpha>0$. Such power-law-decaying interactions feature
in experimental platforms relevant to quantum computation and simulation, such
as Rydberg atoms [28], trapped ion crystals [29, 30], polar molecules [31],
and nitrogen-vacancy color centers in diamond [32]. In all of these platforms,
interactions with a larger environment cannot be neglected, and a Markovian
description of system dynamics is often justified. In such systems, improved
understanding of the fundamental rates of information transfer has spurred the
development of optimal protocols for quantum information processing and state
transfer [16, 18].

In addition to bounding dynamics of open long-range systems, we use these
Lieb-Robinson bounds to constrain the entanglement structure of the
corresponding steady states. For closed systems, Lieb-Robinson bounds have
played an important role in proving rigorous statements on the decay of
correlations in gapped ground states [1]. This justifies the use of finite-
dimensional matrix-product-state representations of the ground state in one-
dimensional systems with local interactions [33]. In this work, we prove the
clustering of correlations in the steady states of open power-law systems,
which may serve as a first step towards establishing an area-law scaling of
entanglement for these systems, similar to what was done in Ref. [34] for the
closed case.

This paper is organized as follows: in Section II, we summarize the existing
Lieb-Robinson bounds for open long-range systems and present two new bounds
that are tighter for particular regimes of the power-law exponent $\alpha$.
The first yields a polynomial light cone for $\alpha>2d$, using a technique
pioneered in Ref. [5]. The second gives a linear light cone for $\alpha>3$ in
1D, using the method from Ref. [13]. In Section III, we also prove the
clustering of correlations in the steady states of open long-range systems.
Specifically, we provide bounds on the extent of the covariance correlations
and mutual information under certain assumptions on the Liouvillian mixing
times. We also prove a stability theorem for the stationary state under local
Liouvillian perturbations, generalizing the results of Ref. [22].

##  II Lieb-Robinson bounds for open long-range systems

In this section, we review the results of the previous best-known Lieb-
Robinson bounds for open long-range systems and state two new Lieb-Robinson
bounds.

As a general set-up, we consider evolution by a long-range Liouvillian
$\mathcal{L}(t)$ that acts on a lattice $\Lambda$ consisting of finite-level
systems at each site. We denote by $\mathcal{H}$ the finite-dimensional
Hilbert space representing all possible states of the system and by
$\mathcal{B(H)}$ the space of all operators on $\mathcal{H}$. For an operator
$O\in\mathcal{B(H)}$, we will be interested in how its expectation value
changes as a function of time: $\langle
O(t)\rangle=\tr[O(t)\rho]=\tr[O\rho(t)]$, where $\rho$ is the initial state of
the system, which evolves (in the Schrödinger picture) via
$\rho(t)=e^{\mathcal{L}t}\rho$. For these purposes, the time-evolution of $O$
can be expressed as $O(t)=e^{\mathcal{L}^{\dagger}t}O$, where
$\mathcal{L}^{\dagger}$ is the adjoint Lindblad superoperator, defined as

$\displaystyle\mathcal{L}^{\dagger}O=+i[H,O]+\sum_{i}\left[L_{i}^{\dagger}OL_{i}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},O\\}\right],$
(1)

where $H$ is the Hamiltonian and $L_{i}$ are Lindblad operators (also referred
to as “jump” operators) [35]. We emphasize that $O(t)$ is not equivalent to
the Heisenberg-Langevin time evolution for the operator $O$. For example, if
the system has a unique steady state, all operators $O(t)$ will be
proportional to the identity at long times:
$\lim_{t\rightarrow\infty}O(t)\sim\mathbb{I}$. Thus two operators that do not
commute at $t=0$ will start to commute at long times.

We will state the Lieb-Robinson-type bounds in this paper in terms of time-
independent Liouvillians. However, we note that the proofs can be generalized
with minor modifications to the case of time-dependent Liouvillians—i.e. those
for which both $H$ and $L_{i}$ are allowed to vary in time.

To impose the long-range condition on $\mathcal{L}$, we decompose it into
$\mathcal{L}=\sum_{Z\subset\Lambda}\mathcal{L}_{Z}$, where for any pair of
sites $i,j$, we have the condition

$\displaystyle\sum_{Z\ni
i,j}\|\mathcal{L}_{Z}\|_{\infty}\coloneqq\sup_{O\in\mathcal{B(H)}}\frac{\|\mathcal{L}_{Z}O\|}{\|O\|}\leq\frac{1}{\text{dist}(i,j)^{\alpha}},$
(2)

where $\|\cdot\|$ denotes the standard operator norm, or $\infty$-norm, and
$\|\cdot\|_{\infty}$ denotes the superoperator, or “$\infty\rightarrow\infty$”
norm (referred to as such because the second term in Eq. 2 uses the operator
$\infty$-norm in both the numerator and the denominator). Here
$\text{dist}(i,j)$ is the distance between $i$ and $j$, and $\alpha$ is the
positive real parameter that controls the long-ranged nature of the
interaction.

###  II.1 Prior work on open-system Lieb-Robinson bounds

In Ref. [23], Sweke _et al._ generalized the Lieb-Robinson bound in Ref. [1]
for $\alpha>d$ to open systems. Letting $A\in\mathcal{B}(X)$ be an operator
supported on $X$, $K_{Y}\in\mathbb{L}_{Y}$ be a Liouvillian supported on $Y$,
and $e^{\mathcal{L}^{\dagger}t}$ be the evolution under the adjoint
Liouvillian superoperator. The corresponding superoperator bound is:

$\displaystyle\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|\absolutevalue{X}\absolutevalue{Y}\left(\frac{e^{vt}-1}{r^{\alpha}}\right),$
(3)

where $r\coloneqq d(X,Y)$, and $C$ and $v$ are $\mathcal{O}\left\lparen
1\right\rparen$ constants. In the closed-system picture, the conventional
Lieb-Robinson-type bound can be recovered by choosing $K_{Y}$ such that
$K_{Y}(O_{X})=i[O_{X},O_{Y}]$ and replacing $\|K_{Y}\|_{\infty}$ with
$2\|O_{Y}\|$.

For this conventional bound, the velocity scales as $v\propto 2^{\alpha}$,
which diverges in the limit $\alpha\rightarrow\infty$. To recover the Lieb-
Robinson bound for short-range interacting systems, an improved bound is
required that uses a slight modification of the technique from Ref. [23]:

$\displaystyle\begin{split}\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|\absolutevalue{X}\absolutevalue{Y}\biggl{(}&\frac{e^{\tilde{v}t}}{[(1-\mu)r]^{\alpha}}\\\
&+e^{\tilde{v}t-\mu r}\biggr{)},\end{split}$ (4)

where $\mu\in(0,1)$ and $\tilde{v}$ are constants, and $\tilde{v}$ is
independent of $\alpha$. The closed-system version of this bound was first
proven in Ref. [36] for two-body interactions and later generalized to
$k$-body interactions in Ref. [8]. In Ref. [23], Sweke _et al._ also prove
Lieb-Robinson-type bounds for $\alpha\leq d$. For this regime of $\alpha$, one
needs to restrict to a finite-sized lattice, due to the energy being (in
general) non-additive for subsystems [37]. Denoting the system size of the
lattice by $N\coloneqq\absolutevalue{\Lambda}$, the combined strength $J$ of
the terms acting on a single site scales as $J=\Theta\left\lparen
N^{1-\alpha/d}\right\rparen$ for $\alpha<d$ and
$J=\Theta\left\lparen\log{N}\right\rparen$ for $\alpha=d$ [27, 10]. The bound
then becomes:

$\displaystyle\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|\absolutevalue{X}\absolutevalue{Y}\left(\frac{e^{Jt}-1}{Jr^{\alpha}}\right).$
(5)

The effective Lieb-Robinson velocity in this case diverges in the
thermodynamic limit, but is finite for all finite $N$.

###  II.2 Power-law light-cone bound for $\alpha>2d$

We prove a Lieb-Robinson bound for $\alpha>2d$ using the truncation-of-
unitaries-approach presented by Tran _et al._ [8]. The technique takes as
input the existing open-systems Lieb-Robinson bound in Eq. 4 and bootstraps it
to obtain a tighter bound:

$\displaystyle\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|\frac{t^{\alpha-d}}{r^{\alpha-2d}}.$ (6)

This bound yields the power-law light-cone contour
$t=r^{\frac{\alpha-2d}{\alpha-d}}$. The proof of this bound involves
approximating the time evolution of the operators by a sequence of operators
that span successively larger and larger subsets of the lattice, and bounding
the error of each successive approximation by the existing Lieb-Robinson
bound. We provide the full details of the derivation in Appendix E.

###  II.3 Linear light-cone bound for $d=1$, $\alpha>3$

Finally, we prove a bound with a linear light cone for open-long-range systems
with $\alpha>3$ in $d=1$ dimensions based on the techniques developed in Ref.
[13]. In the process, we tighten the tail of the Lieb-Robinson bound given in
that work from $1/r$ to approximately $1/r^{\alpha-2}$. The authors of Ref.
[13] proved the following bound for the closed-system dynamics of Hamiltonian
$H=\sum_{ij}H_{ij}$ consisting of two-body terms:

$\displaystyle\left\|[e^{iHt}Ae^{-iHt},B]\right\|\leq
C\left\|A\right\|\left\|B\right\|\frac{t}{r},$ (7)

where $B\in\mathcal{B}(Y)$ is an operator supported on $Y$. Likewise assuming
a two-body Liouvillian $\mathcal{L}=\sum_{ij}\mathcal{L}_{ij}$, we obtain the
following open-systems bound:

$\displaystyle\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
C\left\|K_{Y}\right\|_{\infty}\left\|A\right\|\frac{t}{r^{\alpha-2-o(1)}},$
(8)

where the $o(1)$ denotes some constant that can be made arbitrarily small. The
result yields a linear light cone $t\gtrsim r$ for all $\alpha>3$. The proof
roughly proceeds by expanding out the evolution operator
$e^{\mathcal{L}^{\dagger}t}$ into a series of products of Liouvillian terms
$\mathcal{L}_{ij}$. For each term in the series, we select out a subsequence
of terms that move the operator forward (i.e. towards $Y$) and integrate out
the other terms. By only taking into account the contributions from the terms
in the subsequences, we are able to obtain a tighter Lieb-Robinson bound. We
provide the mathematical details of the proof in Appendix F.

##  III Bounds on correlations in the steady states of open long-range systems

In this section, we prove the clustering of correlations in the steady states
of open long-range systems. We first state a lemma that describes how to use a
modified version of the Lieb-Robinson bounds stated in the previous section to
bound how far operators can spread under evolution by the (adjoint)
Liouvillian $\mathcal{L}^{\dagger}$. Specifically, we give a bound on the
error of approximating the time-evolution of an operator $A$ supported on a
site $X\in\Lambda$ by a truncated adjoint Liouvillian that only acts on ball
of radius $r$ centered on a site $X\in\Lambda$ (see Fig. 1).

Figure 1: The evolution of an operator $A$ initially supported on $X$ by an
adjoint Liouvillian $\mathcal{L}^{\dagger}$ can be approximated by the same
operator evolved by the truncated version of the Liouvillian,
$\tilde{\mathcal{L}}^{\dagger}$, supported on a ball of radius $r$ around $X$,
up to an error given by $C(r,t)$.

######  Lemma 1 (Bounds on the error incurred by approximating of time-evolved
operators by local ones).

Let $A$ be an operator initially supported on a site $X\in\Lambda$ and let
$\tilde{\mathcal{L}}$ be the restriction of the long-range Liouvillian
$\mathcal{L}$ to the ball of radius $r$ centered on $X$. Let $\tilde{A}(t)$ be
the time-evolved version of $A$ under $\tilde{\mathcal{L}}^{\dagger}$. Then
the error in the approximation of $A(t)$ by $\tilde{A}(t)$ is bounded by

$\displaystyle\|A(t)-\tilde{A}(t)\|\leq K\|A\|\,\mathcal{C}(r,t),$ (9)

where $K$ is some constant, and $\mathcal{C}(r,t)$ is a modified version of
the standard Lieb-Robinson bound adapted to the problem of locally
approximating time-evolved operators. In the large-$r$ limit, the tightest-
known bounds for open systems with long-range interactions with $\alpha>d$
scale asymptotically as

$\displaystyle\mathcal{C}(r,t)\propto\begin{cases}\displaystyle
e^{vt}/r^{\alpha-d},&\alpha>d,\\\ \displaystyle
t^{\alpha-d+1}/r^{\alpha-3d},&\alpha>3d,\\\ \displaystyle
t^{2}/r^{\alpha-3},&\alpha>3,d=1.\end{cases}$ (10)

For $\alpha\leq d$, the bounds also depend on the system size of the lattice
$N\coloneqq|\Lambda|$. When $r\propto N^{1/d}$, the bounds scale as follows:

$\mathcal{C}(N,t)\propto\begin{cases}\displaystyle\frac{e^{\Theta\left\lparen
N^{1-\alpha/d}\right\rparen t}-1}{\Theta\left\lparen
N^{1-\alpha/d}\right\rparen},&\alpha<d,\\\
\displaystyle\frac{e^{\Theta\left\lparen\log(N)\right\rparen
t}-1}{\Theta\left\lparen\log(N)\right\rparen},&\alpha=d.\end{cases}$ (11)

This concludes the statement of Lemma 1.

The proof of Lemma 1 follows straightforwardly from the open-system Lieb-
Robinson bounds detailed in Section II. In particular, the three lines of Eq.
10 follow from Eq. 4, Eq. 6, and Eq. 8, respectively, while Eq. 11 comes from
Eq. 5. In Appendix A, we provide the details of the derivation of the bounds
in Lemma 1. We now proceed to derive the bounds on clustering of correlations
in the steady states of gapped, reversible Liouvillians.

###  III.1 Bound on covariance correlations

In this first section, we show how open-system Lieb-Robinson bounds constrain
the correlations in the steady state of a Liouvillian $\mathcal{L}$ with
dissipative gap $\lambda$.

The dissipative gap $\lambda>0$ is defined as the magnitude of the least-
negative non-zero real part of an eigenvalue of $\mathcal{L}$. (Throughout
this work, we shall also assume that the Liouvillian is primitive, i.e. it has
a unique steady state such that $\mathcal{L}$ has one eigenvalue of zero.) In
addition to the Lieb-Robinson bounds, we will also appeal to certain “mixing
bounds” which describe how fast arbitrary initial states (or various
correlation functions) converge to the steady state.

For the mixing bounds, we need to impose “reversibility” on the Liouvillian.
We say that a Liouvillian $\mathcal{L}$ is $s$-reversible if there exists some
operator $\sigma$ such that
$\Gamma_{s}\mathcal{L}^{\dagger}=\mathcal{L}\Gamma_{s}$ is satisfied; the
superoperator $\Gamma_{s}$ is defined via
$\Gamma_{s}(f)=(\sigma^{s}f\sigma^{1-s}+\sigma^{1-s}f\sigma^{s})/2$ and
$s\in[0,1]$. For $s$-reversible Liouvillians, it is easy to see that $\sigma$
is the steady state of $\mathcal{L}$ (since
$\mathcal{L}^{\dagger}(\mathbb{I})=0$). A sufficient condition for a
Liouvillian to satisfy $s$-reversibility (for all $s$) is if the dissipators
$L_{i}$ satisfy a detailed-balance condition (and the Hamiltonian is zero,
$H=0$), which is naturally obeyed for systems coupled to a thermal bath [22].
(More explicitly, the detailed-balance condition is satisfied if dissipators
come in energy raising/lowering pairs with respect to some effective
Hamiltonian $\bar{H}$—for example, if $[\bar{H},L_{\pm}]=\pm\omega L_{\pm}$
and $|L_{-}|/|L_{+}|=\exp(2\beta\omega)$ where $\beta^{-1}$ is an effective
temperature.)

Returning to the topic of correlations, we let $\rho$ be a quantum state
defined on the lattice $\Lambda$. We are interested in the _covariance
correlation_ between non-overlapping $X,Y\in\Lambda$:

$\displaystyle
T_{\rho}(X:Y):=\text{sup}_{\left\|f\right\|=\left\|g\right\|=1}|\text{Tr}[(f\otimes
g)(\rho_{XY}-\rho_{X}\otimes\rho_{Y})]|,$ (12)

where $f$ and $g$ are Hermitian operators with $f$ supported on region $X$ and
$g$ supported on region $Y$, and where $\rho_{X}$ is the reduced density
matrix constructed from $\rho$ by tracing over the complement of $X$. Our goal
is to bound this correlation function in terms of $\lambda$ and the distance
between $X$ and $Y$.

We follow Theorem 9 in Ref. [22]. Let $\sigma$ be the steady state of the
Liouvillian. From the right-hand side of Eq. 12, we define

$\text{Cov}_{\sigma}(f,g)\coloneqq\frac{1}{2}\text{Tr}[(fg+gf)\sigma]-\text{Tr}[f\sigma]\text{Tr}[g\sigma],$
(13)

which is equivalent to the term inside the sup (because $f$ and $g$ commute).
Now we use the bound (which follows directly from the triangle inequality)

$|\text{Cov}_{\sigma}(f,g)|\leq|\text{Cov}_{\sigma}(f_{t},g_{t})|+|\text{Cov}_{\sigma}(f,g)-\text{Cov}_{\sigma}(f_{t},g_{t})|.$
(14)

Here $f_{t}$ and $g_{t}$ are the time-evolved versions of $f$ and $g$ under
$\mathcal{L}^{\dagger}$. This step allows us to relate a static covariance to
time-dependent quantities; we will use dynamical bounds to constrain the form
of the latter, then pick an optimal time which maximally bounds the static
covariance.

The first term on the right is constrained by the variance bound for
$s$-reversible, primitive Liouvillians (see Appendix B)

$|\text{Cov}_{\sigma}(f_{t},g_{t})|\leq
4\left\|f\right\|\left\|g\right\|e^{-2\lambda t},$ (15)

where $\lambda$ is the dissipative gap of $\mathcal{L}$. Intuitively, this
relationship can be understood as follows: the operators $f_{t},g_{t}$ both
evolve (in time) toward an operator that is proportional to the identity, so
the covariance between them will eventually tend to zero as a function of
time. The rate at which this occurs is set by the dissipative gap of the
system.

To bound the second term, we use the relation $\Tr[\sigma f_{t}]=\Tr[\sigma
f]$, which holds for all observables $f$. This gives:

$\displaystyle|\text{Cov}_{\sigma}(f,g)-\text{Cov}_{\sigma}(f_{t},g_{t})|$
(16) $\displaystyle=\frac{1}{2}(|\Tr[(fg-f_{t}g_{t})\sigma]+\Tr[(gf-
g_{t}f_{t})\sigma]|)$ (17)
$\displaystyle=\frac{1}{2}(|\Tr[((fg)_{t}-f_{t}g_{t})\sigma]+\Tr[((gf)_{t}-g_{t}f_{t})\sigma]|)$
(18)
$\displaystyle\leq\frac{1}{2}(\|(fg)_{t}-f_{t}g_{t}\|+\|(gf)_{t}-g_{t}f_{t}\|)$
(19) $\displaystyle\leq K\left\|f\right\|\left\|g\right\|\mathcal{C}(r,t),$
(20)

where $r\coloneqq d(X,Y)$. We obtain the inequality in the final line using
the open-system Lieb-Robinson bounds $\mathcal{C}(r,t)$ given in Lemma 1.
Specifically, we use the following Lemma, which is itself a restatement of
Corollary 7 in Ref. [22]:

######  Lemma 2 (Time-evolution of spatially separated observables).

Take two operators $A$ and $B$ supported on $X,Y\in\Lambda$ respectively such
that $r\coloneqq d(X,Y)$, and let $A(t)=e^{\mathcal{L}^{\dagger}t}A$ and
$B(t)=e^{\mathcal{L}^{\dagger}t}B$ be their time-evolution under the adjoint
Liouvillian $\mathcal{L}^{\dagger}$. We also define
$(AB)(t)=e^{\mathcal{L}^{\dagger}t}(AB)$. Then the following bound holds:

$\displaystyle\|(AB)(t)-A(t)B(t)\|\leq K\|A\|\|B\|\mathcal{C}(r,t),$ (21)

where $\mathcal{C}(r,t)$ is given by Lemma 1 and $K$ is some constant that
depends on lattice parameters.

Lemma 2 bounds the difference between operators that evolve together in the
Heisenberg picture as opposed to evolving separately. Again we emphasize that
$A(t)$, $B(t)$, and $(AB)(t)$ are not equivalent to Heisenberg-Langevin
evolution, a fact that is at the core of this bound. We defer the short proof
of Lemma 2 to Appendix C and move on to proving the bound on the covariance
correlations.

######  Theorem 1 (Bounds on steady-state covariance correlations).

Consider Hermitian operators $f,g$ which are supported on two non-overlapping
subsets $X$ and $Y$ of the $d$-dimensional cubic lattice $\Lambda$, and let
$\mathcal{L}$ be an $s$-reversible Liouvillian with stationary state $\sigma$
and dissipative gap $\lambda$ that satisfies the conditions in Eq. 2. Then
there exists a constant $c>0$ which only depends on $\lambda,v$ such that

$T_{\sigma}(X:Y)\leq\begin{cases}\displaystyle
c\left(r^{\alpha-d}\right)^{\frac{-2\lambda}{v+2\lambda}},&\alpha>d,\\\
\displaystyle c\frac{\log(r)^{\alpha-d+1}}{r^{\alpha-3d}},&\alpha>3d,\\\
\displaystyle c\frac{\log(r)^{2}}{r^{\alpha-3}},&\alpha>3,d=1.\end{cases}$
(22)

###### Proof.

From our previous analysis [see Eqs. (15) and (20)] on the covariance
correlation in Eq. 14, we have

$|\text{Cov}_{\sigma}(f,g)|\leq
4\left\|f\right\|\left\|g\right\|\left(e^{-2\lambda
t}+\frac{K}{4}\mathcal{C}(r,t)\right).$ (23)

To obtain the tightest bound, we minimize with respect to $t$ the function

$h(t)=e^{-\lambda^{\prime}t}+K^{\prime}\mathcal{C}(r,t),$ (24)

where $\lambda^{\prime}=2\lambda,K^{\prime}=K/4$.

We will perform this minimization exactly for the first case in Eq. 22, for
which $\mathcal{C}(r,t)$ is given by the first line of Eq. 10; for the other
cases, we instead use an approximation to the optimal ansatz, which allows us
to obtain an analytical expression for the bound. Setting $dh/dt=0$ in Eq. 24
leads to a minimum at time

$\bar{t}=-\left(\frac{1}{\lambda^{\prime}+v}\right)\log\left(\frac{K^{\prime}v}{\lambda
r^{\alpha-d}}\right).$ (25)

This implies a minimum:

$\displaystyle h(\bar{t})$
$\displaystyle=\left(\frac{K^{\prime}v}{\lambda^{\prime}r^{\alpha-d}}\right)^{\frac{\lambda^{\prime}}{\lambda^{\prime}+v}}+\frac{K^{\prime}}{r^{\alpha-d}}\left(\frac{K^{\prime}v}{\lambda^{\prime}r^{\alpha-d}}\right)^{\frac{-v}{\lambda^{\prime}+v}}$
$\displaystyle\leq c\left(r^{\alpha-d}\right)^{\frac{-2\lambda}{v+2\lambda}}$
(26)

for some constant $c$ which depends on $\lambda,v,K$. Taking the supremum over
$f,g$ gives the bound on $T_{\sigma}(X:Y)$ for $\alpha>d$ in the first line of
Eq. (22).

For the other two cases, we use the ansatz $t^{*}=1+\log(r^{\beta})$. Since
the bound in the second line of Eq. 10 scales as $\mathcal{C}(r,t)\propto
t^{\alpha-d+1}/r^{\alpha-3d}$ for all $t$, we have

$\displaystyle h(t^{*})$
$\displaystyle=e^{-\lambda(1+\log(r^{\beta}))}+K\frac{(1+\log(r^{\beta}))^{\alpha-d+1}}{r^{\alpha-3d}}$
$\displaystyle=\frac{e^{-\lambda}}{r^{\lambda\beta}}+K\frac{(\beta\log(r))^{\alpha-d+1}}{r^{\alpha-3d}}+\mathcal{O}\left\lparen\frac{\log^{\alpha-d}(r)}{r^{\alpha-3d}}\right\rparen.$
(27)

We choose $\beta=(\alpha-3d)/\lambda$, which is positive for $\alpha>3d$. This
gives the ultimate bound of

$\displaystyle h(t^{*})$
$\displaystyle=\frac{e^{-\lambda}+K\left(\frac{\alpha-3d}{\lambda}\log
r\right)^{\alpha-d+1}}{r^{\alpha-3d}}+\mathcal{O}\left\lparen\frac{\log^{\alpha-d}(r)}{r^{\alpha-3d}}\right\rparen$
$\displaystyle=K\left(\frac{\alpha-3d}{\lambda}\right)^{\alpha-d+1}\frac{\log^{\alpha-d+1}(r)}{r^{\alpha-3d}}$
$\displaystyle\quad\qquad+\mathcal{O}\left\lparen\frac{\log^{\alpha-d}(r)}{r^{\alpha-3d}}\right\rparen,$
(28)

which proves the second line of Eq. (22). For the $d=1$ case in the last line
of Eq. (22), the argument proceeds similarly, but we obtain a slightly better
scaling in the logarithmic factor. ∎

Here we discuss the scaling of the bounds in Eq. 22, which is depicted in Fig.
2. The effective exponent of the $1/r$-scaling of the bound for $\alpha>d$ is
$\alpha^{\prime}\equiv(\alpha-d)\frac{2\lambda}{v+2\lambda}$, as compared to
$\tilde{}\alpha\equiv\alpha-3$ for $\alpha>3d$ (neglecting terms doubly
logarithmic in $r$). Since $\alpha^{\prime}$ decreases as a function of $v$,
the former bound becomes looser for larger $v$. In more detail, if we let
$x=\frac{v}{\lambda}$, then $\alpha^{\prime}<\tilde{}\alpha$ for all
$\alpha>\frac{(3x+4)d}{x}$. In the limit of $x\rightarrow\infty$,
$\tilde{}\alpha$ is tighter for all $\alpha>3d$. Thus, for large enough
$\alpha$ and $v$, the power-law light-cone bounds [second line in Eq. (10),
which in turn comes from Eq. (6)] give asymptotically tighter bounds on the
clustering of covariance correlations than the logarithmic light-cone bound
[first line in Eq. (10), which in turn comes from Eq. (4)].

Figure 2: A log-log plot of the tails of the bounds on the various connected
correlation functions in Theorems 1, 2, and 3 for $d=1$. We include the
exponentially decaying tail from the short-range interactions case (red curve)
for comparison. For the power-law decaying bounds, we have different scaling
exponents of the power-law tails for the bound for $\alpha>1$ (blue curve) and
the bound for $\alpha>3$ (green curve). For a given choice of
$x=\frac{v}{\lambda}$, the relative positioning of the curves holds for all
$\alpha>\frac{3x+4}{x}$. In the limit $v\gg\lambda$, the picture holds for all
$\alpha>3$.

###  III.2 Stability result and mutual information bound

In this section, we will use the aforementioned bounds to constrain steady-
state properties of open systems with power-law interactions. In addition to
the newly-derived Lieb-Robinson bounds, we will appeal to a “mixing bound”
which provides an upper bound to how fast an arbitrary initial state will
converge to the steady state. The following mixing bound was derived in Ref.
[38], and generalizes the mixing bound of classical Markov chains to quantum
semigroups:

######  Lemma 3.

Consider a primitive Liouvillian $\mathcal{L}$ that has a full-rank steady
state $\sigma$, and is $\frac{1}{2}$-reversible. Then an arbitrary initial
state $\rho$ will converge to $\sigma$ at a rate bounded by

$\lVert\rho(t)-\sigma\lVert_{1}\leq\sqrt{2\log(\lVert\sigma^{-1}\lVert)}e^{-\beta
t},$ (29)

where $\beta$ is called the log-Sobolev constant associated with
$\mathcal{L}$.

Intuition can be gained by considering an “infinite-temperature” steady state
$\sigma=\mathbb{I}/d_{H}$ where $d_{H}$ is the dimension of the Hilbert space.
The mixing bound above states that the coefficient in front of the exponential
will scale as $\sqrt{\log(d_{H})}$, i.e. it will increase with the dimension
of the Hilbert space. This is because the convergence toward the unique steady
state from an arbitrary initial state can be slow if the dimension of the
Hilbert space is large.

######  Theorem 2 (Effect of perturbations on reduced steady-state density
matrix).

Let $X,Y$ be two non-overlapping subsets of a $d$-dimensional cubic lattice
$\Lambda$. Let $\mathcal{L}_{1}$ be a primitive and $\frac{1}{2}$-reversible
Liouvillian with log-Sobolev constant $\beta$, and let $\mathcal{L}_{2}$ be a
Liouvillian perturbation, acting trivially outside of $X$. Let $\rho$ be the
stationary state of $\mathcal{L}_{1}$, and let $\sigma$ be the stationary
state of $\mathcal{L}_{1}+\mathcal{L}_{2}$. Then,

$\lVert\rho_{Y}-\sigma_{Y}\lVert_{1}\leq\begin{cases}\displaystyle
c\log(\lVert\rho^{-1}\lVert)^{\frac{1}{2}}\left(\frac{1}{r^{\alpha-d}}\right)^{\frac{2\beta}{v+2\beta}},&\alpha>d,\\\
\displaystyle
c\log(\lVert\rho^{-1}\lVert)^{\frac{1}{2}}\frac{\log(r)^{\alpha-d+1}}{r^{\alpha-3d}},&\alpha>3d,\\\
\displaystyle
c\log(\lVert\rho^{-1}\lVert)^{\frac{1}{2}}\frac{\log(r)^{2}}{r^{\alpha-3}},&\alpha>3,d=1,\end{cases}$
(30)

where $c$ is a constant and $r$ is the distance between $X$ and $Y$.

The theorem basically says that the effects of local perturbations in the
Liouvillian will not be felt significantly by the steady state of the system
at sufficiently distant locations. We prove the theorem by first introducing a
time-evolved state to interpolate between the two steady states. This allows
us to use a combination of mixing-time and Lieb-Robinson bounds to restrict
the terms in this bound. Then we apply the same minimization procedure used in
Theorem 1 for the covariance-correlations bound to arrive at the stated bounds
in Eq. 30 [each of which follow directly from the three cases in Eq. 10]. We
defer the proof of this result, which is similar to the proof of Theorem 1, to
Appendix D.

We now prove a bound on the mutual information in the steady state. The mutual
information between two regions $A,B$ is defined as

$I_{\rho}(A:B)=S(\rho_{AB}||\rho_{A}\otimes\rho_{B}),$ (31)

where $S(\rho||\sigma)=\tr[\rho(\log\rho-\log\sigma)]$ is the relative
entropy. The following theorem holds.

######  Theorem 3 (Clustering of mutual information).

Let $A,B$ be two non-overlapping subsets of a $d$-dimensional cubic lattice
$\Lambda$. Let $\mathcal{L}$ be a primitive and $\frac{1}{2}$-reversible
Liouvillian with log-Sobolev constant $\beta$. Let $\rho$ be the stationary
state of $\mathcal{L}$. Then the mutual information between the two regions
$I_{\rho}(A:B)$ is bounded by

$I_{\rho}(A:B)\leq\begin{cases}\displaystyle
c\log(\lVert\rho^{-1}\lVert)^{\frac{3}{2}}\left(\frac{1}{r^{\alpha-d}}\right)^{\frac{2\beta}{v+2\beta}},&\alpha>d,\\\
\displaystyle
c\log(\lVert\rho^{-1}\lVert)^{\frac{3}{2}}\frac{\log(r)^{\alpha-d+1}}{r^{\alpha-3d}},&\alpha>3d,\\\
\displaystyle
c\log(\lVert\rho^{-1}\lVert)^{\frac{3}{2}}\frac{\log(r)^{2}}{r^{\alpha-3}},&\alpha>3,d=1,\end{cases}$
(32)

where $c$ is a constant and $r$ is the distance between $A,B$.

The significance of this result is that the mutual-information correlations in
the steady state of an open long-range system decay as a power-law in the
distance between regions. This bound, which relies on the existence of the
log-Sobolev constant, is tighter than the naive bound that would result from
simply applying the bound on the covariance correlation in Theorem 1 to
$I_{\rho}(A:B)$.

###### Proof.

We define the semi-group $\tilde{\mathcal{L}}$ to be the terms in
$\mathcal{L}$ that act entirely within balls of radius $r/2$ centered around
$A$ and $B$, and let $\sigma$ be the steady state of $\tilde{\mathcal{L}}$.
Simple manipulations imply:

$\displaystyle I_{\rho}(A:B)$
$\displaystyle=-S(\rho_{AB})+S(\rho_{A})+S(\rho_{B})$ (33)
$\displaystyle\leq-S(\rho_{AB})-\tr[\rho_{A}\log\sigma_{A}]-\tr[\rho_{B}\log\sigma_{B}]$
(34) $\displaystyle=S(\rho_{AB}||\sigma_{A}\otimes\sigma_{B}).$ (35)

where we have used $S(\rho||\sigma)\geq 0$ to obtain the inequality. The RHS
further satisfies the inequality:

$S(\rho_{AB}||\sigma_{A}\otimes\sigma_{B})\leq\log(||\rho_{AB}^{-1}||)||\rho_{AB}-\sigma_{A}\otimes\sigma_{B}||_{1},$
(36)

which is a standard result (c.f. Eq. (36) in Ref. [22]). From here, we can
apply the bounds in Theorem 2, using $\mathcal{L}_{1}=\tilde{\mathcal{L}}$,
$\mathcal{L}_{2}=\mathcal{L}-\tilde{\mathcal{L}}$, $Y=A\cup B$, and
$X=\Lambda\setminus Y$. ∎

##  IV Summary and outlook

In this work, we have proven generalized Lieb-Robinson bounds which constrain
the dynamics of open, Markovian systems with power-law interactions and used
them to constrain correlations in the steady state.

We comment briefly on the tightness of the bounds derived in this work.
Intuitively, one might expect that the presence of dissipation should lead to
tighter Lieb-Robinson bounds for open systems than for their closed
counterparts, since the presence of decoherence from a bath might limit the
speed of quantum information transfer. In this work, we have generalized the
proof of Lieb-Robinson bounds from closed system dynamics to Markovian
evolution (a priori, such bounds need not exist for Markovian dynamics).
However, our bounds only depend on interaction range and the dimension of the
lattice. Any bound that only depends on these two inputs cannot be tighter
than the corresponding closed-system Lieb-Robinson bound, since the latter is
a special case of former. As such, the saturating protocols for closed systems
[17, 18] can be used to saturate open Lieb-Robinson bounds such as those
uncovered in Lemma 1. In the future, it would be interesting to add another
degree of freedom into formulations of open Lieb-Robinson bounds: the
dissipative gap. (Some progress has been made in showing that Lieb-Robinson
velocities can get tighter in dissipative systems [20].) In principle, it
might be possible to derive Lieb-Robinson bounds that reduce to closed-system
ones when the dissipative gap is zero, and get tighter in the presence of non-
zero dissipation. Then one can develop protocols that saturate the
dissipative-gap-dependent bounds. Another question in this direction is
whether the conditional evolution generated via a non-Hermitian Hamiltonian
can also exhibit a dissipative-gap-dependent Lieb-Robinson bound that reduces
to the conventional one in the dissipationless limit.

Setting aside the idea of a Lieb-Robinson bound that depends on the
dissipative gap, there is still the question of generalizing the best-known
closed-system bounds to Markovian evolution. In particular, the recent Lieb-
Robinson bounds in Refs. [14] and [15] both provide opportunities for
generalization to open systems. Such a result would likely require a
modification of the interaction-picture technique first developed in [39] and
used in both subsequent works to open-system dynamics. Generalizing these
bounds would directly lead to tighter bounds on operator spreading in Lemma 1
and allow us to prove tighter bounds on correlation clustering in steady
states (Theorem 1 and Theorem 3)

Another way to probe the tightness of the steady-state correlation bounds
derived in this work would be to improve the mixing bounds, which currently
require the open system to be in thermal equilibrium. It would be interesting
to derive more general mixing bounds which also apply to systems that are out
of thermal equilibrium.

One of the salient applications of Lieb-Robinson bounds is in rigorous proofs
on the stability of the spectral gap in topologically ordered quantum matter.
For example, Ref. [3] used closed-system Lieb-Robinson bounds to show that
spatially local perturbations will not close energy gaps in the toric code,
thus leading to phase stability against arbitrary local noise. Can we use a
similar approach to show that local perturbations will not close the
dissipative gap in a topologically-ordered open system? A robust qubit steady-
state structure would be useful toward the quest of passive quantum error
correction [40].

Lieb-Robinson bounds can be used to prove area-law entanglement scaling in the
ground state of one-dimensional systems with local interactions [33]. This
result helps to rigorously justify the validity of the matrix-product state
ansatz for the ground state of such systems. For closed systems with power-law
interactions, Lieb-Robinson bounds can be used to further extend area-law
scaling to certain broad classes of systems [34]. Do the results presented in
this paper have similar implications for area-law scaling of the steady state?
This would have direct implications for the matrix-product operator ansatz in
modeling open systems.

Finally, the Lieb-Robinson-type bounds we proved apply for the operator, or
$\infty$-norm. However, there exists a hierarchy of Lieb-Robinson-like bounds
that have the potential to be tighter for certain information processing tasks
such as scrambling and transferring a quantum state of a local subsystem
without knowledge of the initial state of the rest of the system. These bounds
can use other norms such as the Frobenius norm defined by
$\|O\|_{F}=\sqrt{\Tr{O^{\dagger}O}}$ [17, 41, 42, 43] or apply to free-
particle systems [10, 17]. It would be interesting to generalize these bounds
to open systems as well.

###### Acknowledgements.

We would like to thank Andrew Lucas, Adam Ehrenberg, and Chris Baldwin for
helpful discussions. AYG is supported by the NSF Graduate Research Fellowship
Program under Grant No. DGE-1840340. SL was supported by the NIST NRC Research
Postdoctoral Associateship Award. AYG, MCT, and AVG acknowledge funding by the
DoE ASCR Quantum Testbed Pathfinder program (award No. DE-SC0019040), DoE
Quantum Systems Accelerator, NSF PFCQC program, AFOSR MURI, DoE ASCR
Accelerated Research in Quantum Computing program (award No. DE-SC0020312),
ARO MURI, AFOSR, U.S. Department of Energy Award No. DE-SC0019449, and DARPA
SAVaNT ADVENT. MCT acknowledges additional support from the Department of
Energy, Office of Advanced Scientific Computing Research through the QOALAS
program (grant 17-020469) and the Quantum Algorithms and Machine Learning
grant from NTT, number AGMT DTD 9/24/20.

## References

  * Hastings and Koma [2006] M. B. Hastings and T. Koma, Spectral Gap and Exponential Decay of Correlations, Comm. Math. Phys. 265, 781 (2006).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-Robinson Bounds and the Generation of Correlations and Topological Quantum Order, Phys. Rev. Lett. 97, 050401 (2006).
  * Bravyi  _et al._ [2010] S. Bravyi, M. B. Hastings, and S. Michalakis, Topological quantum order: Stability under local perturbations, J. Math. Phys. 51, 093512 (2010).
  * Osborne [2006] T. J. Osborne, Efficient Approximation of the Dynamics of One-Dimensional Quantum Spin Systems, Phys. Rev. Lett. 97, 157202 (2006).
  * Tran  _et al._ [2019a] M. C. Tran, A. Y. Guo, Y. Su, J. R. Garrison, Z. Eldredge, M. Foss-Feig, A. M. Childs, and A. V. Gorshkov, Locality and digital quantum simulation of power-law interactions, Phys. Rev. X 9, 031006 (2019a).
  * Deshpande  _et al._ [2018] A. Deshpande, B. Fefferman, M. C. Tran, M. Foss-Feig, and A. V. Gorshkov, Dynamical phase transitions in sampling complexity, Phys. Rev. Lett. 121, 030501 (2018).
  * Abanin  _et al._ [2015] D. A. Abanin, W. De Roeck, and F. Huveneers, Exponentially Slow Heating in Periodically Driven Many-Body Systems, PRL 115, 256803 (2015).
  * Tran  _et al._ [2019b] M. C. Tran, A. Ehrenberg, A. Y. Guo, P. Titum, D. A. Abanin, and A. V. Gorshkov, Locality and heating in periodically driven, power-law-interacting systems, Phys. Rev. A 100, 052103 (2019b).
  * Lashkari  _et al._ [2013] N. Lashkari, D. Stanford, M. Hastings, T. Osborne, and P. Hayden, Towards the fast scrambling conjecture, J. High Energy Phys. 2013 (4), 22.
  * Guo  _et al._ [2020] A. Y. Guo, M. C. Tran, A. M. Childs, A. V. Gorshkov, and Z.-X. Gong, Signaling and Scrambling with Strongly Long-Range Interactions, Phys. Rev. A 102, 010401 (2020).
  * Chen and Lucas [2021] C.-F. Chen and A. Lucas, Operator growth bounds from graph theory, Communications in Mathematical Physics 385, 1273 (2021).
  * Wang and Hazzard [2020] Z. Wang and K. R. Hazzard, Tightening the lieb-robinson bound in locally interacting systems, PRX Quantum 1, 010303 (2020).
  * Chen and Lucas [2019] C.-F. Chen and A. Lucas, Finite speed of quantum scrambling with long range interactions, Phys. Rev. Lett. 123, 250605 (2019).
  * Kuwahara and Saito [2020] T. Kuwahara and K. Saito, Strictly linear light cones in long-range interacting systems of arbitrary dimensions, Phys. Rev. X 10, 031010 (2020).
  * Tran  _et al._ [2021a] M. C. Tran, A. Y. Guo, C. L. Baldwin, A. Ehrenberg, A. V. Gorshkov, and A. Lucas, The lieb-robinson light cone for power-law interactions (2021a), arXiv:2103.15828 .
  * Eldredge  _et al._ [2017] Z. Eldredge, Z.-X. Gong, J. T. Young, A. H. Moosavian, M. Foss-Feig, and A. V. Gorshkov, Fast Quantum State Transfer and Entanglement Renormalization Using Long-Range Interactions, Phys. Rev. Lett. 119, 170503 (2017).
  * Tran  _et al._ [2020] M. C. Tran, C.-F. Chen, A. Ehrenberg, A. Y. Guo, A. Deshpande, Y. Hong, Z.-X. Gong, A. V. Gorshkov, and A. Lucas, Hierarchy of Linear Light Cones with Long-Range Interactions, Phys. Rev. X 10, 031009 (2020).
  * Tran  _et al._ [2021b] M. C. Tran, A. Y. Guo, A. Deshpande, A. Lucas, and A. V. Gorshkov, Optimal state transfer and entanglement generation in power-law interacting systems, Phys. Rev. X 11, 031016 (2021b).
  * Poulin [2010] D. Poulin, Lieb-robinson bound and locality for general markovian quantum dynamics, Phys. Rev. Lett. 104, 190401 (2010).
  * Descamps [2013] B. Descamps, Asymptotically decreasing lieb-robinson velocity for a class of dissipative quantum dynamics, J. Math. Phys. 54, 092202 (2013).
  * Cubitt  _et al._ [2015] T. S. Cubitt, A. Lucia, S. Michalakis, and D. Perez-Garcia, Stability of local quantum dissipative systems, Commun. Math. Phys. 337, 1275 (2015).
  * Kastoryano and Eisert [2013] M. J. Kastoryano and J. Eisert, Rapid mixing implies exponential decay of correlations, J. Math. Phys. 54, 102201 (2013).
  * Sweke  _et al._ [2019] R. Sweke, J. Eisert, and M. Kastner, Lieb–robinson bounds for open quantum systems with long-ranged interactions, J. Phys. A Math. Theor. 52, 424003 (2019).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, The quantum-jump approach to dissipative dynamics in quantum optics, Rev. Mod. Phys. 70, 101 (1998).
  * Ashida and Ueda [2018] Y. Ashida and M. Ueda, Full-counting many-particle dynamics: Nonlocal and chiral propagation of correlations, Phys. Rev. Lett. 120, 185301 (2018).
  * Pham and Svore [2013] P. Pham and K. M. Svore, A 2d nearest-neighbor quantum architecture for factoring in polylogarithmic depth, Quantum Info. Comput. 13, 937–962 (2013).
  * [27] In this paper, we use the following standard asymptotic notations: the “big-$\mathcal{O}$” notation $f(x)=\mathcal{O}\left\lparen g(x)\right\rparen$ indicates the existence of constants $c$ and $N$ such that $0\leq f(x)\leq cg(x)$ for all $x\geq N$; the “big-$\Omega$” notation $f(x)=\Omega\left\lparen g(x)\right\rparen$ indicates that $g(x)=\mathcal{O}\left\lparen f(x)\right\rparen$, and the “big-$\Theta$” notation $f(x)=\Theta(g(x))$ when both $f(x)=\mathcal{O}\left\lparen g(x)\right\rparen$ and $f(x)=\Omega\left\lparen g(x)\right\rparen$.
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Britton  _et al._ [2012] J. W. Britton, B. C. Sawyer, A. C. Keith, C. C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Engineered two-dimensional Ising interactions in a trapped-ion quantum simulator with hundreds of spins, Nature 484, 489 (2012).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Yan  _et al._ [2013] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Observation of dipolar spin-exchange interactions with lattice-confined polar molecules, Nature 501, 521 (2013).
  * Yao  _et al._ [2012] N. Y. Yao, L. Jiang, A. V. Gorshkov, P. C. Maurer, G. Giedke, J. I. Cirac, and M. D. Lukin, Scalable architecture for a room temperature solid-state quantum information processor, Nat. Commun. 3, 800 EP (2012).
  * Hastings [2007] M. B. Hastings, An area law for one-dimensional quantum systems, J. Stat. Mech: Theory Exp. 2007, P08024 (2007).
  * Gong  _et al._ [2017] Z.-X. Gong, M. Foss-Feig, F. G. S. L. Brandão, and A. V. Gorshkov, Entanglement area laws for long-range interacting systems, Phys. Rev. Lett. 119, 050501 (2017).
  * Breuer and Petruccione [2010] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford Univ. Press, 2010).
  * Gong  _et al._ [2014] Z.-X. Gong, M. Foss-Feig, S. Michalakis, and A. V. Gorshkov, Persistence of locality in systems with power-law interactions, Phys. Rev. Lett. 113, 030602 (2014).
  * Dauxois  _et al._ [2002] T. Dauxois, S. Ruffo, E. Arimondo, and M. Wilkens,  _Dynamics and Thermodynamics of Systems with Long-Range Interactions: An Introduction_ (Springer Berlin Heidelberg, 2002) pp. 1–19.
  * Kastoryano and Temme [2013] M. J. Kastoryano and K. Temme, Quantum logarithmic sobolev inequalities and rapid mixing, J. Math. Phys. 54, 052202 (2013).
  * Foss-Feig  _et al._ [2015] M. Foss-Feig, Z.-X. Gong, C. W. Clark, and A. V. Gorshkov, Nearly linear light cones in long-range interacting quantum systems, Phys. Rev. Lett. 114, 157201 (2015).
  * Lidar  _et al._ [1998] D. A. Lidar, I. L. Chuang, and K. B. Whaley, Decoherence-free subspaces for quantum computation, Phys. Rev. Lett. 81, 2594 (1998).
  * Kuwahara and Saito [2021] T. Kuwahara and K. Saito, Absence of fast scrambling in thermodynamically stable long-range interacting systems, Phys. Rev. Lett. 126, 030604 (2021).
  * Yin and Lucas [2020] C. Yin and A. Lucas, Bound on quantum scrambling with all-to-all interactions, Phys. Rev. A 102, 022402 (2020).
  * Chen and Lucas [2021] C.-F. Chen and A. Lucas, Optimal Frobenius light cone in spin chains with power-law interactions (2021), arXiv:2105.09960 [quant-ph] .
  * Mostafazadeh [2002] A. Mostafazadeh, Pseudo-hermiticity versus pt symmetry: The necessary condition for the reality of the spectrum of a non-hermitian hamiltonian, J. Math. Phys. 43, 205 (2002).

##  Appendix A Bounds on the error incurred by approximating time-evolved
operators by local ones

Here we use the open-system Lieb-Robinson bounds described in Section II of
the main text to derive the scalings in Lemma 1. Recall that $\tilde{A}(t)$ is
the evolution of the operator $A$ under the Liouvillian $\tilde{\mathcal{L}}$,
the restriction of the long-range Liouvillian $\mathcal{L}$ to
$\mathcal{B}_{r}(X)$, the ball of radius $r$ centered on $X$, for time $t$. We
bound the difference between $A(t)$, which is $A$ evolved by the full
Liouvillian, and $\tilde{A}(t)$ as follows:

$\displaystyle\left\|A(t)-\tilde{A}(t)\right\|$
$\displaystyle=\left\|\int_{0}^{t}\frac{d}{ds}\left[e^{\mathcal{L}^{\dagger}(t-s)}e^{\tilde{\mathcal{L}}^{\dagger}s}A\right]\,\text{d}s\right\|$
(37)
$\displaystyle=\left\|\int_{0}^{t}e^{\mathcal{L}^{\dagger}(t-s)}(\mathcal{L}^{\dagger}-\tilde{\mathcal{L}}^{\dagger})\tilde{A}(s)\,\text{d}s\right\|$
(38)
$\displaystyle\leq\int_{0}^{t}\sum_{j:\text{dist}(j,X)>r}\sum_{i:\text{dist}(i,X)\leq
r}\left\|\mathcal{L}^{\dagger}_{ij}\tilde{A}(s)\right\|\,\text{d}s.$ (39)

In order to bound $\left\|\mathcal{L}^{\dagger}_{ij}\tilde{A}(s)\right\|$, we
turn to the open-system Lieb-Robinson bounds discussed in Section II. Each
line of Eq. 10 and Eq. 11 will correspond to plugging in one of those bounds.
For ease of reference, we reproduce the scalings here:

$\mathcal{C}(r,t)\propto\begin{cases}\displaystyle\frac{e^{\Theta\left\lparen
N^{1-\alpha/d}\right\rparen t}-1}{\Theta\left\lparen
N^{1-\alpha/d}\right\rparen},&\alpha<d,\\\
\displaystyle\frac{e^{\Theta\left\lparen\log(N)\right\rparen
t}-1}{\Theta\left\lparen\log(N)\right\rparen},&\alpha=d.\\\
\displaystyle\frac{e^{vt}}{r^{\alpha-d}},&\alpha>d,\\\
\displaystyle\frac{t^{\alpha-d+1}}{r^{\alpha-3d}},&\alpha>3d,\\\
\displaystyle\frac{t^{2}}{r^{\alpha-3}},&\alpha>3,d=1.\end{cases}$ (40)

The calculations will be similar for each bound, so we will only demonstrate
the result of inserting the power-law light cone bound from Eq. 6 into Eq. 39:

$\displaystyle\left\|A(t)-\tilde{A}(t)\right\|$ $\displaystyle\leq
C\left\|A\right\|\int_{0}^{t}\text{d}s\sum_{j:\text{dist}(j,X)>r}\sum_{i:d(i,X)\leq
r}\left\|\mathcal{L}^{\dagger}_{ij}\right\|\frac{s^{\alpha-d}}{\text{dist}(i,X)^{\alpha-2d}}$
(41) $\displaystyle\leq
C\left\|A\right\|\int_{0}^{t}\text{d}s\sum_{j:\text{dist}(j,X)>r}\sum_{i:\text{dist}(i,X)\leq
r}\frac{1}{\text{dist}(i,j)^{\alpha}}\frac{s^{\alpha-d}}{\text{dist}(i,X)^{\alpha-2d}}$
(42) $\displaystyle\leq
C^{\prime}\left\|A\right\|\int_{0}^{t}\text{d}s\sum_{j:\text{dist}(j,X)>r}\frac{s^{\alpha-d}}{\text{dist}(j,X)^{\alpha-2d}}$
(43) $\displaystyle\leq
C^{\prime\prime}\left\|A\right\|\frac{t^{\alpha-d+1}}{r^{\alpha-3d}}.$ (44)

This yields the expression in the fourth line of Eq. 40. Performing the same
operations for the other bounds gives the other terms in Eq. 40: the first and
second lines come from Eq. 5; the third line comes from Eq. 4, and the last
line comes from Eq. 8.

##  Appendix B Variance bound for reversible Liouvillians

Here we provide a derivation of the covariance bound used in Eq. (15). We show
that $s$-reversibility is important for this bound to hold. We define the
variance of an observable $f$ in the steady state $\sigma$ as
$\text{Var}[f]=\text{Tr}[f^{2}\sigma]-\text{Tr}[f\sigma]^{2}$, which is real
and positive. We wish to find a bound for $\text{Var}[f_{t}]$ for the time-
evolved observable $f_{t}=e^{\mathcal{L}^{\dagger}t}f$.

The Liouvillian is a non-Hermitian superoperator, which means that each
eigenvalue has right and left eigenoperators:

$\mathcal{L}(r_{i})=\lambda_{i}r_{i},\qquad\mathcal{L}^{\dagger}(l_{i})=\lambda_{i}^{*}l_{i}.$
(45)

From the structure of the adjoint Liouvillian ($\mathcal{L}^{\dagger}$), it is
clear that $\mathcal{L}^{\dagger}(\mathbb{I})=0$, where $\mathbb{I}$ is the
identity operator. This implies that one of the eigenvalues $\lambda_{0}$ is
zero, and the corresponding right eigenoperator $\sigma$ is called the steady
state and satisfies $\mathcal{L}(\sigma)=0$ and
$e^{\mathcal{L}t}(\sigma)=\sigma$. The eigenoperators are “bi-orthonormal” via
the Hilbert-Schmidt inner product:
$\text{Tr}[l_{i}^{\dagger}r_{j}]=\text{Tr}[r_{i}^{\dagger}l_{j}]=\delta_{ij}$.

We define the superoperator
$\Gamma_{s}(f)=(\sigma^{s}f\sigma^{1-s}+\sigma^{1-s}f\sigma^{s})/2$ where
$s\in[0,1]$ and $\sigma$ is a full-rank, Hermitian operator with positive
eigenvalues. We say that a Liouvillian is $s$-reversible for some $s\in[0,1]$
if $\Gamma_{s}\mathcal{L}^{\dagger}=\mathcal{L}\Gamma_{s}$. By acting both
sides on the operator $\mathbb{I}$, we see that $\sigma$ is the steady state,
i.e. that $\mathcal{L}(\sigma)=0$. Imposing reversibility implies that the
spectrum must be real because the Liouvillian is pseudo-Hermitian with a
positive-definite metric [44].

The dynamics preserves Hermiticity of a density matrix, which implies that
$\mathcal{L}(f^{\dagger})=[\mathcal{L}(f)]^{\dagger}$, and the same for the
adjoint:
$\mathcal{L}^{\dagger}(f^{\dagger})=[\mathcal{L}^{\dagger}(f)]^{\dagger}$,
where $f$ is an arbitrary operator. This implies that (right and left)
eigenoperators with real eigenvalues must be Hermitian. For $s$-reversible
Liouvillians, the entire spectrum is real, which implies that all
eigenoperators are Hermitian.

Ref. [22] derives a bound for the time-evolved variance in a $s$-reversible
system:

$\text{Var}[f_{t}]\leq e^{-2\lambda_{1}t}\text{Var}[f(t=0)],$ (46)

where $\\{-\lambda_{i}\\}$ is the real, non-positive spectrum of
$\mathcal{L}$, sorted from smallest to largest magnitude with
$\lambda_{0}=0,\lambda_{1}>0$. ($\lambda=\lambda_{1}$, i.e. the dissipative
gap.) Here we derive this bound using the properties of the eigenoperators of
$\mathcal{L}$.

Consider a general Hermitian operator which we write in terms of left
eigenoperators

$f=\sum_{j}c_{j}l_{j}\Rightarrow f_{t}=\sum_{j}c_{j}e^{-\lambda_{j}t}l_{j},$
(47)

where $c_{j}$ are real because $f$ is Hermitian. Noting that
$\text{Tr}[l_{j}\sigma]=0$ for $j\neq 0$, we find

$\text{Var}[f_{t}]=\text{Tr}\left[\left(\sum_{j\neq
0}c_{j}e^{-\lambda_{j}t}l_{j}\right)^{2}\sigma\right]=\sum_{j\neq
0}c_{j}^{2}e^{-2\lambda_{j}t},$ (48)

where in the last equality we have used
$\text{Tr}[l_{i}r_{j}]=\text{Tr}[l_{i}\Gamma_{s}(l_{j})]=\delta_{ij}$. From
this, it is easy to see that $\lambda_{i>1}\geq\lambda_{1}$ implies the bound
Eq. (46). For the more general case of a complex spectrum, it is not clear how
to repeat the derivation above. We therefore find that $s$-reversibility is
sufficient for the bound to hold. (It is unclear whether $s$-reversibility is
necessary for the bound.)

Given the bound Eq. (46), one can repeat the steps outlined in Eqs. [49-55] in
Ref. [22] to obtain the bound used in Eq. (15) of the main text. For
completeness, we include these steps below:

$\displaystyle|\text{Cov}_{\sigma}(f_{t},g_{t})|$
$\displaystyle\leq\sqrt{\text{Var}(f_{t})\text{Var}(g_{t})}$ (49)
$\displaystyle\leq e^{-2t\lambda_{1}}\sqrt{\text{Var}(f)\text{Var}(g)}.$ (50)

The inequality in (49) is due to Holder’s inequality. The variance can be
bounded by

$\displaystyle\sqrt{\text{Var}(f)}$
$\displaystyle=\sqrt{\Tr[\sigma(f-\Tr[\sigma f])^{2}]}$ (51)
$\displaystyle\leq\sqrt{\lVert(f-\Tr[\sigma f])^{2}\lVert}$ (52)
$\displaystyle\leq\lVert f-\Tr[\sigma f]\lVert$ (53) $\displaystyle\leq\lVert
f\lVert+|\Tr[\sigma f]|$ (54) $\displaystyle\leq 2\lVert f\lVert.$ (55)

Putting together (50) and (55) leads to the desired bound (where
$\lambda=\lambda_{1}$, i.e. the dissipative gap):

$|\text{Cov}_{\sigma}(f_{t},g_{t})|\leq
4\left\|f\right\|\left\|g\right\|e^{-2\lambda t},$ (56)

which matches Eq. 15 in the main text.

##  Appendix C Bound on the difference between two operators evolving
separately versus evolving together

In this section, we provide the proof of the bound in Lemma 2. We restate the
lemma here for convenience:

######  Lemma 2.

Take two operators $A$ and $B$ supported on single sites $X,Y\in\Lambda$
respectively such that $r\coloneqq d(X,Y)$, and let
$A(t)=e^{\mathcal{L}^{\dagger}t}A$ and $B(t)=e^{\mathcal{L}^{\dagger}t}B$ be
their time-evolution under the Liouvillian superoperator
$\mathcal{L}^{\dagger}$. We also define
$(AB)(t)=e^{\mathcal{L}^{\dagger}t}(AB)$. Then the following bound holds:

$\displaystyle\|(AB)(t)-A(t)B(t)\|\leq K^{\prime}\|A\|\|B\|\mathcal{C}(r,t),$
(57)

where $\mathcal{C}(r,t)$ is given by the Lieb-Robinson-type bound
corresponding to the system in question (see Lemma 1) and $K^{\prime}$ is some
constant that depends on lattice parameters.

###### Proof.

We define the semi-group $\tilde{\mathcal{L}}^{\dagger}$ to be the terms in
$\mathcal{L}^{\dagger}$ that act entirely within balls of radius $r/2$
centered around $X$ and $Y$. Then, let $\tilde{A}(t)$ be the time-evolved
version of $A$ under $\tilde{\mathcal{L}}^{\dagger}$ and likewise for
$\tilde{B}(t)$. By definition, this implies that
$\tilde{A}(t)\tilde{B}(t)=(\widetilde{AB})(t)$. We then get

$\displaystyle\|(AB)(t)-A(t)B(t)\|\leq\|(AB)(t)-(\widetilde{AB})(t)\|+\|A(t)B(t)-\tilde{A}(t)\tilde{B}(t)\|.$
(58)

The first term on the RHS of Eq. 58 may be bounded by the Lieb-Robinson bound
stated in Lemma 1 (for an operator that is initially supported on two sites
instead of one). The second term can be bounded by

$\displaystyle\|A(t)B(t)-\tilde{A}(t)\tilde{B}(t)\|$
$\displaystyle\leq\|A(t)(B(t)-\tilde{B}(t))\|+\|(A(t)-\tilde{A}(t))\tilde{B}(t)\|$
(59) $\displaystyle\leq\|A\|\|B(t)-\tilde{B}(t)\|+\|A(t)-\tilde{A}(t)\|\|B\|,$
(60)

using $\|A(t)\|\leq\|A\|$ and the submultiplicativity of the operator norm.
Using the Lieb-Robinson bound again, we get

$\|(AB)(t)-A(t)B(t)\|\leq 2K\|A\|\|B\|\mathcal{C}(r,t),$ (61)

which is the same as Eq. 21 in the main text. ∎

##  Appendix D Effect of perturbations on reduced steady-state density matrix

In this section, we provide the proof of Theorem 2. The argument hews closely
to that of Lemma 11 in Ref. [22], but uses the Lieb-Robinson bounds for open
long-range systems given in the main text.

######  Theorem 2.

Let $X,Y$ be two non-overlapping subsets of a $d$-dimensional cubic lattice
$\Lambda$. Let $\mathcal{L}$ be a primitive and $s$-reversible Liouvillian
with log-Sobolev constant $\beta$, and let $\mathcal{Q}$ be a local
Liouvillian perturbation, acting trivially outside of $X$. Let $\rho$ be the
stationary state of $\mathcal{L}$, and let $\sigma$ be the stationary state of
$\mathcal{L}+\mathcal{Q}$. Then,

$\lVert\rho_{Y}-\sigma_{Y}\lVert_{1}\leq\begin{cases}c\log(\lVert\rho^{-1}\lVert)^{\frac{1}{2}}\left(\frac{1}{r^{\alpha-d}}\right)^{\frac{2\beta}{v+2\beta}},&\alpha>d,\\\
c\log(\lVert\rho^{-1}\lVert)^{\frac{1}{2}}\frac{\log(r)^{\alpha-d+1}}{r^{\alpha-3d}},&\alpha>3d,\\\
c\log(\lVert\rho^{-1}\lVert)^{\frac{1}{2}}\frac{\log(r)^{2}}{r^{\alpha-3}},&\alpha>3,\end{cases}$
(62)

where $c$ is some constant, and $r$ is the distance between $X$ and $Y$.

###### Proof.

We use the following definition of the trace norm:

$\frac{1}{2}\lVert\rho-\sigma\lVert_{1}=\max_{0\leq
A\leq\mathbb{I}}\tr[A(\rho-\sigma)],$ (63)

for positive semi-definite $A$. This implies

$\lVert\rho_{Y}-\sigma_{Y}\lVert_{1}=2\tr[(A_{Y}\otimes\mathbb{I}_{Y^{c}})(\rho-\sigma)],$
(64)

where $A_{Y}=\tr_{Y^{c}}[\text{argmax}_{0\leq
A\leq\mathbb{I}}\tr[A(\rho-\sigma)]]$. We use the triangle inequality

$\displaystyle\tr[(A_{Y}\otimes\mathbb{I}_{Y^{c}})(\rho-\sigma)]$
$\displaystyle=\tr[(A_{Y}\otimes\mathbb{I}_{Y^{c}})\left[(e^{\mathcal{L}t}-e^{(\mathcal{L}+\mathcal{Q})t})(\phi)+(\sigma-e^{(\mathcal{L}+\mathcal{Q})t}(\phi))+(e^{\mathcal{L}t}(\phi)-\rho)\right]]$
(65)
$\displaystyle\leq\tr[(A_{Y}\otimes\mathbb{I}_{Y^{c}})(e^{\mathcal{L}t}-e^{(\mathcal{L}+\mathcal{Q})t})(\phi)]+\frac{1}{2}\lVert\tr_{Y^{c}}[\sigma-e^{(\mathcal{L}+\mathcal{Q})t}(\phi)]\lVert_{1}+\frac{1}{2}\lVert\tr_{Y^{c}}[e^{\mathcal{L}t}(\phi)-\rho]\lVert_{1},$
(66)

where $\phi$ is an arbitrary state. Note that we have introduced two time-
evolved operators in this step. We will now use a combination of mixing bounds
and Lieb-Robinson bounds to restrict the RHS. The last term is bounded via the
log-Sobolev bound:

$\frac{1}{2}\lVert\tr_{Y^{c}}[e^{\mathcal{L}t}(\phi)-\rho]\lVert_{1}\leq\left(\frac{1}{2}\log(\lVert\rho^{-1}\lVert)\right)^{\frac{1}{2}}e^{-\beta
t}.$ (67)

This is basically an upper bound on how fast an arbitrary initial state must
converge towards the steady state. The second term in Eq. (65) can be bounded
using a combination of Lieb-Robinson bounds and the log-Sobolev bound:

$\displaystyle\frac{1}{2}\lVert\tr_{Y^{c}}[\sigma-e^{(\mathcal{L}+\mathcal{Q})t}(\phi)]\lVert_{1}$
$\displaystyle=\tr[A_{Y}e^{(\mathcal{L}+\mathcal{Q})t}(\sigma-\phi)]$ (68)
$\displaystyle=\tr[e^{(\mathcal{L}^{\dagger}+\mathcal{Q}^{\dagger})t}(A_{Y})(\sigma-\phi)]$
(69)
$\displaystyle\leq\tr[(e^{(\mathcal{L}^{\dagger}+\mathcal{Q}^{\dagger})t}-e^{\mathcal{L}^{\dagger}t})(A_{Y})(\sigma-\phi)]+\tr[e^{\mathcal{L}^{\dagger}t}(A_{Y})(\sigma-\phi)].$
(70)

The last term can again be bounded via the log-Sobolev bound:

$\displaystyle\tr[e^{\mathcal{L}^{\dagger}t}(A_{Y})(\sigma-\phi)]$
$\displaystyle\leq\frac{1}{2}\lVert
e^{\mathcal{L}t}(\sigma-\phi)\lVert_{1}\leq\left(2\log(\lVert\rho^{-1}\lVert)\right)^{\frac{1}{2}}e^{-\beta
t}.$ (71)

The first term can be bounded via the Lieb-Robinson bound:

$\displaystyle\tr[(e^{(\mathcal{L}^{\dagger}+\mathcal{Q}^{\dagger})t}-e^{\mathcal{L}^{\dagger}t})(A_{Y})(\sigma-\phi)]$
$\displaystyle\leq\tr[(e^{(\mathcal{L}^{\dagger}+\mathcal{Q}^{\dagger})t}-e^{\mathcal{L}^{\dagger}t})(A_{Y})]\lVert\sigma-\phi\lVert_{1}$
(72) $\displaystyle\leq
2\tr[(e^{(\mathcal{L}^{\dagger}+\mathcal{Q}^{\dagger})t}-e^{\mathcal{L}^{\dagger}t})(A_{Y})]$
(73) $\displaystyle\leq
2\tr[(e^{(\mathcal{L}^{\dagger}+\mathcal{Q}^{\dagger})t}-e^{\mathcal{L}_{X^{c}}^{\dagger}t})(A_{Y})]+2\tr[(e^{\mathcal{L}_{X^{c}}^{\dagger}t}-e^{\mathcal{L}^{\dagger}t})(A_{Y})]$
(74) $\displaystyle\leq K\lVert
A_{Y}\lVert\,\mathcal{\mathcal{missing}}{C}(r,t),$ (75)

where $\mathcal{L}_{X^{c}}$ is the Liouvillian restricted to terms that do not
intersect $X$. $K$ is an arbitrary constant, and
$\mathcal{\mathcal{missing}}{C}(r,t)$ is the Lieb-Robinson bound stated in
Lemma 1.

The first term in Eq. (65) can be bounded using the Lieb-Robinson approach
above. Gathering all the bounds together leads to

$\lVert\rho_{Y}-\sigma_{Y}\lVert_{1}\leq
K_{1}\left(\log(\lVert\rho^{-1}\lVert)\right)^{\frac{1}{2}}e^{-\beta
t}+K_{2}\,\mathcal{\mathcal{missing}}{C}(r,t)$ (76)

for arbitrary constants $K_{1},K_{2}$. We wish to pick a time $t$ that
minimizes the RHS. We now note that the RHS has the same functional form as
the function that we needed to minimize for the covariance correlation bound.
Repeating the minimization procedure outlined in Theorem 1, we arrive at the
stated bounds in Eq. 30 of the main text. ∎

##  Appendix E Generalization of the Tran _et al._ bound to open long-range
systems

Here we provide the derivation of the open-systems Lieb-Robinson bound in Eq.
6. We use the generalization of the Hastings & Koma bound to open systems, as
described in [23]. Let $K_{Y}\in\mathbb{L}_{Y}$ be a Liouvillian with support
contained in $Y$ and $\tau(t)\equiv e^{\mathcal{L}^{\dagger}t}$ be the
backwards time-evolution operator. The corresponding superoperator bound is

$\displaystyle\mathcal{C}(r,t)\equiv\left\|K_{Y}(\tau(t)A)\right\|\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|\absolutevalue{X}\absolutevalue{Y}\frac{e^{vt}-1}{r^{\alpha}},$
(77)

If the supports of operators $K_{Y}$ and $A$ are not constant, then summing
Eq. 77 over the sites in those supports gives a bound of

$\displaystyle\mathcal{C}(r,t)\leq\|K_{Y}\|_{\infty}\left\|A\right\|\phi(Y)\frac{e^{vt}}{r^{\alpha-d-1}},$
(78)

where $\phi(Y)$ denotes the boundary of $Y$. For simplicity, we will later
write this bound in the form

$\displaystyle C(r,t)\leq\|K_{Y}\|_{\infty}\left\|A\right\|\phi(Y)f(r,t).$
(79)

To derive the open-systems Lieb-Robinson bound in Eq. 6, we follow the proof
in Tran _et al._ [8]. We first divide up the time interval $[0,t]$ into $M$
timesteps of size $\Delta t\equiv t/M$ and let $t_{i}=it/M$ for $i=0,\dots,M$.
For brevity, we denote by $\tau_{i}\equiv\tau(t_{M-i},t_{M-i+1})$ the time-
evolution operator from time $t_{M-i}$ to $t_{M-i+1}$. We can decompose the
evolution of $A$ by $\tau(t)$ into $M$ timesteps:

$\displaystyle\tau(t)A=\tau_{M}\tau_{M-1}\dots\tau_{1}A.$ (80)

We then approximate the evolution by $\tau_{1}$ by a truncated operator
$A_{1}$ such that

$\displaystyle\left\|\tau_{1}A-A_{1}\right\|=\varepsilon_{1},$ (81)

where $A_{1}$ is supported on sites at most a distance $\ell$ from the support
of $A$. We repeat the above approximation for the other time intervals to get

$\displaystyle\left\|\tau_{2}A_{1}-A_{2}\right\|=\varepsilon_{2},$ (82)
$\displaystyle\left\|\tau_{3}A_{2}-A_{3}\right\|=\varepsilon_{3},$ (83)
$\displaystyle\dots$
$\displaystyle\left\|\tau_{M}A_{M-1}-A_{M}\right\|=\varepsilon_{M}.$ (84)

At the end of this process, we have approximated $\tau(t)A$ by an operator
$A_{M}$ supported on sites located a distance of $M\ell$ from the support of
$A$. We bound the error of this approximation using the triangle inequality:

$\displaystyle\left\|\tau_{M}\dots\tau_{1}A-A_{M}\right\|\leq\varepsilon_{1}+\dots+\varepsilon_{M}.$
(85)

By choosing $M\ell$ slightly less than $r$, we guarantee that the support of
$A_{M}$ does not overlap with $X$, which implies that $K_{Y}(A_{M})=0$ and
therefore that the commutator

$\mathcal{C}(r,t)=\left\|K_{Y}(\tau A)\right\|\leq\left\|K_{Y}(\tau
A-A_{M})\right\|+\left\|K_{Y}(A_{M})\right\|=\left\|K_{Y}(\tau
A-A_{M})\right\|$ (86)

is at most the error of the approximation:
$\varepsilon\equiv\varepsilon_{1}+\dots+\varepsilon_{M}$. To find a bound on
$\varepsilon_{1}$, we trace out the part of $\tau_{1}A$ that lies outside of
$\mathcal{A}_{\ell}(Y)$, the ball of radius $\ell$ around the support of $A$:

$\displaystyle
A_{1}\equiv\frac{1}{\Tr(\mathbb{I}_{\mathcal{A}_{\ell}(Y)^{c}})}\Tr_{\mathcal{A}_{\ell}(Y)^{c}}(\tau_{1}A)\otimes\mathbb{I}_{\mathcal{A}_{\ell}(Y)^{c}}=\int_{\mathcal{A}_{\ell}(Y)^{c}}d\mu(W)W(\tau_{1}A)W^{\dagger},$
(87)

where $S^{c}$ denotes the complement of the set $S$ and the trace is rewritten
as an integral over Haar unitaries $W$ supported on
$\mathcal{A}_{\ell}(Y)^{c}$, and $\mu(W)$ denotes the Haar measure.

Now the error from approximating $\tau_{1}A$ with $A_{1}$ is given by

$\displaystyle\varepsilon_{1}=\left\|\tau_{1}A-A_{1}\right\|$
$\displaystyle=\left\|\tau_{1}A-\int_{\mathcal{A}_{\ell}(Y)^{c}}d\mu(W)W(\tau_{1}A)W^{\dagger}\right\|$
(88)
$\displaystyle=\left\|\int_{\mathcal{A}_{\ell}(Y)^{c}}d\mu(W)\left[\tau_{1}A-W(\tau_{1}A)W^{\dagger}\right]\right\|$
(89)
$\displaystyle\leq\int_{\mathcal{A}_{\ell}(Y)^{c}}d\mu(W)\left\|\left[\tau_{1}A,W\right]\right\|.$
(90)

Plugging this into Eq. 90 gives

$\displaystyle\varepsilon_{1}$
$\displaystyle=\left\|\tau_{1}A-A_{1}\right\|\leq\int_{\mathcal{A}_{\ell}(Y)^{c}}d\mu(W)\left\|A\right\|\phi(Y)f(\ell,\Delta
t)=\absolutevalue{A}\phi(Y)f(\ell,\Delta t),$ (91)

where $\Delta t=t/M$ is the size of each timestep. Applying this to all of the
errors yields

$\varepsilon_{j}\leq\absolutevalue{A}\phi(X_{j})f(\ell,\Delta t),$ (92)

where $X_{j}$ is the support of $A_{j}$. Thus the new bound is

$\displaystyle\mathcal{C}(r,t)\leq 2\|K_{Y}\|_{\infty}\varepsilon$
$\displaystyle\leq
2M\|K_{Y}\|_{\infty}\absolutevalue{A}\phi_{\text{max}}f(\ell,\Delta t)$ (93)
$\displaystyle=2\|K_{Y}\|_{\infty}\absolutevalue{A}\frac{t}{\Delta
t}\phi_{\text{max}}f(\ell,\Delta t),$ (94)

where $\phi_{\text{max}}=\max_{j}\phi(X_{j})$, and we replaced $M$ with
$t/\Delta t$. Without loss of generality, we may set $\Delta t=1$. Using the
form of $f(r,t)$ given in Eq. 78, this yields the bound

$\displaystyle\mathcal{C}(r,t)$ $\displaystyle\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|t\phi_{\max}\frac{e^{v}}{\left(\frac{r}{t}\right)^{\alpha-d-1}}$
(95) $\displaystyle\leq
C\|K_{Y}\|_{\infty}\left\|A\right\|\frac{t^{\alpha-d}}{r^{\alpha-2d}},$ (96)

which matches Eq. 6 in the main text.

##  Appendix F Generalization of the Chen & Lucas bound to open long-range
systems

In this section, we provide the proof of the bound in Eq. 8, which generalizes
the closed-system Lieb-Robinson bound from [13] to open systems. In the
process, we improve the tail of the bound from $1/r$ to $1/r^{\alpha-2-o(1)}$.
Our goal is to prove that, for an operator $A\in\mathcal{B}(X)$ supported on
$X$, for $K_{Y}\in\mathbb{L}_{Y}$ a superoperator supported on $Y$, and for
backward time-evolution operator $e^{\mathcal{L}^{\dagger}t}$, we have

$\displaystyle\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
C\left\|K_{Y}\right\|_{\infty}\left\|A\right\|\frac{t}{r^{\alpha-2}}.$ (97)

To do that, we use a trivial bound

$\displaystyle\left\|K_{Y}(e^{\mathcal{L}^{\dagger}t}A)\right\|\leq
2\left\|K_{Y}\right\|_{\infty}\left\|\mathbb{P}_{Y}e^{\mathcal{L}^{\dagger}t}A\right\|,$
(98)

where $\mathbb{P}_{Y}$ is the projector onto operators supported on sites at
distance $Y$ and beyond. We will now represent the operator $A$ by its
vectorized form $\left|A\right)$, so that $\mathbb{P}_{Y}$ acting on $A$ can
be viewed as a superoperator acting on the vectorized operator:
$\mathbb{P}_{Y}(A)=\mathbb{P}_{Y}\left|A\right)$. Also, from here on out, we
will represent $\mathcal{L}^{\dagger}$ by $\mathcal{L}$ for notational
convenience.

The quantity that we wish to bound is
$\left\|\mathbb{P}_{Y}e^{\mathcal{L}t}\left|A\right)\right\|$, which can be
expanded in a series

$\displaystyle\left\|\mathbb{P}_{Y}e^{\mathcal{L}t}\left|A\right)\right\|$
$\displaystyle=\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\mathcal{L}^{n}\left|A\right)=\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\sum_{\beta_{1},\beta_{2},\dots,\beta_{n}}\mathcal{L}_{\beta_{n}}\dots\mathcal{L}_{\beta_{2}}\mathcal{L}_{\beta_{1}}\left|A\right),$
(99)

where the $\beta_{i}$ correspond either to single-site terms or two-body
couplings, which we will refer to as “jumps.”

###  F.1 More definitions

We need a few more definitions before we can proceed. Consider a sequence of
jumps $\bm{\beta}=(\beta_{n},\dots,\beta_{1})$. First, we denote by
$\nu(\bm{\beta})$ the number of jumps in $\bm{\beta}$ and
$\nu_{q}(\bm{\beta})$ the number of order-$q$ jumps in $\bm{\beta}$. By
“order-$q$” jumps, we mean jumps that are of length at least $2^{q-1}$ and
less than $2^{q}$. For example, $\nu_{1}(\bm{\beta})$ is the number of
nearest-neighbor jumps in $\bm{\beta}$. $\nu_{2}(\bm{\beta})$ counts the
number of jumps of length $2,3$. Given a jump $\beta$,
$\textrm{dist}(\beta,y)$ is the minimum distance from the support of $\beta$
to $y$. The distance between a sequence of jumps $\bm{\beta}$ to $y$ is the
minimum distance between each jump and $y$. We also define a number $N_{q}$
for each $q$ as follows:

$\displaystyle N_{q}=\lceil\frac{\mu}{2^{q\gamma}}\frac{r}{2^{q}}\rceil,$
(100)

where $\gamma\in(0,1)$ is a parameter to be chosen later, and where $\mu<2$ is
a constant chosen to be small enough that

$\displaystyle\sum_{q=1}^{\infty}(N_{q}-1)2^{q}\leq\mu
r\sum_{q=1}^{\infty}2^{-q\gamma}<r.$ (101)

We list the other definitions below (see Fig. 3 for a diagram):

  * •

Given a sequence of jumps $\bm{\beta}$, we define its $q$-forward subsequence
according to Definition 1.

######  Definition 1.

Given a sequence of jumps $\bm{\beta}=(\beta_{n},\dots,\beta_{1})$, its
$q$-forward subsequence $\bm{\lambda}^{(q)}$ is constructed as followed:

    * –

Set $\bm{\lambda}^{(q)}=\\{\\}$ to be an empty sequence and define
$\textrm{dist}(\\{\\},y)=\textrm{dist}(x,y)$.

    * –

For $j=1,\dots,m$:

      * *

If $\textrm{dist}(\beta_{j},y)<\textrm{dist}(\bm{\lambda}^{(q)},y)$ and
$\beta_{j}$ is an order-$q$ jump, add $\beta_{j}$ to $\bm{\lambda}^{(q)}$.

We denote by $\mathcal{F}$ the map from $\bm{\beta}$ to its set of $q$-forward
subsequences $\Lambda=\\{\bm{\lambda}^{(q)}:q=1,\dots,r\\}$. This map is many-
to-one.

  * •

If the $q$-forward subsequence $\bm{\lambda}^{(q)}$ has at least $N_{q}$
jumps, we construct the irreducible $q$-forward subsequence
$\bm{\lambda}^{\prime(q)}$ by taking exactly the first $N_{q}$ jumps in
$\bm{\lambda}^{(q)}$. Otherwise, we say that there is no irreducible
$q$-forward subsequences.

  * •

We denote the map from $\Lambda=\\{\bm{\lambda}^{(q)}\\}$ to the set of
irreducible $q$-forward subsequences
$\Lambda^{\prime}=\\{\bm{\lambda}^{\prime(q)}\\}$ by $\mathcal{T}$. Note that
$\absolutevalue{\Lambda^{\prime}}$ can be less than $\absolutevalue{\Lambda}$
because the length of $\bm{\lambda}^{(q)}$ may be less than $N_{q}$ for some
$q$.

  * •

From a set
$\Lambda^{\prime}=\\{\bm{\lambda}^{\prime(q_{1})},\dots,\bm{\lambda}^{\prime(q_{k})}\\}$
of irreducible $q$-forward subsequences, we define
$\mathcal{I}(\Lambda^{\prime})=\\{\bm{\beta}:\mathcal{T}(\mathcal{F}(\bm{\beta}))\supseteq\Lambda^{\prime}\\}$
to be the set of sequences $\bm{\beta}$ that has $\Lambda^{\prime}$ in its set
of irreducible $q$-forward subsequences.

A sequence of jumps $\bm{\beta}$ $\Lambda=\\{\bm{\lambda}^{(q)}\\}$
$q$-forward subsequences $\Lambda^{\prime}=\\{\bm{\lambda}^{\prime(q)}\\}$
irreducible $q$-forward subsequences $\mathcal{F}$ $\mathcal{T}$ $\mathcal{I}$
Figure 3: A summary of the definitions regarding sequences and subsequences.

###  F.2 Proof

Lemma 3 below guarantees that, for each sequence $\bm{\beta}$ that contributes
to Eq. 99, there exists at least one irreducible $q$-forward subsequence
$\bm{\lambda}^{\prime(q)}$ for some $q$.

######  Lemma 3.

For each sequence $\bm{\beta}$, if
$\mathbb{P}_{Y}\mathcal{L}_{\bm{\beta}}\left|A\right)\neq 0$, then there
exists at least one $q$-forward subsequence such that
$\nu_{q}(\bm{\lambda}^{(q)})\geq N_{q}$.

The proof of this lemma is straightforward. If there exists no such $q$, then
$\nu_{\ell}(\bm{\lambda}^{(q)})\leq N_{q}-1$ for all $q$. By the construction
of $\bm{\lambda}$:

$\displaystyle
r\leq\sum_{q=1}^{r}\nu_{q}(\bm{\lambda}^{(q)})2^{q}\leq\sum_{q=1}^{r}(N_{q}-1)2^{q}<r,$
(102)

which is a contradiction.

In the following, we use the notation $\chi_{q}$ to denote whether
$\bm{\beta}$ has an irreducible $q$-forward subsequence:

$\displaystyle\chi_{q}\mathcal{L}_{\bm{\beta}}\left|A\right)=\begin{cases}\mathcal{L}_{\bm{\beta}}\left|A\right)&\text{if
}\exists\bm{\lambda}^{\prime(q)}\in\mathcal{T}(\mathcal{F}(\bm{\beta}))),\\\
0&\text{otherwise.}\end{cases}$ (103)

We can rewrite the series expansion of Eq. 99 as

$\displaystyle\mathbb{P}_{Y}e^{\mathcal{L}t}\left|A\right)$
$\displaystyle=\mathbb{P}_{Y}\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\sum_{\bm{\beta}}\mathcal{L}_{\bm{\beta}}\left|A\right)$
(104)
$\displaystyle=\mathbb{P}_{Y}\left[1-\prod_{q=1}^{\infty}(1-\chi_{q})\right]\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\sum_{\bm{\beta}}\mathcal{L}_{\bm{\beta}}\left|A\right),$
(105)

where Lemma 3 ensures that $1-\prod_{\ell}(1-\chi_{\ell})=1$ for all sequences
that contribute to Eq. 99. Expanding the product over $\ell$, we will get
terms of the form

$\displaystyle\mathcal{S}(q_{1},\dots,q_{k})$
$\displaystyle=(-1)^{k+1}\mathbb{P}_{Y}\chi_{q_{1}}\chi_{q_{2}}\dots\chi_{q_{k}}\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\sum_{\bm{\beta}}\mathcal{L}_{\bm{\beta}}\left|A\right)$
(106)
$\displaystyle=(-1)^{k+1}\mathbb{P}_{Y}\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\sum_{\bm{\lambda}^{\prime(q_{1})}}\dots\sum_{\bm{\lambda}^{\prime(q_{k})}}\sum_{\begin{subarray}{c}\bm{\beta}\in\mathcal{I}(\\{\bm{\lambda}^{\prime(q_{1})},\dots,\bm{\lambda}^{\prime(q_{k})}\\})\\\
\text{length}(\bm{\beta})=n\end{subarray}}\mathcal{L}_{\bm{\beta}}\left|A\right),$
(107)

for some distinct integers $q_{1},\dots,q_{k}$. In the last line, we sum over
all possible irreducible $q$-forward subsequences $\bm{\lambda}^{(q)}$, for
$q=q_{1},\dots,q_{k}$, then sum over all sequences $\bm{\beta}$ which contains
$\\{\bm{\lambda}^{\prime(q_{1})},\dots,\bm{\lambda}^{\prime(q_{k})}\\}$ in its
set of irreducible $q$-forward subsequences.

We will now upper-bound $\left\|\mathcal{S}(q_{1},\dots,q_{k})\right\|$.
First, let $\bm{\lambda}^{\prime}$ be a sequence consisting of all jumps in
$\bm{\lambda}^{\prime(q_{1})},\dots,\bm{\lambda}^{\prime(q_{k})}$ such that
the set of irreducible $\ell$-forward subsequences of $\bm{\lambda}^{\prime}$
is exactly
$\\{\bm{\lambda}^{\prime(q_{1})},\dots,\bm{\lambda}^{\prime(q_{k})}\\}$. From
$\bm{\lambda}^{\prime}$, we construct $\bm{\beta}$:

$\displaystyle\bm{\beta}=\left(\beta_{m+1,j_{m+1}},\dots,\beta_{m+1,1}\lambda^{\prime}_{m},\dots,\lambda_{2}\beta_{2,j_{2}},\dots,\beta_{2,1},\lambda^{\prime}_{1},\beta_{1,j_{1}},\dots,\beta_{1,1}\right),$
(108)

where
$(\lambda^{\prime}_{m},\dots,\lambda^{\prime}_{1})=\bm{\lambda}^{\prime}$,
$j_{1},\dots,j_{m+1}$ are nonnegative integers, $\beta_{i,j}\in\Gamma_{i}$,
and the sets $\Gamma_{i}$ are constructed recursively for $i=1,\dots,m+1$ as
follows:

  * •

$\Gamma_{1}=\\{(x^{\prime},y^{\prime}):\textrm{dist}((x^{\prime},y^{\prime}),y)<\textrm{dist}(x,y)\text{
if }(x^{\prime},y^{\prime})\text{ is an order-$q$ jump, where
$q=q_{1},\dots,q_{k}$}\\}.$

  * •

Set $c_{q}=r$ for all $q=q_{1},\dots,q_{k}$. Each $c_{q}$ will remember the
distance from $y$ to the last length-$q$ jump. For the sake of the proof, let
$c_{q}=\infty$ for all other $q$.

  * •

For $i=2$ to $m$:

    * –

$\Gamma_{i}=\\{(x^{\prime},y^{\prime}):\textrm{dist}((x^{\prime},y^{\prime}),y)<c_{q(x^{\prime},y^{\prime})}.$

    * –

Update $c_{q}=\textrm{dist}(\lambda^{\prime}_{i},y)$, where $q$ is the order
of the jump $\lambda^{\prime}_{i}$.

  * •

$\Gamma_{m+1}=\\{(x^{\prime},y^{\prime})\\}$ is the set of all possible jumps.

The point of this construction is that each sequence $\bm{\beta}$ appears
exactly once. We can then rewrite

$\displaystyle\mathcal{S}(q_{1},\dots,q_{k})$
$\displaystyle=(-1)^{k+1}\mathbb{P}_{Y}\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\sum_{\bm{\lambda}^{\prime(q_{1})}}\dots\sum_{\bm{\lambda}^{\prime(q_{k})}}\sum_{\begin{subarray}{c}\bm{\beta}\in\mathcal{I}(\\{\bm{\lambda}^{\prime(q_{1})},\dots,\bm{\lambda}^{\prime(q_{k})}\\})\\\
\text{length}(\bm{\beta})=n\end{subarray}}\mathcal{L}_{\bm{\beta}}\left|A\right),$
(109)
$\displaystyle=(-1)^{k+1}\mathbb{P}_{Y}\sum_{\bm{\lambda}^{\prime(q_{1})}}\dots\sum_{\bm{\lambda}^{\prime(q_{k})}}\sum_{\bm{\lambda}^{\prime}}\sum_{j_{m+1}=0}^{\infty}\dots\sum_{j_{1}=0}^{\infty}\frac{t^{m+\sum_{l=1}^{m}{j_{l}}}}{(m+\sum_{l=1}^{m}{j_{l}})!}\mathcal{L}_{\Gamma_{m+1}}^{j_{m+1}}\mathcal{L}_{\lambda_{m+1}}\dots\mathcal{L}_{\lambda_{1}}\mathcal{L}_{\Gamma_{1}}^{j_{1}}\left|A\right),$
(110)
$\displaystyle=(-1)^{k+1}\mathbb{P}_{Y}\sum_{\bm{\lambda}^{\prime(q_{1})}}\dots\sum_{\bm{\lambda}^{\prime(q_{k})}}\sum_{\bm{\lambda}^{\prime}}\int_{\Delta^{m}(t)}dt_{1}\dots
dt_{m}e^{\mathcal{L}_{\Gamma_{m+1}}^{j_{m+1}}(t-t_{m})}\mathcal{L}_{\lambda_{m+1}}\dots\mathcal{L}_{\lambda_{1}}e^{\mathcal{L}_{\Gamma_{1}}^{j_{1}}t_{1}}\left|A\right),$
(111)

where $\Delta^{m}(t)$ is the simplex defined by $0\leq t_{1}\leq\dots\leq
t_{m}\leq t$. Now, we use the triangle inequality:

$\displaystyle\left\|\mathcal{S}(q_{1},\dots,q_{k})\right\|$
$\displaystyle\leq\frac{3}{2}\sum_{\bm{\lambda}^{\prime(q_{1})}}\dots\sum_{\bm{\lambda}^{\prime(q_{k})}}\sum_{\bm{\lambda}^{\prime}}\frac{t^{m}}{m!}\frac{1}{q_{1}^{\alpha
N_{q_{1}}}}\dots\frac{1}{q_{k}^{\alpha N_{q_{k}}}}$ (112)
$\displaystyle\leq\frac{3}{2}\binom{r2^{q}}{N_{q_{1}}}\dots\binom{r2^{q}}{N_{q_{k}}}\binom{m}{N_{q_{1}},\dots,N_{q_{k}}}\frac{t^{m}}{m!}\frac{1}{2^{\alpha
q_{1}N_{q_{1}}}}\dots\frac{1}{2^{\alpha q_{k}N_{q_{k}}}}$ (113)
$\displaystyle=\frac{3}{2}\prod_{i=1,\dots,k}\left[\binom{r2^{q}_{i}}{N_{q_{i}}}\frac{t^{N_{q_{i}}}}{N_{q_{i}}!}\frac{1}{2^{\alpha
q_{i}N_{q_{i}}}}\right],$ (114)

where in the last two lines we use the fact that
$m=N_{q_{1}}+\dots+N_{q_{k}}$. Plugging this bound into Eq. 105, we have

$\displaystyle\left\|\mathbb{P}_{Y}e^{\mathcal{L}t}\left|A\right)\right\|\leq-1+\prod_{q}\left[1+\frac{3}{2}\binom{r2^{q}}{N_{q}}\frac{t^{N_{q}}}{N_{q}!}\frac{1}{q^{\alpha
N_{q}}}\right].$ (115)

Now we use $1+x\leq e^{x}$ to bound

$\displaystyle\left\|\mathbb{P}_{Y}e^{\mathcal{L}t}\left|A\right)\right\|\leq-1+\exp\left[\frac{3}{2}\sum_{q}\binom{r2^{q}}{N_{q}}\frac{t^{N_{q}}}{N_{q}!}\frac{1}{q^{\alpha
N_{q}}}\right].$ (116)

Let $q_{*}$ be the largest integer such that $2^{q_{*}(\gamma+1)}\leq(\mu
r)^{1-\gamma}$. Note that $\mu r/2^{q(\gamma+1)}>1$ for all $q\leq q_{*}$. We
divide the sum in Eq. 116 into two parts:

$\displaystyle\sum_{q}\binom{r2^{q}}{N_{q}}\frac{t^{N_{q}}}{N_{q}!}\frac{1}{2^{\alpha
qN_{q}}}$
$\displaystyle\leq\underbrace{\sum_{q=1}^{q_{*}-1}\binom{r2^{q}}{N_{q}}\frac{t^{N_{q}}}{N_{q}!}\frac{1}{2^{\alpha
qN_{q}}}}_{=S_{1}}+\underbrace{\sum_{q=q_{*}}^{r}\frac{rt}{2^{(\alpha-1)q}}}_{=S_{2}}.$
(117)

First, we estimate $S_{2}$:

$\displaystyle
S_{2}\leq\frac{1}{1-2^{-\alpha}}\frac{rt}{2^{q_{*}(\alpha-1)}}\leq\underbrace{\frac{1}{1-2^{-\alpha}}\mu^{(1-\alpha)/(\gamma+1)}}_{=c_{3}}\frac{t}{r^{\frac{\alpha-1}{\gamma+1}-1}}.$
(118)

Next, we estimate $S_{1}$. Note that $N_{q}\geq\frac{\mu r}{2^{q(\gamma+1)}}$
for all $q$:

$\displaystyle S_{1}$
$\displaystyle\leq\sum_{q=1}^{q_{*}-1}\left(\frac{e^{2}rt}{N^{2}_{q}2^{(\alpha-1)q}}\right)^{N_{q}}$
(119)
$\displaystyle\leq\sum_{q=1}^{q_{*}}\left(\frac{e^{2}t}{\mu^{2}r}2^{q(2\gamma+3-\alpha)}\right)^{N_{q}},$
(120)

where we have used the Stirling’s approximation $x!>x^{x}e^{-x}$. When
$q\rightarrow 1$, $N_{q}\propto r$. The corresponding term in $S_{1}$ decays
with $r$ at least exponentially as $(t/r)^{r}$. On the other hand, when
$q\rightarrow q_{*}$, $N_{q}\rightarrow 1$ and the corresponding term in
$S_{1}$ is instead suppressed by $2^{q(2\gamma+3-\alpha)}$ for all
$\alpha>3+2\gamma$. This limit analysis suggests that we should use two
different bounds on $S_{1}$ for small $q$ and large $q$. For that, we define

$\displaystyle q_{0}\equiv\lfloor\frac{1}{1+\gamma}\log_{2}(\mu
r^{\kappa})\rfloor\leq\frac{1}{1+\gamma}\log_{2}(\mu r^{\kappa})$ (121)

and divide up $S_{1}$ into two sums over $q\leq q_{0}$ and $q_{0}<q\leq
q_{*}$:

$\displaystyle
S_{1}\leq\underbrace{\sum_{q=1}^{q_{0}-1}\left(\frac{e^{2}t}{\mu^{2}r}2^{q(2\gamma+3-\alpha)}\right)^{N_{q}}}_{=S_{1a}}+\underbrace{\sum_{q=q_{0}}^{q_{*}}\left(\frac{e^{2}t}{\mu^{2}r}2^{q(2\gamma+3-\alpha)}\right)^{N_{q}}}_{=S_{1b}}.$
(122)

First, we take the sum over $q\leq q_{0}$. We assume that $\alpha>2\gamma+3$
and $t\leq\mu^{2}r/e^{2}$, so that the inner summand satisfies

$\displaystyle\left(\frac{e^{2}t}{\mu^{2}r}2^{q(2\gamma+3-\alpha)}\right)\leq
1$ (123)

for all $q\leq q_{0}$. Because $N_{q}$ decreases with $q$, we upper bound

$\displaystyle
S_{1a}=\sum_{q=1}^{q_{0}}\left(\frac{e^{2}t}{\mu^{2}r}2^{q(2\gamma+3-\alpha)}\right)^{N_{q}}$
$\displaystyle\leq\left(\frac{e^{2}t}{\mu^{2}r}\right)^{N_{q_{0}}}\sum_{q=1}^{q_{0}}2^{q(2\gamma+3-\alpha)N_{q}}$
(124) $\displaystyle\lesssim\left(\frac{e^{2}t}{\mu^{2}r}\right)^{\frac{\mu
r}{2^{q_{0}(\gamma+1)}}}$ (125)
$\displaystyle\leq\left(\frac{e^{2}t}{\mu^{2}r}\right)^{r^{1-\kappa}}$ (126)
$\displaystyle\leq\frac{t}{r}e^{-r^{1-\kappa}},$ (127)

where in the last line we further assume $t\leq\mu^{2}r/e^{2}$. This gives the
sum over $q\leq q_{0}$ in the term $S_{1}$. To bound the sum over $q_{0}<q\leq
q_{*}$, we note that $N_{q-1}\geq N_{q}+1$ for all $q<q_{*}$. To prove this,
suppose $N_{q-1}=N_{q}$. That means

$\displaystyle\frac{\mu r}{2^{(q-1)(\gamma+1)}}<N_{q-1}=N_{q}\leq\frac{\mu
r}{2^{q(\gamma+1)}}+1$ (128) $\displaystyle\Leftrightarrow$ $\displaystyle
1>(2^{\gamma+1}-1)\frac{\mu r}{2^{q(\gamma+1)}}>\frac{\mu
r}{2^{q(\gamma+1)}},$ (129)

which contradicts with $\mu r/2^{q(\gamma+1)}>1$ for all $q<q_{*}$. Therefore,
$N_{q-1}\geq N_{q}+1$ for all $q<q_{*}$. Since $N_{q_{*}}=1$, it follows that
$N_{q_{*}-n}\geq n+1>n$ for all $n\geq 1$. We make the substitution
$n=q_{*}-q$ to obtain

$\displaystyle S_{1b}$
$\displaystyle=\sum_{q=q_{0}}^{q_{*}}\left(\frac{e^{2}t}{\mu^{2}r}2^{q(2\gamma+3-\alpha)}\right)^{N_{q}}\leq\sum_{n=1}^{q_{*}-q_{0}}\left(\frac{e^{2}t}{\mu^{2}r}2^{(q_{*}-n)(2\gamma+3-\alpha)}\right)^{n},$
(130)

again assuming that $\alpha>3+2\gamma$ and $e^{2}t/(\mu^{2}r)<1$. Now, using
the fact that $q_{*}-n\geq q_{0}$, we have

$\displaystyle 2^{(q_{*}-n)(2\gamma+3-\alpha)}\leq
2^{q_{0}(2\gamma+3-\alpha)}\leq r^{\kappa(2\gamma+3-\alpha)}.$ (132)

Plugging this into the sum yields

$\displaystyle\sum_{n=1}^{q_{*}-q_{0}}\left(\frac{e^{2}t}{\mu^{2}r}2^{(q_{*}-n)(2\gamma+3-\alpha)}\right)^{n}$
$\displaystyle\leq\sum_{n=1}^{q_{*}-q_{0}}\left(\frac{e^{2}t}{\mu^{2}r}r^{\kappa(2\gamma+3-\alpha)}\right)^{n}$
(133)
$\displaystyle=\frac{e^{2}t}{\mu^{2}r^{1-\kappa(2\gamma+3-\alpha)}}\sum_{n=0}^{q_{*}-q_{0}-1}\left(\frac{e^{2}t}{\mu^{2}r^{1-\kappa(2\gamma+3-\alpha)}}\right)^{n}$
(134)
$\displaystyle\leq\frac{e^{2}t}{\mu^{2}r^{1-\kappa(2\gamma+3-\alpha)}}\frac{1}{1-\frac{e^{2}t}{\mu^{2}r^{1-\kappa(2\gamma+3-\alpha)}}}$
(135)
$\displaystyle\leq\underbrace{2\frac{e^{2}}{\mu^{2}}}_{=c_{2}}\frac{t}{r^{1-\kappa(2\gamma+3-\alpha)}},$
(136)

assuming that
$\frac{e^{2}t}{\mu^{2}r^{1-\kappa(2\gamma+3-\alpha)}}\leq\frac{1}{2}$.
Combining everything, we have

$\displaystyle S_{1}+S_{2}\leq
c_{1}\left(\frac{t}{r}e^{-r^{1-\kappa}}\right)+c_{2}\frac{t}{r^{1-\kappa(2\gamma+3-\alpha)}}+c_{3}\frac{t}{r^{\frac{\alpha-1}{1+\gamma}-1}}.$
(137)

We make the simplification that $\kappa=1-\gamma$, so that

$\displaystyle
1-\kappa(2\gamma+3-\alpha)=1-(1-\gamma)(\alpha-3-2\gamma)=\alpha-2\underbrace{-2\gamma-\gamma\alpha+3\gamma+2\gamma^{2}}_{=o(1)}.$
(138)

In addition, for all $\gamma>0$, there exists a constant $c_{\gamma}$ that may
depend on $\alpha$ such that

$\displaystyle e^{-r^{\gamma}}\leq c_{\gamma}\frac{1}{r^{\alpha-3}}$ (139)

for all $r>0$. Therefore,

$\displaystyle\frac{t}{r}e^{-r^{\gamma}}\leq
c_{\gamma}\frac{t}{r^{\alpha-2}}.$ (140)

Substituting Eqs. 140 and 138 into Eq. 137 and letting
$c=c_{1}c_{\gamma}+c_{2}+c_{3}$, we have the desired bound:

$\displaystyle\left\|\mathbb{P}_{Y}e^{\mathcal{L}t}\left|A\right)\right\|\leq
c\frac{t}{r^{\alpha-2-o(1)}},$ (141)

which is exactly Eq. 8 in the main text.