Treelike interactions and fast scrambling with cold atoms

TL;DR

This paper introduces a cold atom quantum spin model with tunable geometry that transitions from linear to treelike structures, enabling the study of fast scrambling and entanglement dynamics in different geometries.

Contribution

It presents a feasible experimental setup for observing emergent geometries and fast information spreading in quantum many-body systems.

Findings

Peak entanglement at the transition point

Exponential quantum information spreading in treelike regime

Observable geometry transition in quench dynamics

Abstract

We propose an experimentally realizable quantum spin model that exhibits fast scrambling, based on non-local interactions which couple sites whose separation is a power of 2. By controlling the relative strengths of deterministic, non-random couplings, we can continuously tune from the linear geometry of a nearest-neighbor spin chain to an ultrametric geometry in which the effective distance between spins is governed by their positions on a tree graph. The transition in geometry can be observed in quench dynamics, and is furthermore manifest in calculations of the entanglement entropy. Between the linear and treelike regimes, we find a peak in entanglement and exponentially fast spreading of quantum information across the system. Our proposed implementation, harnessing photon-mediated interactions among cold atoms in an optical cavity, offers a test case for experimentally observing the emergent geometry of a quantum many-body system.

Equations

$$H=\frac{1}{2S}\sum _{i,j}J(i-j)S_i^{+}{S}_j^{-},$$

$$J(i-j)=\{\begin{array}{ll}{J}_s{2}^{\ell s}& \text{when}|i-j|=2^{\ell },\ell =0,1,2,3\dots \\ 0& \text{otherwise},\end{array}$$

$$E(k)=2J_s\sum _{\ell =0}^{{\log }_2(N/2)}{2}^{\ell s}\cos (2^{\ell }k),$$

$$N\frac{k_{\mathcal{M}}}{2\pi }=\mathcal{M}\left(N\frac{k}{2\pi }\right).$$

$$C(i,j;t)=⟨|[S_i^z(0),S_j^z(t)]{|}^2⟩/S^2,$$

$$C_{\mathrm{cl}}(i,j;t)=\frac{1}{S^2}⟨{\left(\frac{d{S}_j^z(t)}{d{\varphi }_i}\right)}^2⟩,$$

$$J(d=2^{\ell })\propto \frac{n{\Omega }_{\ell }{\Omega }_0{g}^2}{{\Delta }^2\delta },$$

$${\Gamma }_{\mathrm{sc},0}=\frac{{\Omega }_0^2}{4{\Delta }^2}\Gamma ,$$

$$J(d=2^{\ell })\sim {\beta }_{\ell }n\eta \frac{{\Gamma }_{\mathrm{sc},0}\kappa }{\delta },$$

$$L_{k,\pm }=\sqrt{{\gamma }_{k,\pm }}{S}_k^{\pm },$$

$${\mathbf{S}}_k=\frac{1}{\sqrt{N}}\sum _{j=1}^N{e}^{ik{x}_j}{\mathbf{S}}_j.$$

$${\gamma }_{k,\pm }\sim \frac{\kappa }{\delta }{E}_k.$$

$$\gamma \lesssim {\log }_2(N/2)\frac{J\kappa }{\delta }.$$

$$\rho \equiv \frac{J}{\gamma +B{\Gamma }_{\mathrm{sc},0}}\sim \frac{\beta }{\frac{M\beta \kappa }{\delta }+\frac{B\delta }{n\eta \kappa }},$$

$$\rho =\frac{1}{2}\sqrt{\frac{n\eta \beta }{MB}}=\frac{1}{2}\sqrt{\frac{n\eta \beta /M}{1+2M{\beta }^2}}.$$

$$\rho =\frac{\sqrt{n\eta }/2}{(2M{)}^{3/4}}.$$

$$\frac{n\eta }{4}\gtrsim [2{\log }_2(N/2){]}^{3/2},$$

$$⟨n_j(t)⟩=⟨0|c_i{c}_j^{†}(t)c_j(t)c_i^{†}|0⟩=⟨0|[c_i,c_j^{†}(t)][c_j(t),c_i^{†}]|0⟩$$

$$a{\left[d_{(\mathcal{M})}\right]}^b\le t_{ϵ}$$

$$t_{ϵ}\le a^{\prime }{\left[d_{(\mathcal{M})}\right]}^{b^{\prime }}{\left[\log d_{(\mathcal{M})}\right]}^{c^{\prime }}$$

$${\mathcal{O}}_i(t)={\mathcal{O}}_i(0)-it[{\mathcal{O}}_i,H]-\frac{1}{2}{t}^2[[{\mathcal{O}}_i,H],H]+\cdots$$

$$⟨[{\mathcal{O}}_i(t),{\mathcal{O}}_j(0){]}^2⟩\approx t^2⟨[[{\mathcal{O}}_i,H_{ij}],{\mathcal{O}}_j{]}^2⟩$$

$$C(i,j;t)\propto (Jt{)}^{2r_{ij}}$$

$$\left\{S_i^{\alpha },S_j^{\beta }\right\}=\sum _{\gamma }{S}_i^{\gamma }{ϵ}^{\alpha \beta \gamma }{\delta }_{ij}$$

$${\dot{S}}_i^{\alpha }=\left\{S_i^{\alpha },H\right\}.$$

$${\dot{\mathbf{x}}}_i=f({\mathbf{x}}_i)$$

$$P_n={\mathcal{B}}_n^{[1]}{\mathcal{B}}_n^{[2]}...{\mathcal{B}}_n^{[N]},$$

$$\begin{array}[]{l}{\mathcal{B}}_{n}^{[1]}=\left(\begin{array}[]{ccccc}p_{\downarrow}&p_{\uparrow}&0&0&\ldots\end{array}\right),\\ {\mathcal{B}}_{n}^{[1<i<N]}=\left(\begin{array}[]{ccccc}p_{\downarrow}&p_{\uparrow}&0&0&\ldots\\ 0&p_{\downarrow}&p_{\uparrow}&0&\\ 0&0&p_{\downarrow}&p_{\uparrow}&\\ 0&0&0&p_{\downarrow}&\\ \vdots&&&&\ddots\end{array}\right),\\ {\mathcal{B}}_{n}^{[N]}=\left(\begin{array}[]{c}\vdots\\ 0\\ 0\\ p_{\uparrow}\\ p_{\downarrow}\end{array}\right),\end{array}$$

Full text

The fast scrambling conjecture—inspired by studies of the black-hole information problem—predicts a lower bound on the time for information to spread from one to all degrees of freedom of an $N$-body quantum system, scaling as $t_{*}\propto\log(N)$ Sekino and Susskind (2008); Hayden and Preskill (2007). Fast scrambling is conceptually important as a putative signature of the quantum physics of black holes Maldacena et al. (2016). It is also of practical importance because a fast scrambler is an efficient quantum encoder, capable of quickly entangling quantum information across many physical qubits Preskill (2016). While information scrambling has been probed in several pioneering experiments Gärttner et al. (2017); Li et al. (2017); Wei et al. (2018); Landsman et al. (2019), observing fast scrambling remains an outstanding challenge.

In typical quantum systems found in nature, fast scrambling is precluded by the locality of interactions. In the presence of only short-range interactions, information propagation is bounded by a linear “light cone” Lieb and Robinson (1972); Cheneau et al. (2012); Richerme et al. (2014); Bravyi et al. (2006) (Fig. 1a), such that the scrambling time is at least $t_{*}\propto N$. Even with long-range power-law interactions Hastings and Koma (2006); Foss-Feig et al. (2015); Else et al. (2018), the time for quantum correlations to propagate scales polynomially with distance Foss-Feig et al. (2015); Else et al. (2018) and not logarithmically. Potential work-arounds include engineering interactions among spatially overlapped modes Danshita et al. (2017); Chen et al. (2018); Chew et al. (2017) or using light to mediate effectively non-local interactions Swingle et al. (2016); Marino and Rey (2019). Proposals to date have focused on emulating solvable models for fast scrambling featuring random all-to-all couplings Danshita et al. (2017); Chen et al. (2018); Chew et al. (2017); Marino and Rey (2019).

In this Letter, we propose an experimentally realizable quantum spin model featuring a sparse, deterministic graph of non-local interactions that enables fast scrambling. A modest generalization of the simplest model provides a smooth interpolation between two radically different notions of geometry: linear or treelike, as shown in Fig. 1b. Our simplest model is, in a certain sense, exactly between these two limits. Instead of exhibiting extended spatial geometry, it is fully characterized by a coupling graph whose diameter grows logarithmically with the number of sites, resulting in exponentially fast spreading of perturbations that decreases to power-law spreading as we tune away from the original model. Thus, loosely speaking, we see the exponential spread of perturbations as a form of criticality encountered between two incompatible geometric notions of locality (Fig. 1b).

In studying this transition in geometry, we are also motivated by holography, in particular the recently proposed $p$-adic version Gubser et al. (2017); Heydeman et al. (2016) of the anti-de Sitter / Conformal Field Theory correspondence (AdS/CFT). Here, the role of the gravitational bulk is played by the Bruhat-Tits tree, which is an infinite regular tree with $p+1$ edges leading into each vertex, similar to Fig. 1b for $p=2$. Past works have suggested that diffusion on the Bruhat-Tits tree can be used to understand aspects of fast scrambling near black hole horizons Barbón and Magán (2012, 2013).

The models we consider here are inspired by an elegant proposal Hung et al. (2016) for engineering translation-invariant long-range interactions by encoding spins in cold atoms coupled to light in a waveguide or cavity Gopalakrishnan et al. (2009); Leroux et al. (2010); Strack and Sachdev (2011); Hosten et al. (2016); Welte et al. (2018); Davis et al. (2019); Kroeze et al. (2018); Norcia et al. (2018); Landini et al. (2018); Braverman et al. (2019); Borjans et al. (2019); Mivehvar et al. (2019). This scheme can generate Hamiltonians of the form

[TABLE]

where $i,j=1,\ldots,N$, and $\mathbf{S}_{i}$ is a spin-$S$ operator representing either an individual atom or the collective spin of a small atomic ensemble at site $i$. We consider a sparse graph of couplings between sites separated by powers of two,

[TABLE]

where by tuning the exponent $s$ from $-\infty$ to $+\infty$ we interpolate between the linear and treelike limits. The simplest choice $s=0$, where all non-zero couplings are equal, permits a perturbation to spread exponentially because the number of pairwise interactions required to get from site $i$ to site $j$ is bounded above by the number of binary digits in $i$ and $j$ that differ.

The couplings in Eq. (2) can be generated by a pair of Raman processes wherein one atom flips its spin by virtually scattering a photon from a control field into a cavity, and a second atom rescatters this photon such that its spin flops Davis et al. (2019). In a magnetic field gradient Hung et al. (2016), the energy cost of a flip-flop is proportional to the distance $d=i-j$ between spins. Thus, inducing resonant flip-flops at distance $|d|=2^{\ell}$ requires modulating the control field at a frequency $\omega_{\ell}\propto d$. More generally, weakly modulating the control field at multiple frequencies $\omega_{\ell}$ for $\ell=0,1,2,3,\dots$ produces a set of sidebands whose amplitudes dictate the hopping amplitudes $J(d)$ (Fig. 1).

Below, we take $N$ a power of 2 and assume periodic boundary conditions unless otherwise stated, letting $|x|$ denote the minimum value of $\sqrt{(x+qN)^{2}}$ over all integers $q$. We normalize the couplings in Eq. (2) such that the largest is always a constant $J_{0}$, letting $J_{s}=J_{0}$ for $s\leq 0$ and $J_{s}=J_{0}(N/2)^{-s}$ for $s>0$.

The key features of our spin model are evident already in the spreading of a single initially localized spin excitation (Fig. 2a), which is governed by the magnon dispersion relation (Fig. 2b)

[TABLE]

where $k\in[0,2\pi)$ is the wavenumber. Intriguingly, as $N$ increases, the dispersion relation shows fractal behavior converging to a Weierstrass function for $0<s<1$ 111Related fractal behavior was recently found in a simplified statistical mechanical model with sparse couplings Gubser et al. (2018) of essentially the same form as in Eq. (2).. For arbitrary $s$, Eq. (3) allows us to analytically compute the single-magnon dynamics as shown in Fig. 2a, where we introduce a single excitation on site $i=N/2$ and plot the mean occupation $\langle n_{j}(t)\rangle$ as a function of site number $j$ and time $t$.

The single-particle dynamics reveals the geometry of the interaction graph and its dependence on the exponent $s$. When $s$ is large and negative, the excitation spreads ballistically, as expected for a nearest-neighbor spin chain, producing the linear light cone in Fig. 2a(i). By contrast, for $s>0$, where interactions grow with physical distance, the excitation jumps discontinuously between distant sites (Fig. 2a(iii)). Rather than interpreting the apparent absence of a light cone for $s>0$ as the absence of locality, we argue that a new version of locality emerges based on the $2$-adic norm $|x|_{2}=2^{-v(x)}$, where $2^{v(x)}$ is the largest power of 2 that divides $x$. The distance $|i-j|_{2}$ between sites $i$ and $j$ is called ultrametric because the distance of the sum of two steps is never greater than the larger of the two steps’ distance; by contrast the usual distance $|i-j|$ is called Archimedean because many small steps can be combined into a large jump.

We can understand the $2$-adic norm as a tree-like measure of distance because $|i-j|_{2}=2^{d_{\text{tree}}(i,j)/2}/N$, where $d_{\text{tree}}(i,j)$ is the number of edges between sites $i$ and $j$ along the regular tree in Fig. 1b 222The tree appears to break translation invariance, but in fact it does not: $d_{\text{tree}}(i+1,j+1)=d_{\text{tree}}(i,j)$.. The leaves are numbered in order of increasing ${\cal M}(i)$, where the discrete Monna map ${\cal M}$ reverses the bit order in the site number. For example, for $N=8$ sites, ${\cal M}(1)=4$ because in binary, ${\cal M}(001_{2})=100_{2}$. Noting that $Nk/2\pi$ is an integer, we may likewise define a Monna-mapped wavenumber $k_{\mathcal{M}}$ by

[TABLE]

For large positive $s$, we rearrange the spins according to the Monna map and find that a light cone reappears (Fig. 2a(v)) and the dispersion relation is smoothed out (Fig. 2b), corroborating the transformation to the treelike geometry defined by the $2$-adic norm.

The radical difference between Archimedean and ultrametric geometry raises the question of what happens near $s=0$, where short- and long-range couplings are equally strong. Here, it is useful to think of the sites as arranged on a hypercube 333We thank P. Hayden for discussions on this point., where each site number in binary specifies a corner. The coupling graph consists of the edges of the hypercube plus some diagonals, as shown in Fig. 1b for $N=8$ sites. For any $N$, the graph distance $r_{ij}$ counts the minimum number of edges required to connect sites $i$ and $j$, and is at most $\lceil\frac{1}{2}\log_{2}(N)\rceil$, where $\lceil x\rceil$ denotes the smallest integer greater than or equal to $x$.

We expect the logarithmic diameter of the interaction graph at $s=0$ to enable a localized perturbation to spread over the entire system in a logarithmic time $t\propto\log(N)$. As a first test, in Fig. 2a we examine single-magnon transport as a function of both physical distance $d=|i-j|$ and graph distance $r_{ij}$. We observe spin transport on a timescale that is roughly linear in graph distance (Fig. 2a(iv)), and thus at most logarithmic in physical distance $d$ and system size $N$.

To permit this logarithmic timescale, the polynomial light cone in Fig. 2a(i) must break down as $s$ approaches 0. To examine this breakdown, we evaluate the time $t_{\epsilon}$ for the magnon occupation at a distance $d$ from the initial site to reach a threshold value $\langle n_{i+d}\rangle=\epsilon=1/N^{2}$. For arbitrary $s$, we may always bound $t_{\epsilon}$ from below by a polynomial of the form $a\left[d_{(\mathcal{M})}\right]^{b}\leq t_{\epsilon}$ that depends on physical distance $d$ for $s\leq 0$ or Monna distance $d_{\mathcal{M}}\equiv\mathcal{M}(d)$ for $s>0$, with non-negative constants $a,b$. We determine the dependence of the exponent $b$ on $s$ from a fit to the points of fastest growth SM in Fig. 2a, finding a direct proportionality $b\propto|s|$ (red points in Fig. 2c). The vanishing exponent $b$ at $s=0$ signifies the breakdown of the polynomial light cone.

Furthermore, we verify that even the slowest-growing occupations take only a logarithmic time to reach the threshold $\epsilon$. For arbitrary $s$, the time $t_{\epsilon}$ is bounded above by $t_{\epsilon}\leq a^{\prime}\left[d_{(\mathcal{M})}\right]^{b^{\prime}}\left[\log d_{(\mathcal{M})}\right]^{c^{\prime}}$, with non-negative constants $a^{\prime},b^{\prime},c^{\prime}$ (dashed red lines in Fig. 2a) SM . Here, as $s\rightarrow 0$ the polynomial exponent $b^{\prime}$ again disappears, but the logarithmic exponent $c^{\prime}$ remains (blue points in Fig. 2c). For both the fastest- and slowest-growing occupations we observe an approximate symmetry about $s=0$ that suggests a duality between the Archimedean and non-Archimedean regimes Gubser et al. (2018).

While the spin dynamics provides strong evidence for the underlying geometry of the model as a function of $s$, the structure of entanglement provides an even sharper proxy for locality. We analyze the growth of entanglement for a system initialized in a product state of spins polarized along $\mathbf{\hat{x}}$. Figure 3a shows the entanglement entropy $S_{A}$ of a subsystem $A$ as a function of its size and time, with the spins partitioned either according to their physical position (top) or to their Monna-mapped ordering (bottom). In either of the two limits where $|s|$ is large, there is a natural way of partitioning the system such that entanglement remains low, as summarized in Fig. 3b: for $s<0$, the cut must be between nearest neighbors in the linear chain (blue curve), whereas for $s>0$ the cut must be between branches of the tree (red curve).

Near $s=0$, however, entanglement is high no matter how we cut the system. We verify this by plotting the minimum entanglement entropy over all possible bipartitions, without regard to locality. The resulting gray curves in Fig. 3b, for different partition sizes $L$, show a sharp peak in entanglement at $s=0$. Thus, at the crossover between the Archimedean and $2$-adic geometries, there is no good notion of locality, and all spins are strongly coupled to one another.

Our non-local spin models generically exhibit quantum chaos. One indicator is that energy level spacings in the $s=0$ model at half filling ($N/2$ magnons for $S=1/2$) exhibit random-matrix statistics SM ; Weinberg and Bukov (2017); Wimberger (2014). But is the highly connected model at $s=0$ a fast scrambler? To probe this question, we consider the out-of-time-order correlation function (OTOC) Larkin and Ovchinnikov (1969); Shenker and Stanford (2014); Swingle et al. (2016); Swingle and Yunger Halpern (2018); Zhang et al. (2019)

[TABLE]

where $S^{z}_{i,j}$ are local spin operators at sites $i,j$ and $S^{z}_{j}(t)=e^{iHt}S^{z}_{j}e^{-iHt}$. In a typical fast scrambler, $C$ approaches its saturation value on a timescale $t_{*}\propto\log(N)/\lambda$, where $\lambda>0$ is a Lyapunov exponent quantifying the system’s exponential sensitivity to perturbations.

In systems with local interactions, fast scrambling is precluded by the Lieb-Robinson bound Lieb and Robinson (1972); Bentsen et al. (2019), which restricts $C\lesssim e^{-(d-vt)}$ to exponentially small values at distances $d=|i-j|$ outside a light cone with Lieb-Robinson velocity $v$, thereby preventing saturation of OTOCs until a time $t_{*}\propto N$. By contrast, known models for fast scrambling feature random all-to-all couplings Lashkari et al. (2013); Maldacena and Stanford (2016) and have no sense of spatial locality. Our model at $s=0$ offers an alternative route to fast scrambling: despite its effective light cone vs. graph distance $r_{ij}$ (Fig. 2a(iv)), the early-time growth of OTOCs is permitted to reach values $C\sim e^{-(r_{ij}-vt)}\sim 1/N^{\alpha}$ Bentsen et al. (2019) due to the logarithmic graph diameter $r_{\mathrm{max}}\approx\frac{1}{2}\log_{2}(N)$, where $\alpha$ is a constant of order unity and $v\propto J_{0}\log(N)$ because each spin has $\log(N)$ couplings. Subsequent Lyapunov growth $C\sim e^{\lambda t}/N^{\alpha}$ therefore allows OTOCs to reach saturation in a time $t_{*}\sim\alpha\log(N)/\lambda$.

A Lyapunov regime, however, is not guaranteed or even expected in a finite-size system with small local Hilbert-space dimension, e.g., at spin $S=1/2$ Khemani et al. (2018). We thus first analyze the $s=0$ model semiclassically in the limit where each site contains a spin $S\gg 1$, which is natural to implement experimentally by letting each $\mathbf{S}_{i}$ represent the collective spin of an ensemble. We consider the averaged sensitivity

[TABLE]

for a small initial rotation $\phi_{i}$ of spin $i$ about the $z$-axis, whose correspondence to the OTOC can be observed by replacing the commutator in Eq. (5) with a Poisson bracket Polkovnikov (2010); Cotler et al. (2018). We calculate the average in Eq. (6) for an ensemble in which each spin has a random initial orientation in the $xy$-plane, the classical limit of an infinite-temperature state at half filling.

The growth in sensitivity $C_{\mathrm{cl}}(t)$ generically exhibits two distinct regimes Marino and Rey (2019), visible in Fig. 4. The first is a rapid power-law growth $C_{\mathrm{cl}}(t)\propto(J_{0}t)^{2r_{ij}}$ for $J_{0}t\lesssim 1$. A transition to exponential growth occurs for $J_{0}t\gtrsim 1$. Crucial to fast scrambling is that, by the time the exponential growth begins, the OTOCs have already reached values $C_{\mathrm{cl}}\sim 1/N^{\alpha}$. We expect the subsequent exponential growth $C_{\mathrm{cl}}(t)\sim e^{\lambda t}/N^{\alpha}$ to yield a value $C_{\mathrm{cl}}(t_{*})\sim 1$ at time $t_{*}\sim\alpha\log(N)/\lambda$. To verify this behavior, we first fit the exponential growth for a range of system sizes to obtain the dependence of the Lyapunov exponent $\lambda$ on $N$, and then evaluate the time $t_{*}$ to reach $C_{\mathrm{cl}}(t_{*})=1$ in terms of $\lambda$ (Fig. 4b). Fitting the dependence of this semiclassical scrambling time on $N$ yields $\lambda t_{*}=\alpha\log(N)+\beta$, with $\alpha=1.1(1)$.

Our semiclassical analysis demonstrates that the correlations developing at early times $J_{0}t\lesssim 1$, before the onset of exponential growth, are crucial for enabling fast scrambling: the weakest correlations should be only algebraically small ($C\propto N^{-\alpha}$) and not exponentially small in $N$. We now apply this insight to investigate whether the quantum model at $S=1/2$ can be a fast scrambler.

As the essence of fast scrambling is efficient spreading of information across an exponentially large Hilbert space, the process is intrinsically difficult to study numerically. Nevertheless, we calculate the early-time dynamics for $N=64$ sites at half filling using Matrix Product State (MPS) techniques Schollwoeck (2011); Haegeman et al. (2011); Koffel et al. (2012); Haegeman et al. (2013, 2016); SM , for a system at infinite temperature with open boundary conditions. We see no Lyapunov regime, which is not surprising since there is no small parameter for $S=1/2$ at finite $N$ to prevent rapid saturation of the OTOC (Fig. 4c). However, Fig. 4c(ii) indicates early-time correlations that fall off exponentially with graph distance $r_{ij}$, consistent with correlations across all sites that are algebraically rather than exponentially small in $N$—the aforementioned necessary condition for fast scrambling.

Our results indicate that fast scrambling is accessible in sparsely coupled models without disorder, and might generalize to a wider range of coupling patterns in which the graph diameter grows logarithmically with system size. Near-term cavity-QED experiments offer promise for observing the logarithmic timescale for spin transport at $s=0$, as well as the linear-to-treelike transition. In the large-$S$ regime, where each site contains an atomic ensemble, experiments will benefit from a collective enhancement in the coherence of interactions. Ultimately, implementations of the spin-1/2 model—for sufficiently strong atom-light coupling Colombe et al. (2007); SM —could test for fast scrambling SM ; Swingle et al. (2016); Swingle and Yunger Halpern (2018); Zhang et al. (2019) at large $N\sim 10^{3}$ in the quantum regime.

Complementarily, future theoretical work may investigate whether a Lyapunov regime can emerge in the spin-1/2 quantum model through a coarse-graining procedure that exploits the self-similarity of the coupling pattern. Measures of entanglement throughout the linear-to-treelike transition also merit further study, and may enable a more explicit connection to holography via entanglement wedges or tensor networks Swingle (2012); Pastawski et al. (2015); Bhattacharyya et al. (2018); Heydeman et al. (2016, 2018); Vidal (2008); Shi et al. (2006). Future work may also explore prospects for harnessing the rapid and deterministic generation of entanglement for quantum information processing.

I Treelike Interactions and Fast Scrambling with Cold Atoms:

Supplemental Material

In this supplement, we elaborate on the proposed experimental implementation and on numerical and analytical techniques used to obtain the results in the main text. In Section I we discuss details of the experimental implementation including the effects of dissipation. In Section II we explain the methods used to find upper and lower bounds on magnon occupation growth times $t_{\epsilon}$, present evidence for chaotic level statistics, discuss the short-time expansion of OTOCs, and elaborate on the semiclassical and Matrix-Product State numerical methods used in the main text.

II I. Experimental Implementation

II.1 A. Engineering the Couplings

We propose an experimental implementation in a one-dimensional array of atoms uniformly coupled to the mode of an optical cavity, with vacuum Rabi frequency $2g$. At each of $N$ sites are $n$ atoms that encode a spin of length $S=nF$, either in the Zeeman states of an individual spin-$F$ atom ($n=1$) or in the collective magnetization of a localized ensemble. A magnetic field gradient produces a Zeeman splitting $\omega_{Z}+j\omega$ that depends linearly on the site number $j\in\{1,\dots,N\}$. Spin-exchange interactions are introduced by an optical control field at large detuning $\Delta$ from atomic excited states, as in Ref. Davis et al. (2019). This control field, with center frequency $\omega_{0}$, drives virtual Raman processes wherein an atom absorbs a photon and emits it into the cavity mode at frequency $\omega_{c}$, whence it is rescattered by another atom.

To turn on interactions at each distance $d=2^{\ell}$, the control field is amplitude-modulated at frequency $\omega_{\ell}$ and modulation index $\beta_{\ell}\ll 1$ to generate sidebands of Rabi frequencies $\Omega_{\ell}=\beta_{\ell}\Omega_{0}$, where $\Omega_{0}$ is the Rabi frequency of the carrier. We assume that the detuning $\delta=\omega_{0}-\omega_{c}-\omega_{Z}$ of the carrier from Raman resonance is large compared to all modulation frequencies and to the cavity linewidth $\kappa$, i.e., we assume $\delta>(\log_{2}N)\omega$ and $\delta>\kappa$. The spin-exchange couplings are then given by

[TABLE]

where $n$ is the number of atoms per site.

II.2 B. Effects of Dissipation

To calculate how coherent the interactions can be, we would like to compare the couplings $J$ with the rates of two decay processes: collective decay via photon loss from the cavity mode and single-atom decay due to spontaneous emission. To this end, we start by re-expressing the interaction strengths in terms of the free-space scattering rate due to the control field carrier,

[TABLE]

where $\Gamma$ is the atomic excited-state linewidth. In terms of $\Gamma_{\mathrm{sc,0}}$ and the single-atom cooperativity $\eta=4g^{2}/(\kappa\Gamma)$, we have

[TABLE]

which shows that interactions can dominate over spontaneous emission for sufficiently large cooperativity $n\eta$. Below, in identifying requirements on the cooperativity, we will keep in mind that the total spontaneous emission rate is augmented by the power in the sidebands: $\Gamma_{\mathrm{sc}}=B\Gamma_{\mathrm{sc,0}}$, where $B=1+2\sum_{\ell}|\beta_{\ell}|^{2}$.

A second dissipation mechanism is photon loss from the cavity. Cavity decay leads to collective dissipation described by a set of Lindblad operators

[TABLE]

where we have defined the spin-wave operators

[TABLE]

The decay rates $\gamma_{k,\pm}$ are related to the dispersion relation $E_{k}$ in Eq. (3) of the main text by

[TABLE]

For example, for our $s=0$ model, which has $\log_{2}(N/2)$ couplings of equal strength $J_{\ell}\equiv J$, the decay rates $\gamma_{k,\pm}$ are bounded above by

[TABLE]

II.3 C. Interaction-to-Decay Ratio

We calculate an overall interaction-to-decay ratio at $s=0$ by comparing the couplings $J$ with the rates of collective decay $\gamma$ via the cavity and spontaneous emission $\Gamma_{\mathrm{sc}}$. We thus obtain a ratio

[TABLE]

where $\beta$ is the modulation index for each sideband, $M\equiv\log_{2}(N/2)$ is the number of modulation frequencies, and $B=1+2M\beta^{2}$. At an optimal detuning $\delta/\kappa\sim\sqrt{n\eta M\beta/B}$, we obtain

[TABLE]

Eq. (S9) is maximized by choosing the modulation index such that total power in the sidebands equals the power in the carrier, $2M\beta^{2}=1$, which is consistent with the weak-modulation requirement $\beta\ll 1$ provided that $N$ is large. The interaction-to-decay ratio is then

[TABLE]

Reaching the scrambling time $Jt\sim 1$ requires an interaction-to-decay ratio $\rho\gtrsim 1$. By Eq. (S10), we thus require a large collective cooperativity of the subensemble at each site:

[TABLE]

where the weak scaling of the right-hand side with $N$ is convenient for scaling to large system sizes. For example, Eq. (S10) can be satisfied for a system of $N=2^{10}$ sites with $n=300$ atoms per site and single-atom cooperativity $\eta\sim 1$ (Fig. S2). A higher but demonstrated single-atom cooperativity $\eta\sim 150$ Colombe et al. (2007) suffices to begin exploring fully quantum spin-1/2 models with a single atom per site.

Even when these conditions on the cooperativity are satisfied, dissipation may modify shape of the light cone, as pointed out in Ref. Zhang et al. (2019). This effect can be mitigated by measuring ‘corrected’ OTOCs Zhang et al. (2019); Swingle and Yunger Halpern (2018), which divide out the effects of dissipation and recover the shape of the true light cone.

II.4 D. Effects of Finite Modulation Depth

Our derivation of the couplings $J_{\ell}$ in Sec. I.A. assumes the conceptually simple limit of small modulation depth $\beta_{\ell}\ll 1$. In practice, for the system sizes most readily accessible in experiments, a finite modulation depth $\beta\sim 0.2$ is advantageous for maximizing the interaction-to-decay ratio (Fig. S2). Operating at finite modulation depth adds weak additional couplings (a factor $\beta$ smaller than the intended couplings) at distances that are not powers of two. While these couplings are unlikely to appreciably alter the fast scrambling dynamics or the geometry of the coupling graph, it is worth noting that the modulation waveform can be designed to produce only the desired couplings even at finite modulation depth. In the most general case, for a control field $\Omega(t)$, the amplitude modulation waveform $|\Omega(t)|^{2}$ dictates the magnon dispersion relation $E(k)\propto|\Omega(k/\omega)|^{2}$. Thus, while the field $\Omega(t)\propto 1+2\sum_{\ell}\beta_{\ell}\cos(2^{\ell}\omega t)$ produces precisely the desired coupling pattern in the limit $\beta_{\ell}\ll 1$, a more general class of waveforms $\Omega(t)\propto\sqrt{E(\omega t)+\mathrm{const.}}$ produce the same coupling pattern at finite modulation depth.

III II. Numerical and Analytical Techniques

III.1 A. Magnon Occupation and Information Spreading

To characterize the spreading of local perturbations in the single-magnon sector, we introduce a single excitation at site $i$ and evaluate the time $t_{\epsilon}$ at which the magnon occupation at another site $j$ reaches a threshold value $\langle n_{j}\rangle=\epsilon=1/N^{2}$. We note that the magnon occupation $\langle n_{j}(t)\rangle$, when the system is initialized with one excitation at site $i$, is equivalent to the OTOC

[TABLE]

since $c_{j}(t){\left|0\right\rangle}=0$. (Analogous relations hold for more general lowering operators $A^{-},B^{-}$, so long as $B^{-}{\left|0\right\rangle}=0$ and ${\left|0\right\rangle}$ is an eigenstate of $H$.) The occupation $\langle n_{j}(t)\rangle$ therefore encodes similar information about the spread of information in the system and is bounded by the Lieb-Robinson bound that controls the OTOC.

III.2 B. Characterizing Fastest- and Slowest-Growing Magnon Occupations

In Fig. S3, we plot the threshold time $t_{\epsilon}$ versus physical distance $d=|i-j|$ or Monna-mapped distance $d_{\mathcal{M}}=\mathcal{M}(d)$ from the original site for different values of $s$. While these plots show significant scatter, we observe that one can always bound the values $t_{\epsilon}$ from below by a polynomial light cone of the form

[TABLE]

and from above by a function of the form

[TABLE]

for some non-negative constants $a,b,a^{\prime},b^{\prime},c^{\prime}$. To choose the lower-bound coefficients $a,b$ optimally, we first fit Eq. (S13) to the fastest-growing occupations (black circles in Fig. S3), which always occur at distances $d_{(\mathcal{M})}$ that are powers of 2. The resulting lower bounds are shown in black in Fig. S3.

To find optimal upper-bound coefficients $a^{\prime},b^{\prime},c^{\prime}$, we first isolate the slowest-growing occupations by splitting the scatter plot into bins $d\in[2^{n},2^{n+1})$ for integers $0\leq n\leq\log_{2}N-2$ and picking the largest value of $t_{\epsilon}$ within each bin (red diamonds in Fig. S3). We then use a constrained minimization search to find coefficients $a^{\prime},b^{\prime},c^{\prime}$ that minimize the sum of the squares of the differences between the slowest data points in each bin and the upper-bound function of Eq. (S14). In addition, to ensure that the resulting function bounds the slowest points from above, we place a constraint on the minimization search such that the optimal function must always lie above the slowest datapoints. Finally, to ensure that the upper and lower bounds Eqs. (S13) and (S14) hold simultaneously for large $d_{(\mathcal{M})}$, we demand that $a^{\prime}\geq a$ and $b^{\prime}\geq b$. The resulting optimal upper-bound curves are plotted in red in Fig. S3.

III.3 C. Energy Level Statistics

To verify that the $s=0$ model is chaotic in the interacting quantum regime, we calculate the energy level statistics for a system of $N=16$ sites at half filling by exact diagonalization Weinberg and Bukov (2017). Figure S4 shows the distribution of spacings between neighboring energy levels, after unfolding the spectrum in each momentum sector according to the procedure in Ref. Wimberger (2014). The distribution is consistent with the random-matrix statistics of the Gaussian Orthogonal Ensemble.

III.4 D. Short-Time Expansion of OTOCs

Here we show that OTOCs display power-law growth $C(t)\propto(Jt)^{2r_{ij}}$ at early times $Jt\ll 1$ Marino and Rey (2019), where $J$ represents the typical energy scale of the Hamiltonian $H$. Consider the OTOC $C(i,j;t)=\langle[\mathcal{O}_{i}(t),\mathcal{O}_{j}(0)]^{2}\rangle$, where $\mathcal{O}_{i}$ is an arbitrary operator at site $i$. At early times, the growth of OTOCs can be understood via a short-time expansion:

[TABLE]

If there is a direct coupling between sites $i$ and $j$, then at short times $t\ll J^{-1}$ the above expansion is simply dominated by the first two terms, leading to:

[TABLE]

where $H_{ij}$ is the term in the Hamiltonian coupling sites $i\neq j$. If there is no direct coupling between $i,j$, then this lowest-order quadratic term vanishes and we must consider higher-order terms in the expansion. More generally, if $r_{ij}$ is the minimum number of couplings $H_{nm}$ required to hop from site $i$ to site $j$, then all lowest-order terms $\sim t^{2n}$ for $n<r_{ij}$ vanish, and the early-time growth is dominated by

[TABLE]

where the constant of proportionality includes contributions from nested commutators $\langle[[[[\mathcal{O}_{i},H],H],\ldots,H],\mathcal{O}_{j}]^{2}\rangle$. For example, for the interaction graph at $s=0$, we find that $C(0,1;t)$ and $C(0,2;t)$ grow like $(J_{0}t)^{2}$ at early times since there is a direct coupling between sites $0,1$ and sites $0,2$. Similarly, $C(0,3;t)$ grows like $(J_{0}t)^{4}$ because it requires two couplings to hop between sites 0 and 3, and $C(0,11;t)$ grows like $(J_{0}t)^{6}$ because it requires three couplings to hop between sites 0 and 11.

III.5 E. Semiclassical Numerics

While exact diagonalization and Matrix Product State techniques can give us access to the full quantum dynamics of the Hamiltonian (1), these techniques limit us to relatively small system sizes. To examine information spreading in larger system sizes while retaining access to strong interactions and chaos, we must resort to a semiclassical treatment by considering the limit where the spin lengths $S\rightarrow\infty$ Polkovnikov (2010). In this limit, each spin operator $\mathbf{S}_{i}$ can be understood as a classical angular momentum vector obeying the Poisson bracket algebra:

[TABLE]

and whose equations of motion are given by

[TABLE]

Defining unit-length spin vectors $\mathbf{x}_{i}=\mathbf{S}_{i}/S$ and taking the limit $S\rightarrow\infty$, we obtain classical equations of motion of the form:

[TABLE]

where $f$ is a quadratic function of the spin vectors $\mathbf{x}_{i}$ that depends on the coupling graph $J(i-j)$ but is independent of $S$.

To mimic an infinite-temperature quantum state consisting of exactly $N/2$ magnons, we prepare the spins along random directions in the $x-y$ plane. (Preparing spins with additional out-of-plane fluctuations in the $z$-direction, or with uniformly random directions, has no significant effect on the dynamics.) We then estimate the classical sensitivity $C_{\mathrm{cl}}(t)$ by perturbing a single spin at site $i$ by a small rotation $\phi_{i}$ about the $z$ axis, and measure the resulting change $\Delta\mathbf{x}_{j}(t)$ in spin $j$ at a later time $t$.

III.6 F. Simulating OTOCs with Matrix Product Operators

Here we describe how OTOCs are calculated using tensor network methods. We perform simulations with Matrix Product States and Operators (MPS/MPO) Schollwoeck (2011), which a regularly used to successfully simulate quantum many-body dynamics in one dimensional systems García-Ripoll (2006); Paeckel et al. (2019). Although here we have only made use of MPS and MPO techniques, in future work it might be beneficial to try other generalized forms of tensor networks for treating these sparse models, such as the Tree Tensor Network Shi et al. (2006) or the Multi-scale Entanglement Renormalization Ansatz Vidal (2008). These may account more naturally for the geometry of long-range sparse interactions between spins.

In our calculations, we consider the symmetry sector of the Hilbert space of $N$ spins with a fixed number of magnons, $n$. The initial state of the system $\rho_{0}$ is the infinite-temperature state within this symmetry sector, and may be written as $\rho_{0}=P_{n}/{\mathcal{Z}}$, where the $P_{n}$ is the projector on states with $n$ magnons and ${\mathcal{Z}}$ is the partition function, which in the case of infinite temperature is equal to the number of permutations of $n$ magnons on $N$ sites.

We construct the projector $P_{n}$ as an MPO by writing it as a product of $n\times n$ matrices,

[TABLE]

where

[TABLE]

are constructed out of operators $p_{\downarrow}={\left|\downarrow\right\rangle}{\left\langle\downarrow\right|}$ and $p_{\uparrow}={\left|\uparrow\right\rangle}{\left\langle\uparrow\right|}$ acting on local spins.

Having constructed this infinite temperature state, the next step is to time evolve the operator $S_{j}^{z}(t)$ from Eq. (5) in the Heisenberg picture. In order to do this, we combine the local indices of the operator, thus vectorizing it so that it can be propagated in time using the same techniques we apply for MPS. We find that using a TDVP Haegeman et al. (2011); Koffel et al. (2012); Haegeman et al. (2013, 2016) provides a good balance between speed and accuracy in this calculation (although for some purposes, using the Runge-Kutta method García-Ripoll (2006) could yield higher accuracy). After the time evolution, we separate out the local indices of $S_{j}^{z}(t)$ and contract it with $S_{i}^{z}(0)$ and $\rho_{0}$ to obtain $C(i,j;t)$ in Eq. (5).