Classical simulation of quantum circuits by dynamical localization: analytic results for Pauli-observable scrambling in time-dependent disorder

TL;DR

This paper demonstrates that certain quantum circuits exhibit dynamical localization, allowing classical simulation of quantum information scrambling through analytic formulas for out-of-time-ordered correlators in disordered spin systems.

Contribution

It extends Anderson localization to quantum circuits, providing an analytic approach to simulate quantum dynamics efficiently in localized phases.

Findings

Analytic formula for out-of-time-ordered correlator in matchgate circuits

Efficient classical evaluation of local observable dynamics

Identification of localized phases under disordered Hamiltonians

Abstract

We extend the concept of Anderson localization, the confinement of quantum information in a spatially irregular potential, to quantum circuits. Considering matchgate circuits, generated by time-dependent spin-1/2 XY Hamiltonians, we give an analytic formula for the out-of-time-ordered correlator of a local observable, and show that it can be efficiently evaluated by a classical computer even when the explicit Heisenberg time evolution cannot. Because this quantity bounds the average error incurred by truncating the evolution to a spatially limited region, we demonstrate dynamical localization as a means for classically simulating quantum computation and give examples of localized phases under certain spatio-temporal disordered Hamiltonians.

Tables

Table 1. Table 1: Summary of the time-dependent models considered in this work, corresponding to the numerical phase diagram shown in Fig.  2 . ν j ∼ [ − ν , ν ] similar-to subscript 𝜈 𝑗 𝜈 𝜈 \nu_{j}\sim[-\nu,\nu] are chosen uniformly randomly, and the κ j ​ ( t ) ∼ [ − Δ , Δ ] similar-to subscript 𝜅 𝑗 𝑡 Δ Δ \kappa_{j}(t)\sim[-\Delta,\Delta] are uniformly randomly sampled every δ ​ t = 0.25 𝛿 𝑡 0.25 \delta t=0.25

Parent Hamiltonian (8)
$H\left(t\right)={\sum }_{j=1}^{n-1}2{\mu }_j\left(t\right)\left(X_j{X}_{j+1}+Y_j{Y}_{j+1}\right)+{\sum }_{j=1}^{n}2{\nu }_j\left(t\right)Z_j$
Model	Fluctuation	Fixed	Varying
1	${\nu }_j\left(t\right)={\nu }_j+\left(2/\delta t\right){\kappa }_j\left(t\right)$	$\mu =1$	$\nu$, $\Delta$
2	${\mu }_j\left(t\right)=\mu +\left(2/\delta t\right){\kappa }_j\left(t\right)$	$\nu =1$	$\mu$, $\Delta$

$\overrightarrow{\alpha }$	$\equiv$	$({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k)$, where , for all ${\alpha }_j\in \{1,2,\dots ,N\}$ for some $N\in \mathbb{N}$, and $j\in \{1,2,\dots ,k-1\}$; i.e. an ordered tuple of indices.
$\left|\overrightarrow{\alpha }\right|$	$\equiv$	$k$, the number of elements in the tuple $\overrightarrow{\alpha }=({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k)$.
$\left[k\right]$	$\equiv$	$(1,2,\dots ,k)$, the particular tuple consisting of consecutive integers from $1$ to $k$.
$\overline{\alpha }$	$\equiv$	Ordered tuple of indices in the set-complement of $\overrightarrow{\alpha }$ in $\left[N\right]$.
$\left(\right)$	$\equiv$	The empty tuple, for which $\left|\left(\right)\right|=0$.
$\overrightarrow{\alpha }\cup \overrightarrow{\beta }$	$\equiv$	Ordered tuple of all indices in either $\overrightarrow{\alpha }$ or $\overrightarrow{\beta }$.
$\overrightarrow{\alpha }\cap \overrightarrow{\beta }$	$\equiv$	Ordered tuple of all indices in both $\overrightarrow{\alpha }$ and $\overrightarrow{\beta }$.
$\overrightarrow{\alpha }\subset \overrightarrow{\beta }$	$\equiv$	$\overrightarrow{\alpha }$ is an ordered tuple of indices ${\alpha }_j\in \overrightarrow{\beta }$ for all .
$\overrightarrow{\alpha }\backslash {\alpha }_j$	$\equiv$	Ordered tuple of indices from $\overrightarrow{\alpha }$ with the particular element ${\alpha }_j$ removed.
$(\overrightarrow{\alpha },\overrightarrow{\beta })$	$\equiv$	$({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k,{\beta }_1,{\beta }_2,\dots )$; i.e. the unordered concatenation of $\overrightarrow{\alpha }$ and $\overrightarrow{\beta }$.
$C_{\overrightarrow{\alpha }}$	$\equiv$	${\prod }_{j=1}^{\left|\overrightarrow{\alpha }\right|}{c}_{{\alpha }_j}$, the product of Majorana operators with indices in ascending order, called a Majorana configuration operator.
${𝐮}_{\overrightarrow{\alpha }\overrightarrow{\beta }}$	$\equiv$	Submatrix ${𝐮}_{\overrightarrow{\alpha }\overrightarrow{\beta }}$ of $𝐮$, with row indices from $\overrightarrow{\alpha }$ and column indices from $\overrightarrow{\beta }$; i.e. ${\left({𝐮}_{\overrightarrow{\alpha }\overrightarrow{\beta }}\right)}_{jk}=u_{{\alpha }_j{\beta }_k}$.

Equations

$$||[A,B(t)]||\lesssim \min (|t|,1)e^{-\eta d(A,B)}.$$

$$G(V,W)=\left(\begin{array}{cccc}{V}_{00}& 0& 0& V_{01}\\ 0& W_{00}& W_{01}& 0\\ 0& W_{10}& W_{11}& 0\\ V_{10}& 0& 0& V_{11}\end{array}\right),$$

$$c_{2k-1}\equiv {\otimes }_{j=1}^{k-1}{Z}_j\otimes X_k\text{mbox}{c}_{2k}\equiv {\otimes }_{j=1}^{k-1}{Z}_j\otimes Y_k,$$

$$U^{†}{c}_{\mu }U=\sum _{\nu =1}^{2n}{u}_{\mu \nu }{c}_{\nu },$$

$$U^{†}{C}_{\vec{\alpha }}U$$

$$Z_k=-i{C}_{\left(2k-1,2k\right)}\text{mbox}\mathrm{and}\text{mbox}{X}_k=(-i{)}^{k-1}{C}_{\left(1,\dots ,2k-1\right)}.$$

$$\frac{1}{2^{n+2}}||[A,U^{†}BU]|{|}^2=\frac{1}{2}\left\{1\pm \det [{\mathbf{u}}_{\vec{\alpha }[2n]}(\mathbf{I}-2{\mathbf{P}}_{\vec{\eta }}){\mathbf{u}}_{[2n]\vec{\alpha }}^{\mathrm{T}}]\right\},$$

$$\frac{1}{2^{n+2}}||[{\mathbf{n}}_s\cdot {\mathbf{\sigma }}_s,U^{†}BU]|{|}^2\equiv {\mathbf{n}}_s^{\ast }\cdot {\mathbf{M}}_s\cdot {\mathbf{n}}_s$$

$${\mathcal{E}}_s(O)=\frac{1}{4}\left(O+\sum _{k\in \{x,y,z\}}{\sigma }_s^{k}O{\sigma }_s^k\right),$$

$$\overline{|⟨U^{†}BU⟩-⟨({\otimes }_{s\in S}{\mathcal{E}}_s)[U^{†}BU]⟩|}\le \sum _{s\in S}\sqrt{{\mathbf{n}}_s^{\ast }\cdot {\mathbf{M}}_s\cdot {\mathbf{n}}_s},$$

$$H(t)=\sum _{j=1}^{n-1}2{\mu }_j(t)\left(X_j{X}_{j+1}+Y_j{Y}_{j+1}\right)+\sum _{j=1}^{n}2{\nu }_j(t)Z_j$$

$$U^{†}{C}_{\vec{\alpha }}U=\sum _{\{\vec{\beta }||\vec{\beta }|=|\vec{\alpha }|\}}\det \left({\mathbf{u}}_{\vec{\alpha }\vec{\beta }}\right)C_{\vec{\beta }},$$

$$U^{†}{c}_{{\alpha }_1}{c}_{{\alpha }_2}U$$

$$=\left(\sum _{{\beta }_1}{u}_{{\alpha }_1{\beta }_1}{c}_{{\beta }_1}\right)\left(\sum _{{\beta }_2}{u}_{{\alpha }_2{\beta }_2}{c}_{{\beta }_2}\right)$$

$$=\left(\sum _{\{({\beta }_1,{\beta }_2)|{\beta }_1={\beta }_2\}}+\sum _{\{({\beta }_1,{\beta }_2)|{\beta }_1<{\beta }_2\}}+\sum _{\{({\beta }_1,{\beta }_2)|{\beta }_1>{\beta }_2\}}\right)u_{{\alpha }_1{\beta }_1}{u}_{{\alpha }_2{\beta }_2}{c}_{{\beta }_1}{c}_{{\beta }_2}$$

$$=\left(\sum _{{\beta }_1}{u}_{{\alpha }_1{\beta }_1}{u}_{{\alpha }_2{\beta }_1}\right)I+\sum _{\{({\beta }_1,{\beta }_2)|{\beta }_1<{\beta }_2\}}\left(u_{{\alpha }_1{\beta }_1}{u}_{{\alpha }_2{\beta }_2}-u_{{\alpha }_1{\beta }_2}{u}_{{\alpha }_2{\beta }_1}\right)c_{{\beta }_1}{c}_{{\beta }_2}$$

$$=\sum _{\{({\beta }_1,{\beta }_2)|{\beta }_1<{\beta }_2\}}\det \left[{\mathbf{u}}_{\left({\alpha }_1,{\alpha }_2\right)\left({\beta }_1,{\beta }_2\right)}\right]c_{{\beta }_1}{c}_{{\beta }_2}$$

$$\{c_{\mu },c_{\nu }\}=2{\delta }_{\mu \nu }I,$$

$$U^{†}{C}_{\vec{\alpha }}{c}_{{\alpha }_{k+1}}U$$

$$=\left[\sum _{\{\vec{\beta }||\vec{\beta }|=|\vec{\alpha }|\}}\det \left({\mathbf{u}}_{\vec{\alpha }\vec{\beta }}\right)C_{\vec{\beta }}\right]\left(\sum _{{\beta }_{k+1}}{u}_{{\alpha }_{k+1}{\beta }_{k+1}}{c}_{{\beta }_{k+1}}\right)$$

$$U^{†}{C}_{\vec{\alpha }}{c}_{{\alpha }_{k+1}}U$$

$$\sum _{\gamma \notin \vec{\beta }\backslash {\beta }_{k+1}}{u}_{{\alpha }_{k+1}\gamma }\det \left[{\mathbf{u}}_{\vec{\alpha }\left(\vec{\beta }\backslash {\beta }_{k+1},\gamma \right)}\right]=\sum _{\gamma }{u}_{{\alpha }_{k+1}\gamma }\det \left[{\mathbf{u}}_{\vec{\alpha }\left(\vec{\beta }\backslash {\beta }_{k+1},\gamma \right)}\right],$$

$$\sum _{\gamma }{u}_{{\alpha }_{k+1}\gamma }\det \left[{\mathbf{u}}_{\vec{\alpha }\left(\vec{\beta }\backslash {\beta }_{k+1},\gamma \right)}\right]$$

$$\sum _{\gamma }{u}_{{\alpha }_{k+1}\gamma }{u}_{{\alpha }_l\gamma }={\delta }_{{\alpha }_{k+1}{\alpha }_l}$$

$$\sum _{\gamma \in \vec{\beta }\cup {\beta }_{k+1}}(-1{)}^{|\{j\le k+1|{\beta }_j>\gamma \}|}{u}_{{\alpha }_{k+1}\gamma }\det \left[{\mathbf{u}}_{\vec{\alpha },\left(\vec{\beta }\cup {\beta }_{k+1}\right)\backslash \gamma }\right]$$

$$(-1{)}^{|\{j\le k+1|{\beta }_j>\gamma \}|}=(-1{)}^{(k+1)-s}=(-1{)}^{(k+1)+s}$$

$$U^{†}{C}_{\vec{\alpha }}{c}_{{\alpha }_{k+1}}U$$

$$||[A,U^{†}BU]|{|}^2$$

$$=-\mathrm{tr}\{[A,U^{†}BU{]}^2\}$$

$$=-\mathrm{tr}\left[A\left(U^{†}BU\right)A\left(U^{†}BU\right)-A\left(U^{†}BU\right)\left(U^{†}BU\right)A-\left(U^{†}BU\right)AA\left(U^{†}BU\right)+\left(U^{†}BU\right)A\left(U^{†}BU\right)A\right]$$

$$=\mathrm{tr}\left[2I-2A\left(U^{†}BU\right)A\left(U^{†}BU\right)\right]$$

$$||[A,U^{†}BU]|{|}^2$$

$$\frac{1}{2^{n+2}}||[A,U^{†}BU]|{|}^2=\frac{1}{2}\left\{1-\frac{1}{2^n}\mathrm{tr}\left[A\left(U^{†}BU\right)A\left(U^{†}BU\right)\right]\right\}$$

$$A^{†}$$

$$i^a{C}_{\vec{\eta }}$$

Full text

*Introduction.—*A peculiar phenomenon exhibited uniquely by quantum lattice systems is the suppression of conductance in the presence of disorder. This effect, known as Anderson localization Anderson (1958); Lahini et al. (2010); Abrahams (2010); Crespi et al. (2013) in the single-particle setting and many-body localization Gornyi et al. (2005); Basko et al. (2006); Oganesyan and Huse (2007); Pal and Huse (2010); Vosk and Altman (2014); Bar Lev et al. (2015); Bera et al. (2015) in the interacting-multi-particle regime, is a result of interference, which confines a local disturbance to a bounded region near its initial position for a very long time. As a result, these systems do not act as thermal reservoirs for their own subsystems Deutsch (1991); Pal (2012); Nandkishore and Huse (2015); Huse et al. (2013); Kjäll et al. (2014), since local subsystems retain information about their initial conditions forever.

An important consequence of the fact that local quantum information does not mix, or scramble, among nonlocal degrees of freedom in localizing systems is that many properties of these systems can be efficiently simulated classically. Such properties include local integrals of motion Kim et al. (2014); You et al. (2016), Hamiltonians’ eigenbases Bauer and Nayak (2013); You et al. (2016); Friesdorf et al. (2015); Huang (2015); Yu et al. (2017); Pollmann et al. (2016), unitary time-evolution operators Burrell and Osborne (2007); Žnidarič et al. (2008), and samples from their output distributions Deshpande et al. (2017). In light of these many results, we ask the question of whether localization could manifest in time-dependent quantum systems, such as those performing a quantum computation. If this is possible, then it would allow for the efficient simulation of an otherwise apparently complex quantum algorithm by classical means.

However, there are currently very few prior investigations into localization in the time-dependent regime. Initial explorations into fluctuating disorder Burrell et al. (2009); Sacha and Delande (2016) and Floquet circuit ensembles Chandran and Laumann (2015) suggest that a form of localization persists in these time-dependent cases, yet few analytic results are known in general. In this work, we consider the setting of nearest-neighbor matchgate circuits, which are generated by time-dependent spin-1/2 XY Hamiltonians, and can be mapped onto the dynamics of free-fermions with arbitrarily time-dependent single-particle Hamiltonians by the Jordan-Wigner transformation Knill (2001); Terhal and DiVincenzo (2002); Bravyi and Kitaev (2002); Bravyi (2006); Jozsa and Miyake (2008); de Melo et al. (2013); Brod and Childs (2014). These circuits therefore constitute the natural framework in which to study the generalization of Anderson localization to quantum circuits.

Despite the encoding by free-fermion dynamics, some properties of these circuits are not known to be classically simulable. The example we consider is the Pauli-expectation value $\langle U^{\dagger}X_{j}U\rangle$ on an arbitrary qubit $j$ in the output of a matchgate circuit $U$ from an initial arbitrary product state. Despite being local in the qubit picture, this observable takes the form of a long-range correlation function in the fermion picture and requires exponential resources to simulate by brute-force. We solve this problem for circuits $U$ describing localizing dynamics by exploiting the confinement of their measurement observables in the Heisenberg picture. For time-independent Hamiltonians, this confinement is described formally by the so-called zero-velocity Lieb-Robinson bound Hamza et al. (2012)

[TABLE]

This inequality states that the degree of noncommutativity between the local observable $A$ and time-evolved observable $B(t)\equiv U^{\dagger}BU$ initially separated by lattice distance $d(A,B)>0$, such that $[A,B]=0$, is exponentially decaying with decay constant $\eta>0$. It gives an effective speed at which disturbances propagate Lieb and Robinson (1972); Nachtergaele and Sims (2006); Hastings and Koma (2006); Nachtergaele and Sims (2010), which goes to zero with increasing propagation time in localizing systems. Correlations between distant lattice sites take exponential time to develop Bravyi et al. (2006).

In the case where $A$ and $B$ are unitary, and the norm taken is the Frobenius norm $||O||^{2}\equiv\operatorname{tr}{\left(O^{\dagger}O\right)}$, the left-hand side of (1) is known as the infinite-temperature out-of-time-ordered correlation function (OTO correlator). This quantity has arisen as a useful diagnostic tool for studying scrambling in chaotic quantum chaotic systems Maldacena et al. (2016); Roberts and Yoshida (2017); Aleiner et al. (2016); Yunger Halpern (2017); Cotler et al. (2017), including black holes Shenker and Stanford (2014); Kitaev (2014, 2015); Roberts et al. (2015), and recently, for many-body localization Swingle and Chowdhury (2017); Chen (2016); Huang et al. (2017); He and Lu (2017); Fan et al. (2017). As a first result, we provide an analytic formula for this quantity when $A$ and $B$ are Pauli observables, and the time evolution is described by a matchgate circuit. This is surprising considering that the evolution itself cannot even be stored efficiently by a classical computer in general, and so it constitutes an exponential speedup over the brute-force method. We next show that this quantity bounds the average-case change in expectation-value magnitude from truncating the Heisenberg evolution of $B(t)$ to a subset of qubits and thus provides a measure of the expected error incurred by such truncation. Finally, we provide numerical analysis verifying the bound (1) for two natural models of time-dependent disorder and construct phase diagrams demonstrating their transitions to localizing dynamics and subsequent classical simulability.

*Background.—*Define a matchgate $G(V,W)$ to be the following 2-qubit unitary, written in the (ordered) computational basis $\{\left|00\right\rangle,\left|01\right\rangle,\left|10\right\rangle,\left|11\right\rangle\}$ as

[TABLE]

where $V,W\in\mathbf{SU}(2)$ are single-qubit unitaries (crucially, $\det V=\det W$). $G(V,W)$ preserves the eigenspaces of $Z\otimes Z$ and so may be written as $e^{i\mathcal{L}}$, where $\mathcal{L}$ is an element of the vector space spanned by $\{X\otimes X,X\otimes Y,Y\otimes X,\\ Y\otimes Y,Z\otimes I,I\otimes Z\}$. When the 2 qubits on which $G(V,W)$ acts are a nearest-neighboring pair, $\mathcal{L}$ is an instance of the 2-qubit spin-1/2 XY model. Such a Hamiltonian on $n$ qubits is quadratic in the $2n$ Majorana operators $\{c_{\mu}\}_{\mu=1}^{2n}$, given by the Jordan-Wigner transformation

[TABLE]

where $c_{\mu}^{2}=I$ and $c_{\mu}c_{\nu}=-c_{\nu}c_{\mu}$ for all $\mu\neq\nu\in\{1,\dots,2n\}$. It is straightforward to verify that any such unitaries $U$ generated by quadratics in the Majorana modes form a group, and that this group is exactly that of circuits composed of nearest-neighbor matchgates Jozsa and Miyake (2008). Furthermore, such $U$ preserve the number of Majorana modes, as

[TABLE]

where $\mathbf{u}\in\mathbf{SO}(2n)$ is a $2n\times 2n$ orthogonal matrix. We introduce a Majorana configuration as an ordered tuple of indices $\vec{\alpha}\equiv(\alpha_{1},\alpha_{2},\dots,\alpha_{k})$ of degree $|\vec{\alpha}|\equiv k$ with $\alpha_{j}\in\{1,\dots,2n\}$ and $\alpha_{j}<\alpha_{j+1}$ for all $j\in\{1,\dots,k\}$. The corresponding Majorana configuration operator is the ordered product $C_{\vec{\alpha}}\equiv\prod_{j=1}^{|\vec{\alpha}|}c_{\alpha_{j}}$, with Majorana indices ascending from left to right. Finally, denote by $\mathbf{u}_{\vec{\alpha}\vec{\beta}}$ the submatrix of $\mathbf{u}$ given by taking the rows indexed by $\vec{\alpha}$ and the columns indexed by $\vec{\beta}$, i.e. $(\mathbf{u}_{\vec{\alpha}\vec{\beta}})_{jk}\equiv u_{\alpha_{j}\beta_{k}}$.

Majorana configuration operators transform under matchgate evolution as (see Appendix A)

[TABLE]

That is, the degree of a Majorana configuration operator is preserved, and configuration transition amplitudes are given by determinants of the corresponding single-mode transition submatrices. We also note that

[TABLE]

From Eq. (4), we see that the Heisenberg evolution of $Z_{k}$ will always consist of ${{n}\choose{2}}$ terms, regardless of $k$. However, that of $X_{k}$ will consist of ${{n}\choose{2k-1}}$ terms, which may scale exponentially with $n$ if $k$ also scales with $n$, such as for $X_{\lfloor n/2\rfloor}$ in the center of the chain. This is reflected in the fact that $\langle U^{\dagger}Z_{k}U\rangle$ can always be computed efficiently by a classical computer when the expectation values $\langle C_{\vec{\beta}}\rangle$, for $|\vec{\beta}|=2$, can be, such as for product input Jozsa and Miyake (2008). On the other hand, $U^{\dagger}X_{k}U$ cannot even be stored efficiently on a classical computer in the worst case, so the same strategy will not work \bibnoteThough the distribution of such a measurement can be sampled efficiently, a weaker form of simulation, by the method in Ref. Brod (2016). Nevertheless, localization will provide a means to efficiently approximate this quantity, as we state formally below.

*Analytic results.—*We are able to efficiently calculate the left-hand side of (1) in our setting by our first result

[TABLE]

where $(\mathbf{P}_{\vec{\eta}})_{jk}=0$ if $j\notin\vec{\eta}$ or $k\notin\vec{\eta}$, and $(\mathbf{P}_{\vec{\eta}})_{jk}=\delta_{jk}$ otherwise (i.e. $\mathbf{P}_{\vec{\eta}}$ is the projector onto the modes $\vec{\eta}$). The sign factor is simply $(-1)^{|\vec{\alpha}||\vec{\eta}|+1}$.

The result follows from the Cauchy-Binet formula (see Appendix B). In fact, it is possible to modify the Cauchy-Binet formula to obtain an analyltic calculation of the OTO correlator when $A\equiv\bm{n}_{s}\cdot\bm{\sigma}_{s}$ is any single-site Pauli observable (Appendices C - E). This allows us to regard the OTO correlator as a quadratic form $\mathbf{M}_{s}$, as

[TABLE]

When the bound (1) holds, we can efficiently approximate the evolution in Eq. (4) by truncating the sum to those $\vec{\beta}$ whose support lies strictly within a constant subset of qubits. We model this truncation by the action of a completely depolarizing channel

[TABLE]

which takes any single-qubit operator to its identity component. With our second result, we show that this truncation incurs a bounded error in the average case (see Appendix F):

[TABLE]

where $\overline{(\cdot)}$ denotes an average over a product basis whose Bloch axes are orthogonal to the vectors $\{\mathbf{n}_{s}\}_{s}$, and $\mathbf{M}_{s}$ is as defined in Theorem 1.

*Numerical example.—*Theorem 1 is valid for every unitary time-evolution satisfying (3). However, we will narrow the focus of our numerical analysis to two specific Hamiltonian models, each of the form:

[TABLE]

In Model 1, we allow the local disorder to fluctuate in time about some mean static disorder, and in Model 2, we allow interactions to fluctuate in space and time about mean translationally invariant interactions. The flucutations are chosen as independent, identically distributed random samples taken from the interval $[-\Delta,\Delta]$ every period $\delta t=0.25$. We vary the strength of the mean value and the fluctuation strength $\Delta$ for each model, keeping the remaining parameter fixed. This is intended to resemble a discrete-time control setup, wherein some of the parameters are constrained but others may be varied with some control strength. The static limit, for which $\Delta=0$, is well-understood (see e.g. Ref. Burrell and Osborne (2007)) and will provide a convenient reference point. These models are summarized in Table 1.

As the Hamiltonian (8) is quadratic in the Majorana operators, its time-evolution operator may be expressed as a matchgate circuit Jozsa and Miyake (2008) and so may apply our Theorem 1. In Fig. 1, we plot representative profiles of the OTO correlator (hereafter referred to as “light cones”) for $B=Z_{50}$ (left) and $B=X_{50}$ (right) for $n=100$ in Model 1 in its ballistic (top), diffusive (middle), and localized (bottom) phases. In the localized propagation, for which static disorder $\nu=2$ and $\Delta=0$, the bound (1) is satisfied, and the observable support remains confined. As we increase fluctuations relative to static disorder in the middle plots, for which $\nu=0.75$ and $\Delta=1$, we see that time-dependent fluctuations induce a transition to diffusive propagation.

We identify the propagation phase of each profile by the taking principal singular component of its light cone, treated as a numerical matrix (see Appendix G). We argue that this gives an operationally meaningful, robust, and numerically inexpensive means of extracting the envelope and decay profile.

We characterize the propagation phase by fitting the principal temporal component of each lightcone to a polynomial and extract the exponent of the leading-order term $t^{m}$. In Fig. 2, we plot $m$ for $Z_{50}$ and $X_{50}$ for our two models for $n=100$, as in Fig. 1, as phase diagrams. We identify the ballistic phase with regions where $m$ is very nearly one, the localized phase with regions where $m$ is very nearly zero, and the diffusive phase with regions where $m$ is nearly $0.5$. With this identification, we see that as $\Delta\rightarrow 0$, our results agree with the known limit of static local disorder in Model 1. Similarly, as $\Delta$ becomes large, we see the emergence of a diffusive phase, which is consistent with the results put forth in Burrell et al. (2009). Finally, we see that, for small $\Delta\neq 0$, the localized phase survives. This indicates the existence of new matchgate circuits for which localization may be applied to classically simulate $\langle U^{\dagger}X_{k}U\rangle$ for arbitrary $k$ in a general product state input.

*Discussion.—*We have presented examples where localization may be applied as a tool for classically simulating quantum circuits which were a priori believed to be classically intractable. This is achieved by an analytic calculation of the OTO correlator (presented in Appendix E), followed by truncation to a subset of qubits for which this quantity falls below a certain threshold.

One advantage of our method is that it gives the interior of the light cone in addition to its envelope. In each phase, we see that the light cone interior for $X_{50}$ generally has a higher value than that for $Z_{50}$. This is a consistent difference between the profiles of these operators, which may have important consequences for the complexity required to exactly simulate their expectation values, in a similar fashion with Roberts and Yoshida (2017). We attribute the emergence of a near-ballistic region in the $X$ phase diagram of Model 2, which is absent from the $Z$ diagram, to this observation. Though some amplitude propagates ballistically for both observables in this region, this only manifests as a spreading of the exponential tails for $Z$. For $X$ however, this is exhibited as a ballistic spreading of the high-amplitude region due to interference between its many constituent Majorana operators. This indicates that, at least in the presence of fluctuating interactions, propagation behavior between different local operators can be strikingly different.

We empirically observe, however, that the difference between $X$ and $Z$ observables only emerges at late times, in the saturation value of the OTO correlator for each of these operators. By examining the constant-position slices from the light cones in Fig. 1 as functions of time, we observe a characteristic exponential early-time behavior for these values, which is identical between $X$ and $Z$ propagation. Only as the growth of these operators saturate to roughly their constant values do the differences emerge. This suggests that the evolution of low-degree Majorana configurations may be useful as a good heuristic to observe the lightcone envelope, in a similar spirit to the treatment given in Ref. Xu and Swingle (2018) for finding a low bond-dimension matrix product operator approximation to the evolution of such observables in interacting-fermion systems.

Although we chose here so-called matchgate circuits, related to the time evolution of free fermions, because of their correspondence to Anderson localization, our method is expected to have further applications and extensions. On the former, for example, one may apply it to other random-circuit ensembles, such as those with Haar random matchgates, to study scrambling in general. A preliminary analysis indicates that propagation in this case seems to scale logarithmically, rather than polynomially. On the latter, it may be feasible to extend our method to analyze universal quantum computation by considering matchgate circuits acting on certain entangled input states, such as those of Ref. Bravyi (2006), as our method is independent of the input. In a similar way as Anderson localization has been extended to many-body localization, certain perturbative analysis (in analogy to that performed in Bravyi and Gosset (2016) for Clifford circuits) could probe possible dynamical localization in general quantum circuits as well as simulating interacting fermions. We report some progress on this direction in an upcoming work.

This work was supported in part by National Science Foundation grants PHY-1521016.

Note Added: Since the original submission of this work, the authors became aware of several new related works relating to scrambling in random quantum circuits Nahum et al. (2017), matrix-product appproximations to the OTO correlator in the mixed-field quantum Ising chain Xu and Swingle (2018), and exact calculation of the OTO correlator for the transverse-field quantum Ising chain Lin and Motrunich (2018). We expect that ours will contribute alongside these insightful papers to develop a unified picture of this growing field.

I Outline of Supplementary Material

In the following sections, we prove Theorems 1 and 2 in the main text. In the first section, we show that the transition amplitudes between Majorana configuration operators – ordered products of Majorana operators – under matchgate evolution are given by determinants of submatrices of the single-Majorana transition matrix. We use this to simply prove Theorem 1 in the main text. Next, we show a modification to the Cauchy-Binet formula to compute sums of only those configuration transition amplitudes which involve a fixed background of Majorana operators. In Appendix D, we apply this formula to compute a sum of only those configuration transition amplitudes which involve a fixed parity on a given subset. Finally, we use the results of Appendices A-D to arrive at our analytic calculation of the infinite-temperature OTO correlator with respect to any single-qubit observable, Eq. (6) in the main text, in Appendix E. In Appendix F, we prove a bound on the average-case change in expectation value induced by the depolarizing channel in terms of this quantity (Theorem 2). Finally, we elaborate on some of our numerical techniques in Appendix G.

II Summary of Notation

Additionally, when we say we have a matchgate unitary $U$ on $n$ qubits, we associate it with the $2n\times 2n$ orthogonal matrix $\mathbf{u}\in\mathbf{SO}(2n)$ by the corresponding lowercase symbol in bold. We also denote block matrices in the usual way [e.g. \left(\begin{array}[]{c c}\mathbf{A}&\mathbf{B}\\ \mathbf{C}&\mathbf{D}\end{array}\right)].

III Appendix A: Majorana Configuration Operator Transition Amplitudes

Here we prove the formula

[TABLE]

for some matchgate unitary $U$, by induction on the number of Majorana factors $k\equiv|\vec{\alpha}|$. First consider the case where $k=2$, and let $\alpha_{1}<\alpha_{2}$. We have

[TABLE]

where, from the third to the fourth line, we used the Majorana algebra anticommutation relations

[TABLE]

relabeling dummy indices $\beta_{1}\leftrightarrow\beta_{2}$ on the third sum in the third line. From the fourth to the fifth line, we see that the identity term vanishes as its coefficient is the inner product between two distinct row vectors of an orthogonal matrix. This proves the statement for $k=2$. Next, we assume the statement holds for general $k$ and use this assumption to prove the statement for $k+1$. Without loss of generality, assume $\alpha_{j}<\alpha_{k+1}$ for all $j\leq k$. We now have

[TABLE]

Each of the terms in the sum above falls into one of two categories. Either (i) $\beta_{k+1}\in\vec{\beta}$, and $C_{\vec{\beta}}c_{\beta_{k+1}}=\pm C_{\vec{\beta}\backslash\beta_{k+1}}$, or (ii) $\beta_{k+1}\notin\vec{\beta}$, and $C_{\vec{\beta}}c_{\beta_{k+1}}=\pm C_{\vec{\beta}\cup\beta_{k+1}}$, with sign given in both cases by $(-1)^{|\{j\leq k|\beta_{j}>\beta_{k+1}\}|}$. We first proceed to demonstrate that all of the terms in category (i) vanish. Fix a particular configuration operator $C_{\vec{\beta}\backslash\beta_{k+1}}$. The coefficient on this operator in the r.h.s. of Eq. (S2) is given by a sum over all indices $\gamma$ which could have been removed from $\vec{\beta}$ to yield $\vec{\beta}\backslash\beta_{k+1}$

[TABLE]

where we were able to cancel any sign factors on the terms inside the sum by reordering columns in $\mathbf{u}$ so that the $\gamma$ column appears at the rightmost position, using the alternating sign property of the determinant. The equality is due to the fact that if $\gamma\in\vec{\beta}\backslash\beta_{k+1}$, then the determinant in that term evaluates to zero. Finally, we use the multilinearity property of the determinant to bring the sum on to the last column, as

[TABLE]

i.e. the determinant of a matrix whose last column vector is $\sum_{\gamma}u_{\alpha_{k+1}\gamma}\mathbf{u}_{\vec{\alpha}\gamma}$. The $l$th element of this column is given by

[TABLE]

again following from the fact that this sum is the inner product between two column vectors of an orthogonal matrix. However, $\alpha_{k+1}>\alpha_{l}$ by assumption, so this sum is actually always zero and the determinant in (S4) vanishes. Each of the terms in category (i) therefore vanishes, and the only terms in the r.h.s of Eq. (S2) that survive are in category (ii). We examine these terms by next fixing a particular configuration operator $C_{\vec{\beta}\cup\beta_{k+1}}$. The coefficient on this operator in the r.h.s. of (S2) is given by a sum over all indices $\gamma$ that could have been added to $\vec{\beta}$ to yield $\vec{\beta}\cup\beta_{k+1}$ (we cannot cancel sign factors this time)

[TABLE]

To see that this equality indeed holds, relabel indices in $\vec{\beta}\cup\beta_{k+1}$ such that $\beta_{i}<\beta_{i+1}$ for all $i\leq k$, and suppose $\gamma=\beta_{s}$ in this labeling. Then we have

[TABLE]

As $\alpha_{j}<\alpha_{k+1}$ for all $j\leq k$, this is exactly the sign factor that would appear had we expanded along the last [$(k+1)$st] row of the matrix in the r.h.s. above, since $\beta_{s}$ appears as the $s$th column of this matrix. We therefore have

[TABLE]

which proves the statement for $k+1$, given that it holds for $k$. This completes our inductive proof of Eq. (S1). $\square$

IV Appendix B: Proof of Theorem 1

Let $A$ and $B$ be two Pauli operators (i.e. $A^{\dagger}=A$, $B^{\dagger}=B$, and $A^{2}=B^{2}=I$) and $U$ be a matchgate unitary, then we have

[TABLE]

From the first to the second line, we used Hermiticity of $A$ and $B$, and we expanded the expresion in the third line. From the third to the fourth line, we used $A^{2}=B^{2}=I$, and from the fourth to the fifth line, we used $\operatorname{tr}(I)=2^{n}$ on $n$ qubits. Dividing by a normalization of $2^{n+2}$ on both sides, we have

[TABLE]

Next, we assume $A=i^{a}C_{\vec{\eta}}$ and $B=i^{b}C_{\vec{\alpha}}$ for some integers $a$ and $b$, such that $i^{a}$ and $i^{b}$ ensure that $A$ and $B$ are Hermitian, respectively. E.g., for $A$, this implies

[TABLE]

where $C^{\dagger}_{\vec{\eta}}=(-1)^{\frac{1}{2}|\vec{\eta}|(|\vec{\eta}|-1)}C_{\vec{\eta}}$, since the individual Majorana modes are Hermitian and we need $\frac{1}{2}|\vec{\eta}|(|\vec{\eta}|-1)$ anticommutations between individual distinct modes to reverse the order of $|\vec{\eta}|$ modes. Since $A^{\dagger}=A$, Eq. (S6) implies

[TABLE]

and multiplying both sides of this equation by $i^{a}C^{\dagger}_{\vec{\eta}}$ gives

[TABLE]

Thus, $i^{2a}=(-1)^{\frac{1}{2}|\vec{\eta}|(|\vec{\eta}|-1)}$, and similarly, $i^{2b}=(-1)^{\frac{1}{2}|\vec{\alpha}|(|\vec{\alpha}|-1)}$. Continuing from the dynamical term in Eq. (S5), we have

[TABLE]

From the second to the third line, we used

[TABLE]

and similarly from the third to the fourth line, as well as the fact that $|\vec{\beta}|=|\vec{\alpha}|$. We recognize that

[TABLE]

where $\mathbf{P}_{\vec{\eta}}$ is the projector onto modes $\vec{\eta}$. We therefore have

[TABLE]

which follows from the Cauchy-Binet formula. This therefore proves the theorem

[TABLE]

$\square$

V Appendix C: Modified Cauchy-Binet Formula

Here we prove

[TABLE]

for $\vec{S}$ disjoint from $\vec{B}$. We first rearrange rows and columns inside the matrices $\mathbf{u}$ and $\mathbf{v}$ in the l.h.s. of (S10) to bring each of them to a fiducial form, $\mathbf{u}^{\prime}$ and $\mathbf{v}^{\prime}$, respectively. These are such that $\mathbf{u}^{\prime}_{\vec{\alpha}^{\prime}\vec{S}^{\prime}}=\mathbf{u}_{\vec{\alpha}\vec{S}}$ and $\mathbf{u}^{\prime}_{\vec{\alpha}^{\prime}\vec{B}^{\prime}}=\mathbf{u}_{\vec{\alpha}\vec{B}}$ (and similarly for $\mathbf{v}^{\prime}$), for $\vec{\alpha}^{\prime}\equiv[|\vec{\alpha}|]$ and $\vec{S}^{\prime}\equiv[|\vec{S}|]$. That is, we bring the rows $\vec{\alpha}$ to the top and the columns $\vec{S}$ to the left inside the matrices $\mathbf{u}$ and $\mathbf{v}$ without changing the internal ordering of these tuples, nor the ordering of $\vec{B}$. This is done purely for convenience of presentation and will not affect the argument, as we will undo the rearrangement in the end. We will continue to refer to the numbers of elements in these rearranged tuples by those of their unprimed counterparts (i.e. using $|\vec{S}|$ instead of $|\vec{S}^{\prime}|$), as they are equal. Since this rearrangement is done for both $\mathbf{u}$ and $\mathbf{v}$, any resulting sign factor acquired due to the alternating sign property of the determinant will cancel in the product, and we have

[TABLE]

We will next need the Laplace “expansion by complimentary minors” formula

[TABLE]

where $\varepsilon^{\vec{H},\vec{L}}=(-1)^{\sum_{j=1}^{k}\left(H_{j}+L_{j}\right)}$, $\vec{L}$ is a fixed subset of the columns of $\mathbf{u}$ of size $k$, the sum is over all subsets $\vec{H}$ of rows of $\mathbf{u}$ of size $k$, and $\overline{L}$ and $\overline{H}$ are the set-complements of $\vec{L}$ and $\vec{H}$ in the sets of all columns and rows of $\mathbf{u}$, respectively. This is the analogous formula to expanding the determinant by minors of a fixed column, generalized to a subset of columns $\vec{L}$. Applying Eq. (S12) to the columns $\vec{S}^{\prime}$ of $\mathbf{u}^{\prime}$ and $\mathbf{v}^{\prime}$ on the r.h.s. of (S11) gives

[TABLE]

We next apply the Cauchy-Binet formula to the sum in square brackets, as

[TABLE]

Notice that the matrix in the determinant of the r.h.s. above is simply $\mathbf{u}^{\prime}_{\vec{\alpha}^{\prime}\vec{B}^{\prime}}\mathbf{v}^{\prime\rm{T}}_{\vec{B}^{\prime}\vec{\alpha}^{\prime}}$ with the rows $\vec{H}$ and columns $\vec{L}$ removed (i.e. instead of removing the rows and columns and then multiplying, we can multiply and then remove rows and columns from the product). This gives

[TABLE]

The next step is to “put back in” the rows $\vec{H}$ and columns $\vec{L}$. We do this by treating (S14) as an expansion of the determinant of a larger matrix by complimentary minors of rows $\mathbf{v}^{\prime\rm{T}}_{\vec{S}^{\prime}\vec{\alpha}^{\prime}}$ and then columns $\mathbf{u}^{\prime}_{\vec{\alpha}^{\prime}\vec{S}^{\prime}}$. Working backwards, we see that the sum over $\vec{H}$ in Eq. (S14) evaluates to

[TABLE]

To put the rows back in, we note that

[TABLE]

where $\vec{L}^{\prime}$ is related to $\vec{L}$ by $L^{\prime}_{j}=L_{j}+|\vec{S}|$ for all $j\in\{1,2,\dots,|\vec{L}|=|\vec{S}|\}$. Shifting the columns over by $|\vec{S}|$ inside the determinant gives the overall factor of $(-1)^{|\vec{S}|}$. Thus, the sum over $\vec{L}$ in Eq. (S14) evaluates to

[TABLE]

In the first equality, we apply Eq. (S16) together with the fact that, for each term in the sum for which $\vec{L}$ contains any of the first $|\vec{S}|$ columns, the matrix in the first determinant factor of that term contains at least one column of all zeros, and so the term evaluates to zero. This brings the second line to a form which we recognize to be an expansion by complementary minors of the rows $\vec{S}^{\prime}$ in the larger matrix in the third line. Finally, we use $\mathbf{u}^{\prime}_{\vec{\alpha}^{\prime}\vec{S}^{\prime}}=\mathbf{u}_{\vec{\alpha}\vec{S}}$ (and similarly for the other submatrices) to undo our initial row and column rearrangements and therefore obtain the formula, Eq. (S10). $\square$

VI Appendix D: Fixed-Parity Sum

Here we prove

[TABLE]

where $\vec{S}$ is disjoint from $\vec{B}\equiv(\vec{B}_{l},\vec{B}_{r})$. The statement follows simply from

[TABLE]

In the first equality, we used the fact that $(-1)^{|\vec{\beta}_{r}|}=1$ in the former sum, and $(-1)^{|\vec{\beta}_{r}|}=-1$ in the latter. In the second, we simply combined the sums over all even-sized $\vec{\beta}_{r}$ and all odd-sized $\vec{\beta}_{r}$ with the same summand into the sum over all $(\vec{\beta}_{l},\vec{\beta}_{r})$. We next apply the steps we used to prove the modified Cauchy-Binet formula, except now applying

[TABLE]

to evaluate the sum with the sign factor $(-1)^{|\vec{\beta}_{r}|}$ which appears in the place of Eq. (S13). This gives

[TABLE]

from Eq. (S10). Applying this, together with

[TABLE]

we solve for the sums over even $|\vec{\beta}_{r}|$ and odd $|\vec{\beta}_{r}|$ individually to obtain Eq. (S21) above. $\square$

VII Appendix E: Exact Calculation of the OTO Correlator – Proof of Equation 6

Here we give an explicit calculation of the matrix $\mathbf{M}_{s}$, defined implicitly by

[TABLE]

for some single-spin operator $\mathbf{n}_{s}\cdot\bm{\sigma}_{s}$ acting on qubit $s$, matchgate unitary $U$, and Majorana configuration $C_{\vec{\alpha}}$, where $||A||^{2}\equiv\operatorname{tr}{\left(A^{\dagger}A\right)}$. Since we are only considering a single spin $s$, we drop the spin labels on $\mathbf{n}\equiv\mathbf{n}_{s}$ and $\mathbf{M}\equiv\mathbf{M}_{s}$ in this section for convenience, choosing to label the components of these objects by subscripts instead. We begin by expanding the r.h.s., using this definition and Eq. (S1)

[TABLE]

Thus

[TABLE]

We can obtain the diagonal elements of $\mathbf{M}$, for which $k=k^{\prime}$, from Theorem 1. It remains, then, to calculate the off-diagonal elements of $\mathbf{M}$. Clearly, $\mathbf{M}=\mathbf{M}^{\mathrm{T}}$ by the cyclic property of the trace and the fact that the sum over $\vec{\beta}$ and $\vec{\beta}^{\prime}$ is symmetric in these indices. Furthermore, we must have $M_{xz}=M_{yz}=0$ since there are no $\vec{\beta}$ and $\vec{\beta}^{\prime}$ for which $|\vec{\beta}|=|\vec{\beta}^{\prime}|$ and which describe the same Pauli string on every spin except $s$, with a $X_{s}$ (or $Y_{s}$) present on one and a $Z_{s}$ present on the other. This is because the former are always described by exactly one Majorana operator and the latter by either two or zero Majorana operators. Therefore, we are only left to calculate $M_{xy}.$ In this case, we have

[TABLE]

The only nonvanishing terms will be those for which $\vec{\beta}$ and $\vec{\beta}^{\prime}$ describe the same Pauli string on every spin except for $s$ and for which there is a $Y_{s}$ present for $C_{\vec{\beta}}$ and an $X_{s}$ present for $C_{\vec{\beta}^{\prime}}$. We examine the conditions under which an operator $\sigma^{k^{\prime}}_{s}$ will be present in $C_{\vec{\beta}}$. Let $\vec{\beta}=\left(\vec{\beta}_{l},\vec{\beta}_{s},\vec{\beta}_{r}\right)$, where $\vec{\beta}_{l}\subset[2(s-1)]$ consists of all indices less than $2s-1$, and $\vec{\beta}_{r}\subset\overline{[2s]}$ consists of all indices greater than $2s$ (corresponding to the spins to the left and right of $s$, respectively). We have

[TABLE]

When $|\vec{\beta}_{r}|$ is odd, then $C_{\vec{\beta}}$ and $C_{\vec{\beta}^{\prime}}$ will have a relative phase of $-1$ between them since $X_{s}Z_{s}=-iY_{s}$, but $Y_{s}Z_{s}=iX_{s}$. Additionally taking into account the fact that $[X_{s},Y_{s}]=-[Y_{s},X_{s}]$ gives

[TABLE]

therefore completing the calculation. For completeness, we list the elements of $\mathbf{M}$, once again, below

[TABLE]

VIII Appendix F: Proof of Theorem 2

Here we prove the statement

[TABLE]

where the average $\overline{(\cdot)}$ is taken over a product basis $\{\otimes_{j=1}^{n}|(-1)^{\ell_{j}}\mathbf{n}^{\bot}_{j}\rangle\}_{\bm{\ell}\in\{0,1\}^{\times n}}$ whose Bloch axes are orthogonal to the $\{\mathbf{n}_{s}\}_{s}$, and $\mathcal{E}_{s}$ is the depolarizing channel on qubit $s$, given by

[TABLE]

This channel has the effect of projecting $O$ onto its component which acts only as the identity on site $s$. We further assume that $O$ is a traceless Hermitian operator. We first show that the left-hand side of Eq. (S30) is bounded by the trace norm for any averaging orthonormal basis set $\{|j\rangle\}_{j}$ and Hermitian operator $A$, as

[TABLE]

From the first to the second line, we used the fact that we assume $A$ to be a Hermitian operator, and so can be diagonalized as $A=U^{\dagger}DU$, for $D$ a diagonal matrix. From the second to the third line, we applied the triangle inequality, and from the third to the fourth line, we rearranged sums and used the fact that $\sum_{j=1}^{2^{n}}|U_{kj}|^{2}=1$. The resulting quantity is the trace norm.

Without loss of generality, let $S=\{1,2,\dots,|S|\}$. Then we have, for any input state

[TABLE]

where, from the first to the second line, we expressed the difference as a telescoping sum (let $\otimes_{j=1}^{0}\mathcal{E}_{j}(O)\equiv O$) and applied the triangle inequality in the third line. Next, we use the fact that, for any single-spin input $|\mathbf{n}^{\bot}\rangle$, the depolarized operator expectation value is the same as that of the dephased operator in any basis orthogonal to the input. That is, let $\left(\mathbf{n}\cdot\bm{\sigma}\right)|\mathbf{n}^{\bot}\rangle=e^{i\phi}|-\mathbf{n}^{\bot}\rangle$ for some phase $\phi$, and we have

[TABLE]

on input state $|\mathbf{n}^{\bot}\rangle$ on qubit $s$. Furthermore, taking the same quantity for input state $|-\mathbf{n}^{\bot}\rangle$ has the effect of exchanging $O$ and $\left(\mathbf{n}\cdot\bm{\sigma}_{s}\right)O\left(\mathbf{n}\cdot\bm{\sigma}_{s}\right)$ inside the absolute value, leaving it unchanged. Thus, we average over a product basis of $\{\otimes_{j=1}^{n}|(-1)^{\ell_{j}}\mathbf{n}^{\bot}_{j}\rangle\}_{\bm{\ell}\in\{0,1\}^{\times n}}$ on both sides of Eq. (S33) and apply Eq. (S34) to obtain

[TABLE]

using the bound (S32) in the second line and the fact that the depolarizing channel cannot increase the trace distance between operators in the third. Finally, we apply the operator norm inequality

[TABLE]

where $r\leq 2^{n}$ is the rank of $A$, together with the fact that

[TABLE]

for $||\cdot||$ the Frobenius norm, to arrive at our result

[TABLE]

We remark that the bound is a sum of terms which are always between zero and one, and by Markov’s inequality, we can use this result to bound the fraction of these basis states for which the change in expectation value exceeds some threshold $\delta$

[TABLE]

Since we expect each of these terms to be exponentially decaying outside of the lightcone, depolarizing in this region should only induce a change in expectation value greater than the threshold for some exponentially small set of input states. This affords us exponential precision so long as we are willing to tolerate a fixed fraction of these “pathological” states. Finally, since we had a freedom in choosing which basis over which to depolarize in Eq. (S34), we may optimize each of the $\{\mathbf{n}_{s}\}_{s}$ over the subspace for which they are orthogonal to the Bloch axes of the input to make this bound as restrictive as possible.

IX Appendix G: Numerical Analysis

We characterize the propagation of $B(t)$ by taking the singular value decomposition in discrete time

[TABLE]

where $t_{k}\equiv k\delta t$, and $\lambda_{1}>\lambda_{2}>\dots$ for the singular values $\lambda_{j}$. We take the function $u_{1}(t)$ to represent the light cone envelope and $v_{1}(s)$ the decay profile of $||[\mathbf{n}_{s}\cdot\bm{\sigma}_{s},B(t)]||$ outside of this envelope. This method has several advantages: (i) the principal singular value component is the closest product approximation to the light cone, and so has the form of the right-hand-side of Eq. (1) in the main text; (ii) the singular values themselves give the error incurred in the approximation; and (iii) the principal singular value component is robust to fluctuations from specific disorder realizations, greatly reducing the number of samples needed. It therefore gives an operationally meaningful, robust, and numerically inexpensive means of extracting the envelope and decay profile, which is completely general beyond the setting of matchgate circuits considered here.

In Fig. 3, we plot the results of our analysis for the $X$ light cones of Fig. 1 in the main text. We see that the envelopes (linear scale on top left, and log-log scale on top right) $u_{1}(t)$ propagate as polynomials with different exponents $m$, which we take to be indicative of the dynamical phases of these profiles.