Epidemic dynamics in open quantum spin systems

TL;DR

This paper investigates how epidemic-like spreading dynamics manifest in open quantum spin systems, revealing different phase transitions depending on classical or quantum effects, with implications for Rydberg atom experiments.

Contribution

It introduces a quantum model for epidemic dynamics in open many-body systems and analyzes phase transitions influenced by classical and quantum effects.

Findings

Different non-equilibrium phase transitions depending on classical or quantum dominance

Potential observability of phenomena in Rydberg gases with long-range interactions

Identification of conditions for epidemic spreading in quantum regimes

Abstract

We explore the non-equilibrium evolution and stationary states of an open many-body system which displays epidemic spreading dynamics in a classical and a quantum regime. Our study is motivated by recent experiments conducted in strongly interacting gases of highly excited Rydberg atoms where the facilitated excitation of Rydberg states competes with radiative decay. These systems approximately implement open quantum versions of models for population dynamics or disease spreading where species can be in a healthy, infected or immune state. We show that in a two-dimensional lattice, depending on the dominance of either classical or quantum effects, the system may display a different kind of non-equilibrium phase transition. We moreover discuss the observability of our findings in laser driven Rydberg gases with particular focus on the role of long-range interactions.

Figures

Equations

$${\Pi }_k=\sum _{l\in {\Lambda }_k}{r}_l\prod _{m\in {\Lambda }_k\backslash \left\{l\right\}}(1-r_m)=\sum _{l\in {\Lambda }_k}{r}_l+\dots ,$$

$${\partial }_t\mu =\sum _k\left[\alpha {\Pi }_k\left(\mathcal{L}({\sigma }_k^{+})+\mathcal{L}({\sigma }_k^{-})\right)+\kappa \mathcal{L}\left({|n⟩}_k⟨r|\right)\right]\mu ,$$

$${\partial }_{t}G$$

$${\partial }_{t}R$$

$${\Gamma }_k^{-1}\approx 1+R^{2\beta }{\left[1-\sum _{l\in {\Lambda }_k}{r}_l-\eta \text{sideset}{}^{\prime }{\sum }_m\frac{r_m}{|{\mathbf{x}}_k-{\mathbf{x}}_m{|}^{\beta }}\right]}^2,$$

$${\partial }_t{G}_k$$

$${\partial }_t{N}_k$$

$$g_k={|g⟩}_k⟨g|,r_k={|r⟩}_k⟨r|,n_k={|n⟩}_k⟨n|,{\sigma }_k^{+}={|r⟩}_k⟨g|,{\sigma }_k^{-}={|g⟩}_k⟨r|,{\sigma }_k^x={\sigma }_k^{+}+{\sigma }_k^{-}\text{ and }{\sigma }_k^y=-i\left({\sigma }_k^{+}-{\sigma }_k^{-}\right).$$

$$H=\Omega \sum _k{\Pi }_k{\sigma }_k^x,$$

$${\Pi }_k=\sum _{l\in {\Lambda }_k}{r}_l\prod _{m\in {\Lambda }_k\backslash \left\{l\right\}}(g_m+n_m)=\sum _{l\in {\Lambda }_k}{r}_l\prod _{m\in {\Lambda }_k\backslash \left\{l\right\}}(1-r_m)$$

$$J_k^{(\mathrm{imm})}=\sqrt{\kappa }{|n⟩}_k⟨r|\text{ and }{J}_k^{(\mathrm{deph})}=\sqrt{\gamma }{|r⟩}_k⟨r|,$$

$${\partial }_t\rho =-i\left[H,\rho \right]+\sum _k\left[\mathcal{L}(J_k^{(\mathrm{imm})})+\mathcal{L}(J_k^{(\mathrm{deph})})\right]\rho ,$$

$$\mathcal{L}(J)\rho =J\rho J^{†}-\frac{1}{2}\left\{J^{†}J,\rho \right\}.$$

$${\partial }_t\mu =\sum _k\left[\alpha \left(\mathcal{L}({\Pi }_k{\sigma }_k^{+})+\mathcal{L}({\Pi }_k{\sigma }_k^{-})\right)+\kappa \mathcal{L}\left({|n⟩}_k⟨r|\right)\right]\mu ,$$

$${\partial }_{t}Q=\sum _k\left[\alpha \left({\mathcal{L}}^{\ast }({\Pi }_k{\sigma }_k^{+})+{\mathcal{L}}^{\ast }({\Pi }_k{\sigma }_k^{-})\right)+\kappa {\mathcal{L}}^{\ast }({|n⟩}_k⟨r|)\right]Q,$$

$${\mathcal{L}}^{\ast }(J)Q=J^{†}QJ-\frac{1}{2}\left\{J^{†}J,Q\right\}.$$

$${\partial }_{t}Q=\sum _k\left[\alpha \sum _{l\in {\Lambda }_k}\left({\mathcal{L}}^{\ast }(r_l{\sigma }_k^{+})+{\mathcal{L}}^{\ast }(r_l{\sigma }_k^{-})\right)+\kappa {\mathcal{L}}^{\ast }({|n⟩}_k⟨r|)\right]Q.$$

$${\partial }_t{r}_m$$

$${\partial }_t{n}_m$$

$${\partial }_t{r}_m$$

$${\partial }_t{n}_m$$

$${\partial }_t{r}_m$$

$${\partial }_t{n}_m$$

$${\partial }_t{r}_m$$

$${\partial }_t{g}_m$$

$${\partial }_t{n}_m$$

$${\partial }_t{R}_m$$

$${\partial }_t{G}_m$$

$${\partial }_t{N}_m$$

$${\partial }_{t}R$$

$${\partial }_{t}G$$

$${\partial }_{t}N$$

$${\partial }_{t}O=i\left[H,O\right]+\kappa \sum _k{\mathcal{L}}^{\ast }({|n⟩}_k⟨r|)O$$

$${\partial }_t{r}_m$$

$${\partial }_t{n}_m$$

$${\partial }_t{\sigma }_m^{-}$$

Full text

Introduction — Cold atoms and ions are versatile platforms for the exploration of non-equilibrium physics. Recent examples include studies on creation and dynamics of quasi-particles Schirotzek et al. (2009); Jurcevic et al. (2014), spreading of entanglement and correlations Islam et al. (2015); Fukuhara et al. (2015); Kaufman et al. (2016), as well as many-body localization in disordered systems Choi et al. (2016); Schreiber et al. (2015); Smith et al. (2016). In particular, so-called Rydberg gases with their strong long-range interactions Saffman et al. (2010); Löw et al. (2012) permit the exploration of open and closed many-body physics Schauss et al. (2012); Schempp et al. (2014); Urvoy et al. (2015); Weber et al. (2015); Helmrich et al. (2016), with recent experiments probing non-equilibrium dynamics Helmrich et al. (2016); Valado et al. (2016); Letscher et al. (2017); Goldschmidt et al. (2016), phase transitions Carr et al. (2013); Malossi et al. (2014); Gutierrez et al. (2016); Ding et al. (2016) and disorder-induced localization phenomena Marcuzzi et al. (2017).

An intriguing aspect is that in Rydberg gases one can achieve control over the relative strength of quantum fluctuations and classical noise Löw et al. (2012). This permits the exploration of dynamical phenomena in settings which can be regarded as quantum generalizations of classical non-equilibrium systems Hinrichsen (2000); Ódor (2004). A recent example is a quantum version of the so-called contact process Marcuzzi et al. (2016); Buchhold et al. (2017), a simple stochastic model for population dynamics featuring a non-equilibrium phase transition Henkel et al. (2008) whose character changes drastically when moving from a purely-classical to a quantum regime.

In this work we shed light on the collective dynamics of an open quantum system generalization of a general epidemic process (GEP) Grassberger (1983), belonging to the dynamic percolation universality class Cardy (1983); Cardy and Grassberger (1985); Janssen (1985); Henkel et al. (2008). In a Rydberg system, a similar dynamics can be expected by considering atoms with three relevant states, which can be labeled as “healthy”, “infected” and “immune”, where infected sites have the ability to infect their healthy neighbors, or heal and become immune. A scenario similar to this has been recently realized and studied experimentally in Helmrich et al. (2016), where a connection to the GEP was conjectured. The scope of this work is not to propose a quantum simulation protocol for GEP, but to demonstrate that, under dominant classical noise, the system follows the same phenomenology as the GEP and undergoes a continuous transition between two phases: one where the contagion starting from an initially infected site is unable to percolate throughout the system, and one where the initial infected site triggers a self-sustaining wavefront (an “outbreak”) which covers an extensive portion of the system and leaves behind a trail of immunized sites. In the quantum regime, a mean-field treatment suggests that the density of immune sites displays a sequence of jumps resulting from the presence of multiple wavefronts.

Model — We consider atoms with three internal states: a ground state $|g\rangle$ (healthy), a Rydberg state $|r\rangle$ (infected) and a second stable state $|n\rangle$ (immune). These states are coupled as depicted in Fig. 1a: $|g\rangle$ is excited to $|r\rangle$ via a laser with Rabi frequency $\Omega$ and detuning $\Delta$ and the state $|r\rangle$ decays radiatively into $|n\rangle$ at rate $\kappa$. Note that this implicitly assumes that the decay from $|r\rangle$ does not proceed via long-lived intermediate states. Radiative decay from $|r\rangle$ to $|g\rangle$ is neglected for simplicity (see SM for details). The atoms are placed on a two-dimensional square lattice with $L$ sites and spacing $a$, one per site (see e.g. experiments in Refs. Barredo et al. (2015); Schauß et al. (2015)).

Collective behavior emerges when the probability for an atom to undergo the transition $|g\rangle\rightarrow|r\rangle$ (infection of a healthy site) depends on the state of its neighbors. For Rydberg atoms this is achieved by enforcing the so-called “facilitation” (or “anti-blockade”) condition Ates et al. (2007a); Amthor et al. (2010); Gärttner et al. (2013); Schönleber et al. (2014); Lesanovsky and Garrahan (2014); Urvoy et al. (2015); Valado et al. (2016). Here, the detuning $\Delta$ of the excitation laser is set to compensate the interaction $V_{\mathrm{NN}}$ between neighboring atoms, which makes the transition $|g\rangle\rightarrow|r\rangle$ resonant provided a neighbor is already in state $|r\rangle$ (Fig. 1(a)). This situation has already been explored in a two-level setup Schempp et al. (2014); Urvoy et al. (2015); Valado et al. (2016); Thaicharoen et al. (2015), and very recently also in the considered three-state setting Helmrich et al. (2016). We remark that it is crucial for infection to only occur locally, i.e. in the neighborhood of an already infected site. In our case, it is therefore important that the interactions decay sufficiently fast with the distance. A variety of different behaviors are known to emerge in systems where they do not Bouchet and Barre (2005); Chatterjee et al. (2013); Fischer et al. (2015); Fasshauer and Lode (2016).

We now consider a minimal model of the resulting many-body dynamics in which the density matrix $\rho$ of the system evolves under a Lindblad master equation Breuer and Petruccione (2002) $\partial_{t}\rho=-i\left[H,\rho\right]+\sum^{L}_{k=1}\left[\kappa\mathcal{L}\left(\left|n\right>_{k}\!\left<r\right|\right)\rho+\gamma\mathcal{L}\left(r_{k}\right)\rho\right]$, whose terms are sketched in Fig. 1(b). The coherent evolution is governed by the Hamiltonian $H=\Omega\sum_{k}\Pi_{k}\sigma_{k}^{x}$ where $\sigma_{k}^{x}=\left|g\right>_{k}\!\left<r\right|+\left|r\right>_{k}\!\left<g\right|=\sigma^{-}_{k}+\sigma^{+}_{k}$, and $\Pi_{k}$ is a projector onto the subspace in which exactly one of the nearest neighbors of $k$ is in state $\left|r\right>$. Denoting the set of nearest neighbors of $k$ by $\Lambda_{k}$, it reads

[TABLE]

with $r_{k}=\left|r\right>_{k}\!\left<r\right|$, $g_{k}=\left|g\right>_{k}\!\left<g\right|$ and $n_{k}=\left|n\right>_{k}\!\left<n\right|$. The dots denote higher-order terms in the operators $r_{k}$. This projector constrains infection to occur only in the neighborhood of a single infected site. The dissipative dynamics is described by Lindblad terms $\mathcal{L}\left(J\right)\rho=J\rho J^{\dagger}-\{J^{\dagger}J,\rho\}/2$. The first describes decay from $|r\rangle$ to $|n\rangle$ at rate $\kappa$, the second dephasing of quantum coherences at rate $\gamma$. Controlling the dephasing strength, which is achievable by modifying the excitation laser linewidth or the temperature of the atoms, allows switching between classical and quantum regimes Löw et al. (2012); Levi et al. (2016). In the following, the initial state is always a single atom in state $\left|r\right\rangle$ (infected) in the center of the lattice and all the others in state $\left|g\right\rangle$ (healthy).

Classical regime — We first consider the regime of strong dephasing $\gamma\gg\Omega$. Here an effective dynamics can be defined for the diagonal of the density matrix $\mu_{ij}=\delta_{ij}\rho_{ii}$ in the $\left|r,g,n\right\rangle$ basis and the corresponding (classical) master equation reads Degenfeld-Schonburg and Hartmann (2014); Marcuzzi et al. (2014)

[TABLE]

with $\alpha=4\Omega^{2}/\gamma$. This means that atoms undergo incoherent state changes from $|g\rangle$ to $|r\rangle$ and vice versa, with a rate conditioned on their local neighborhood. Furthermore, decay from $|r\rangle$ to $|n\rangle$ is possible. The process (2) is similar to a GEP but differs from it by (i) the presence of facilitated transitions $\left|r\right\rangle\to\left|g\right\rangle$ and (ii) for the stricter constraint that a single neighboring infected site is required for facilitation, whereas in the GEP the infection rate is proportional to the number of infected neighbors.

In a GEP, for fixed $\kappa$, if the facilitation rate lies below its critical value, the initial infection is unable to propagate and the density of immunes $N=\sum_{k}\langle n_{k}\rangle/L$ ($L$ being the number of sites) vanishes in the thermodynamic limit. Conversely, above the critical point there is a finite probability for the infection to percolate (see Fig. 1(c)), propagating as a single travelling wavefront and leaving behind a finite fraction of immune sites $N>0$ Grassberger (1983) (i.e. a single outbreak takes place).

We start from a uniform mean-field approximation where we neglect the higher order terms in (1), effectively relaxing (ii). Introducing the quantities $R=\sum_{k}\langle r_{k}\rangle/L$ and $G=\sum_{k}\langle g_{k}\rangle/L$, this yields

[TABLE]

Analogously to what is found for the GEP Grassberger (1983), these equations feature a constant of motion $\partial_{t}(\log(R-G+\kappa/8\alpha)+8N\alpha/\kappa)=0$ which permits the determination of the stationary phase diagram for different initial conditions, shown in Fig. 2(a) SM . In a uniform approximation, the closest initial condition to the one we start from is a vanishingly small density of infections ($R(t=0)=\epsilon\to 0^{+}$) and in this limit the two phases (which are illustrated in Fig. 1(c)) can be clearly identified, separated by a critical point at $\alpha_{c}=\kappa/4$. The non-uniform mean-field dynamics is shown in Fig. 2(b) and highlights the absence/presence of an outbreak in the two phases.

We have then performed continuous-time Monte Carlo simulations of the classical master equation (2). In Fig. 2(c) we show the stationary density $N_{\mathrm{SS}}=\lim_{t\to\infty}N(t)$ for different system sizes and observe a sharpening crossover from a vanishing to a finite-valued phase when increasing $L$. The dynamics shows the expected GEP behavior: for $\alpha<\alpha_{c}\approx 1.71\,\kappa$ the process fails to percolate and no outbreak is produced. In the supercritical phase, instead, there is a finite probability of a single outbreak immunizing a macroscopic portion of the system. A finite-size scaling analysis performed with the tabulated critical exponents for the GEP provides an excellent collapse of the curves SM , establishing the connection on firm grounds. The full constraint given by (1) leads to the same scaling behavior near the critical point. This could be expected since the differences between our process and the GEP are of higher order in the density of infected sites, which vanishes close to the critical point. More rigorously, one can show that these terms are irrelevant in a renormalization-group sense SM .

Rydberg gases — We now discuss the observability of this physics in Rydberg gases interacting with a power-law potential $V(r)=C_{\beta}/r^{\beta}$ (here we consider the van der Waals case $\beta=6$) under anti-blockade conditions. As shown in Refs. Ates et al. (2007b); Lesanovsky and Garrahan (2014); Marcuzzi et al. (2014, 2015), in this case the constraint $\Pi_{k}$ is replaced by the operator $\Gamma_{k}=1/(1+R^{2\beta}(1-\sum_{m\neq k}\frac{r_{m}}{|{\bf{x}}_{k}-{\bf{x}}_{m}|^{\beta}})^{2}),$ where $R=(2C_{\beta}/a^{\beta}\gamma)^{1/\beta}$ is the so-called dissipative blockade radius and ${\bf{x}}_{k}$ denotes the position of the $k$-th atom. The lattice spacing $a$ is taken as the facilitation distance: if just one nearest neighbor of site $k$ is infected, the rate is maximized, i.e. $\Gamma_{k}=1$, whereas other configurations lead to its suppression. In this context, two main effects alter the physics of the previous model: the possibility of unfacilitated infection ($\left|g\right\rangle\leftrightarrow\left|r\right\rangle$ in the absence of infected neighbors) and the fact that the interactions extend beyond nearest neighbors Gutierrez et al. (2016). To gain insight into their role, we exploit the rapid decay of the tails of $V(r)$ to truncate the interaction beyond a distance of two lattice sites,

[TABLE]

where the primed sum runs only over $m$ such that $1<|{\bf{x}}_{k}-{\bf{x}}_{m}|\leq 2$. Here we have introduced a parameter $\eta$, which allows us to control the strength of the ‘long-range’ part. While this parametrization is used here for convenience, in practice potential shaping techniques can be applied to modify and possibly suppress the potential tails, see e.g. Refs. Büchler et al. (2007); Sevinçli and Pohl (2014); Marcuzzi et al. (2015).

First, we consider $\eta=0$ (nearest-neighbors interactions only), where $\Gamma_{k}$ is well approximated by the constraint $\Pi_{k}$ [Eq. (1)], provided that $R>1$. Unfacilitated (spontaneous) infection can now occur at a rate $\alpha\Gamma_{\mathrm{spont}}=\alpha/(1+R^{2\beta})>0$. Albeit rare for $R\gg 1$, these processes dramatically alter the stationary-state properties of the system, invariably leading to $N_{\mathrm{SS}}=1$. In the renormalization-group language, the spontaneous processes constitute a relevant perturbation. Nevertheless, for sufficiently large $R$ a timescale separation occurs: the outbreak follows the phenomenology observed for the idealized case up to times on the order of $\kappa t\sim(1+R^{2\beta})/(L\alpha/\kappa)$, which is an underestimate of the mean waiting time of the Poisson process producing spontaneous infection. This is illustrated in Fig. 3(a), where we show the stationary density of immune states $N_{\mathrm{SS}}$ of the idealized process (rate function $\alpha\Pi_{k}$), together with the density $N(t)$ at time $\kappa t=100$ for the Rydberg rates with $R=5$, $\beta=6$ and $\eta=0$. These curves display remarkable agreement showing that this phase transition in fact underlies the transient Rydberg dynamics. In the inset of Fig. 3(a) we moreover show how the characteristic sigmoidal profile remains stable for a long period of time. In Fig. 3(b) we show that outbreaks are visible for some time, but spontaneous infection leads to a rising “background” which eventually overcomes the epidemic process and immunizes the entire system.

Long-range interactions ($\eta>0$) counteract the formation of large clusters of infected and immune sites. For instance, an isolated atom in state $|r\rangle$ can infect, say, the site directly below it at the maximal rate, since the nearest-neighbor interaction compensates the laser detuning. To infect a third one at the right of the former, however, one has to include the additional shift due to next-nearest-neighbor interactions, which brings the atomic transition off resonance, hindering the propagation. We illustrate this by showing in Fig. 3(a) the density of immune states at time $\kappa t=100$ for $R=2.5$ and different values of $\eta$ (all red symbols). As $\eta$ is increased, it becomes more difficult for the process to spread. This however is a lattice effect; in dense atomic clouds the transition may reappear. This is suggested by recent experimental works Urvoy et al. (2015); Valado et al. (2016); Helmrich et al. (2016); Gutierrez et al. (2016) that reveal collective phenomena in the presence of thermal motion.

Quantum regime — We now set the dephasing rate $\gamma=0$ and study what we refer to as the “quantum case”. We describe the dynamics via the mean-field equations

[TABLE]

Here, $R_{k}\,(G_{k},N_{k})=\left\langle r_{k}\,(g_{k},n_{k})\right\rangle$ and $\Sigma_{k}=\langle\sigma^{y}_{k}\rangle$ with $\sigma_{k}^{y}=i\left|g\right>_{k}\!\left<r\right|-i\left|r\right>_{k}\!\left<g\right|$. Additionally, $\overline{R}_{k}=\sum_{l\in\Lambda_{k}}\left\langle r_{l}\right\rangle$ and we assume that initially no coherences are present. Meanfield will yield qualitatively reasonable predictions unless long-range correlations develop. As we shall show, the transition becomes discontinuous in the quantum regime and is therefore not associated to a divergent correlation length. Our analysis cannot capture the role of spatial dimensionality, but should identify the correct qualitative behavior in sufficiently high dimension, see Buchhold et al. (2017). A more detailed analysis could be achieved exploiting cluster mean-field methods Jin et al. (2016).

In Fig. 4(a) we show the stationary state density $N_{\mathrm{SS}}$ resulting from Eqs. (5). A striking difference with respect to the classical case is the appearance of an oscillating behavior as a function of the driving parameter $\Omega/\kappa$. From our numerical analysis, it appears that the peaks become sharper as the size of the system is increased, with their positions remaining approximately fixed. This suggests that in the thermodynamic limit a sequence of discontinuous jumps will form at fixed values $\widetilde{\Omega}_{j}/\kappa$, $(j=1,2,\ldots)$. Note that a standard numerical analysis of the problem (exact diagonalization, quantum-jump Monte Carlo) is forbiddingly complex, due to the exponential scaling of the state space dimension with the number of atoms.

These jumps feature an intriguing dynamic counterpart: in the quantum case more than one outbreak can occur. As shown in Fig. 4(b), if $\Omega/\kappa$ is too small (first column), no outbreak takes place, as in the classical case. For $\Omega\gtrsim{\widetilde{\Omega}}_{1}\approx 0.83\kappa$, instead, a single outbreak leaves behind an approximately uniform density of immunes (second column). Increasing $\Omega/\kappa$ further, this residual density decreases until, at a second threshold value ${\widetilde{\Omega}}_{2}\approx 1.643\kappa$, a second outbreak is generated, which causes the final density $N_{\rm{SS}}$ to jump to a higher value. Every new jump in $N_{\rm{SS}}$ appears to be associated to a new outbreak in the dynamics.

This repeating pattern allows us, by analyzing the first jump, to make predictions on the subsequent ones. To this end, we make two simplifying assumptions compatible with the numerically observed behavior: (I) every outbreak leaves behind a uniform density of immunes and (II) stems from the center of the lattice. We focus now on a neighbourhood of the center after the $(j-1)$-th outbreak, at some time $t_{j}$. By (I) there is an immune density $N_{k}(t_{j})=N^{(j)}>0$; by (II), $R_{k}$, $G_{k}$ and $\Sigma_{k}$ at time $t_{j}$ correspond to their initial conditions rescaled by $(1-N^{(j)})$. Therefore, the facilitation rate is bounded by $\Omega\overline{R}_{k}\leq 4\Omega(1-N^{(j)})$; as a first approximation, the last factor can be reabsorbed by $\Omega\to\Omega^{(j)}=\Omega(1-N^{(j)})$. In other words, the process after $t_{j}$ proceeds like at $t=0$ (no immunes present), but with a modified frequency $\Omega^{(j)}<\Omega$.

Hence, if $\Omega^{(j)}<\widetilde{\Omega}_{1}$ the process stops at $N_{\rm{SS}}=N^{(j)}$, up to subextensive additions, meaning that $\Omega<\widetilde{\Omega}_{j}$. If $\widetilde{\Omega}_{1}<\Omega^{(j)}<\widetilde{\Omega}_{2}$, instead, a $j$-th outbreak will be produced, but a $(j+1)$-th will not take place, corresponding to $\widetilde{\Omega}_{j}<\Omega<\widetilde{\Omega}_{j+1}$, and so forth. All the processes sharing the same reference frequency $\Omega^{(j)}$ are thus equivalent to the same reference process occurring in the absence of immunes and their stationary points will lie on a curve $\Omega(1-N_{\rm{SS}})=const$. The extremal curves in this set, passing through the top and bottom of the first jump, are displayed in Fig. 4(a) and bound well the data. With an ansatz $N_{\rm{SS}}\propto 1/\Omega$ for the decrease between jumps, one can formulate a more detailed prediction out of the same considerations, represented by the blue line (see SM for details).

Conclusions — We have analyzed a simple model for epidemic spreading in an open quantum system, which has been inspired by recent experimental work Helmrich et al. (2016), and investigated its connection with the so-called general epidemic process Grassberger (1983). In the presence of strong dephasing, the process has a direct relation to the GEP, displaying a continuos transition in the same universality class. In the quantum limit, instead, an intriguingly different physics emerges featuring a sequence of discontinuous jumps. This susprising behavior warrants further theoretical and experimental investigation.

S1 Derivation of the equations of motion in the open quantum system formalism

We consider an arrangement of quantum 3-level systems on a square lattice. We call the three states $|g\rangle_{k}$, $|r\rangle_{k}$ and $|n\rangle_{k}$, which correspond to the healthy, infected and immune states discussed in the main text, and $k$ denotes the lattice index. For later convenience, we introduce the one-site operators

[TABLE]

As in the main text, we shall also denote by $\Lambda_{k}$ the set of all lattice indices which are nearest neighbors with site $k$.

The coherent part of the dynamics is generated by a Hamiltonian

[TABLE]

where

[TABLE]

is the projector over all possible configurations in which a single neighbor of site $k$ is in state $\left|r\right\rangle$ (infected) and all the others are in superpositions of $\left|g\right\rangle$ (healthy) and $\left|n\right\rangle$ (immune).

The dissipative part is instead given in terms of jump operators

[TABLE]

which describe immunization ($\left|r\right\rangle\to\left|n\right\rangle$) and dephasing (decay of coherence between states $\left|r\right\rangle$ and the remaining ones) at site $k$.

The state $\rho$ of the system evolves under the Lindblad equation

[TABLE]

where

[TABLE]

S1.1 Effectively classical equations of motion in the large dephasing limit

We consider here the large dephasing limit $\gamma\gg\Omega$. This condition induces a separation between the timescales on which coherence is produced ($\propto\Omega^{-1}$) and destroyed ($\propto\gamma^{-1}$) and allows to perform an adiabatic elimination of the coherent terms of the density matrix $\rho$ (see Refs. Degenfeld-Schonburg and Hartmann (2014); Marcuzzi et al. (2014) and the projective Nakajima-Zwanzig method described in Breuer and Petruccione (2002) and references therein). In simpler terms, if the density matrix is written in the basis spanned by the “classical” states $\otimes_{k}\left|g,r,n\right\rangle_{k}$ this allows to write an effective reduced master equation involving just the diagonal part $\mu$, whose elements are simply $\mu_{ij}=\delta_{ij}\rho_{ii}$. This effective equation (see, e.g., Lesanovsky and Garrahan (2013, 2014); Everest et al. (2016)) reads, to leading order in $\Omega/\gamma$,

[TABLE]

with $\alpha=4\Omega^{2}/\gamma$. This equation is equal to Eq. (2) in the main text, since $\Pi_{k}$ is a diagonal projector which only acts on the neighbors of $k$ and therefore $\Pi_{k}^{2}=\Pi_{k}$ and it commutes with all $k$-th-site-local operators (such as $\sigma_{k}^{\pm}$) and $\mu$. This is a classical master equation for the classical probabilities $\mu_{ii}$ and describes a stochastic process with three main mechanisms, which we are going to refer to, in the order in which they appear above, as “classical infection”, “classical facilitation down” and “immunization”. The second process is absent in the general epidemic process and corresponds to the possibility for an infected site to heal a nearby infected site without immunising it at a rate $\alpha$. The corresponding evolution in the Heisenberg picture for a quantity $Q$ (corresponding to a diagonal operator in the classical basis) reads now

[TABLE]

where $\mathcal{L}^{\ast}$ is the adjoint dissipator

[TABLE]

In the following, since higher-order terms in the immune density are not expected to modify the physics of the transitions studied in the main text, we shall replace $\Pi_{k}$ by a softer constraint, which makes the rate proportional to the number of infected neighbors a site has. The corresponding simplified equation reads

[TABLE]

Since this is a linear equation, one can study it term by term. The relevant evolution equations for the one-site operators are

• •

classical infection

[TABLE]

• •

classical facilitation down

[TABLE]

• •

immunization

[TABLE]

All the remaining terms can be obtained via the local probability conservation $r_{m}+g_{m}+n_{m}=1$ (the probability of a site being in any state is $1$). Employing the shorthand $\overline{r}_{m}=\sum_{l\epsilon\Lambda_{m}}r_{l}$ we arrive at

[TABLE]

We now introduce the notation $N_{m}=\left\langle n_{m}\right\rangle$, $R_{m}=\left\langle r_{m}\right\rangle$, $G_{m}=\left\langle g_{m}\right\rangle$ and $\overline{R}_{m}=\left\langle\overline{r}_{m}\right\rangle$ for the expectation value of these observables and perform a mean-field approximation. This means we factorize all correlation functions acting on different sites, i.e., $\left\langle Q_{k}Q_{m}\right\rangle\to\left\langle Q_{k}\right\rangle\left\langle Q_{m}\right\rangle$ for $k\neq m$ and for any one-site-local quantity $Q$. The equations then read

[TABLE]

If we additionally assume translational invariance and introduce the notation $N=\sum_{k}\left\langle n_{k}\right\rangle/L$, $R=\sum_{k}\left\langle r_{k}\right\rangle/L$, $G=\sum_{k}\left\langle g_{k}\right\rangle/L$ and the coordination number $z$ (number of nearest neighbors per site, $4$ for a square lattice) we arrive at

[TABLE]

S1.2 Equations of motion in the quantum case

For $\gamma=0$, the derivation of the equations of motion for the one-site observables proceeds along the same lines, with the adjoint Lindblad equation now reading

[TABLE]

for any observable $O$ and $\mathcal{L}^{\ast}$ as in (S9). For simplicity, we consider here the softer constraint $\Pi_{k}\to\overline{r}_{k}$ as well. We can now write again the various contributions term by term, which read

• •

immunization

[TABLE]

[TABLE]

[TABLE]

• •

quantum facilitation

[TABLE]

[TABLE]

The equations for $\sigma^{+}_{m}$ can be derived by taking the hermitian conjugate of the equations for $\sigma^{-}_{m}$, and those for $\sigma_{m}^{x/y}$ simply by linear combination according to the definitions. Combining all terms and defining $\overline{s}^{x}_{m}=\sum_{l\in\Lambda_{m}}\sigma^{x}_{l}$, in this case the equations read

[TABLE]

The mean-field approximation yields in this case the equations

[TABLE]

where we additionally defined $\overline{S}_{m}^{x}=\left\langle\overline{s}^{x}_{m}\right\rangle$ and $\Sigma^{x/y}_{m}=\left\langle\sigma^{x/y}_{m}\right\rangle$. We see now that, if $\Sigma_{m}^{x}=0~{}\forall m$, then $\partial_{t}\Sigma_{m}^{x}=0~{}\forall m$ and the $x$ component of the coherences remains null in the mean-field dynamics. Restricting to initial conditions with vanishing coherences, one can thus neglect the last equation altogether and work with the reduced set

[TABLE]

and make the substitution $\Sigma^{y}_{m}\to\Sigma_{m}$, which reproduces Eqs. (5) in the main text. Introducing the notation $\Sigma=\sum_{m}\left\langle\sigma_{m}\right\rangle/L$, the uniform case reads now

[TABLE]

S2 Properties of the homogeneous mean-field equations

S2.1 Classical equation limit

In the large dephasing limit $\gamma\gg\Omega$ Eqs. (S20) and (S21) close among themselves. The fact that $R=0$ is sufficient to annihilate both derivatives shows that, in principle, any value of $G$ is admitted in the stationary state, which reflects the fact that any configuration with only healthy and immune sites is absorbing for the stochastic process. The specific value of $G$ that the dynamics asymptotically reaches is then established by the initial condition, which hints at the presence of a constant of motion $F(G,R,N)$ ($\partial_{t}F(G,R,N)=0$), apart from the trivial one $G+R+N$, which carries the corresponding information. This is analogous to the considerations found in Grassberger (1983); Janssen (1985) for the deterministic equations. In our case, the equation for the difference

[TABLE]

can be recast in the form

[TABLE]

which identifies

[TABLE]

and can be implicitly integrated to yield

[TABLE]

where $R_{0}$, $G_{0}$ and $N_{0}$ denote the initial conditions. Exploiting the conservation of probability $R+G+N=1$, we can also write

[TABLE]

and, since $R=\partial_{t}N/\kappa$, after some rearrangements,

[TABLE]

Abbreviating $\xi=2z\alpha/\kappa$ and rescaling time by $\tau=\kappa t$ we find the more compact expression

[TABLE]

Note that this single equation for $N$ is not generally equivalent to Eqs. (S20), (S21), (S22): it only yields the same result as long as the same initial conditions $R_{0}$ and $G_{0}$ for $R$ and $G$ are taken. Finally, setting $N_{0}=0$, $R_{0}=\epsilon$, $G_{0}=1-\epsilon$ (initially a fraction $\epsilon$ of sites is in state $|r\rangle$, while the rest are in state $|g\rangle$), one finds

[TABLE]

In order to study the stationary properties (at initial conditions fixed as above), we now consider the long-time limit and set $\partial_{\tau}N=0$. The initial condition we assume in the main text is of a single infected site placed at the centre of the system at $t=0$, corresponding to a vanishing initial density ($\epsilon\to 0^{+}$) in the thermodynamic limit. We thus set $\epsilon=0$ as well in (S54). Note that the limits $t\to\infty$ and $\epsilon\to 0^{+}$ do not commute: starting the dynamics with $R_{0}=0$ will yield the trivial solution $R=0$, $G=G_{0}$, $N=N_{0}$. The stationary equation now reads

[TABLE]

where $N_{\mathrm{SS}}$ stands for the density of immunes in the stationary state. Since $\xi\geq 0$, the function on the r.h.s. is concave, whereas the one on the l.h.s. is just a line. Hence, at any value of $\xi$ there can be at most two intersections. One is always present at $N_{\mathrm{SS}}=0$. Evaluating the derivatives in [math], it is not difficult to show that the line is tangent to the curve for $\xi=2$. For $\xi>2$ a second point of intersection appears at some value $0<N_{\mathrm{SS}}<1$, while the solution $N_{\mathrm{SS}}=0$ becomes unstable. For $1<\xi<2$ this second point lies at negative values and is thus unphysical. For $\xi\leq 1$ there is only the vanishing solution. Hence, $\xi_{c}=2$ marks a critical point where the system switches from a non-percolating phase ($\xi<2$) where the process always stops at some finite time, and a percolating phase ($\xi>2$) where the process has a finite probability of surviving for arbitrary times.

For $\xi>2$ but close to the critical point, the stationary density of immune sites $N$ remains small and one can Taylor expand Eq. (S55) to find

[TABLE]

Thus, the mean-field prediction for the order parameter critical exponent is $\beta_{MF}=1$.

In Figure S1, we show the meanfield solution (S54) for different values of $\epsilon$, where one can see a crossover between the sub- and super-critical phases that approaches a continuous phase transition as $\epsilon\to 0^{+}$. These solutions have been found numerically, assuming the underlying lattice is square, and therefore $z=4$, which leads to $(\alpha/\kappa)_{c}=\frac{\xi_{c}}{8}=\frac{1}{4}$.

S2.2 Quantum regime

In the quantum regime $(\gamma=0)$ the equations of motion (S43)-(S46) are such that $R=\Sigma=0$ makes all derivatives vanish. As in the classical case, all states with no coherence nor component over the states $\left|r\right\rangle_{k}$ are valid stationary states of the system. The values reached by $G$ and $N$ in the long-time limit are hence once again dependent on the initial conditions. In this case, however, we have been unable to identify a constant of motion (apart from the trivial one $N+R+G=1$) to connect the initial and final values. We thus resort to a numerical analysis. For the initial conditions $G(0)=1-\epsilon$, $R(0)=\epsilon$ and $N(0)=\Sigma(0)=0$, and again taking $z=4$, we find the stationary state density $N_{\mathrm{SS}}$ displayed in Fig. S2. A closer inspection of these curves suggests that, despite the abrupt increases that look like discontinuities due to the numerical resolution, they do not feature any discontinuous behavior.

As pointed out in the main text, a clear difference with respect to the classical case is the appearance of an oscillating behavior of $N_{\rm{SS}}$ in $\Omega/\kappa$. For the uniform mean-field equations, the oscillation peaks tend to shift to larger values of $\Omega/\kappa$ as the initial density of infected sites, $R(0)=\epsilon$, is decreased. For $\epsilon=10^{-5}$ the system remains in a state with $N_{\mathrm{SS}}\approx 0$ over the whole parameter range shown. We have verified that, in the parameter range shown, the position of the first jump scales as $\widetilde{\Omega}_{1}\sim\kappa/(\sqrt{\epsilon})$ for small $\epsilon$. Extrapolating this behaviour to $\epsilon\to 0^{+}$ indicates that, in contrast to the classical case, the homogeneous meanfield equations do not lead to epidemic spreading for any finite $\Omega$ as long as the initial infected density is vanishingly small ($\epsilon\to 0^{+}$).

The reason for the discrepancy between the predictions obtained from the homogeneous meanfield equations and the space-dependent meanfield dynamics explored in the main text seems to be related to the ambiguity existing in the definition of the initial condition. While meant to reproduce the situation of a single infected site among $L$, in fact, the uniform mean-field equations only admit uniform initial conditions. Conceptually, the closest one to reality is one where $R(0)=1/L$ while $G(0)=1-1/L$. In the classical case, this can be interpreted as a uniform probability density of finding a single particle anywhere in the system. In the quantum case, instead, another interpretation in terms of quantum superpositions is possible: one could have, in fact, every site in state $\sqrt{L^{-1}}\,|r\rangle+\sqrt{1-L^{-1}}\,|g\rangle$. It is very much possible that, due to the different structure of the equations, the quantum process is unable to give rise to a self-sustaining infection in the limit $L\to\infty$ if $\Omega$ is kept finite. The numerical solution of the inhomogeneous equations (S39)-(S42) shows a very different behaviour (as discussed in the main text) when instead we set the real initial condition of having a single site in state $\left|r\right\rangle$ in the center and all the others in state $\left|g\right\rangle$.

S3 Critical behavior of the epidemic dynamics in the classical limit

In order to determine the universality class of the classical stochastic dynamics given by Eqs. (1)-(2) in the main text, we have performed a finite size scaling analysis for both the full constraint $\Pi_{k}=\sum_{l\in\Lambda_{k}}r_{l}\prod_{m\in\Lambda_{k}\backslash\left\{l\right\}}(g_{m}+n_{m})$ and its linear approximation ($\Pi_{k}\approx\sum_{l\in\Lambda_{k}}r_{l}$). Performing the same finite size scaling analysis of Ref. Argolo et al. (2011) we have first derived the critical value of the control parameter $\alpha_{c}(L)$ (setting $\kappa=1$ for simplicity) by locating the maximum of the numerical derivative of the order parameter $N_{\text{SS}}$ with respect to $\alpha$ or, in other words, the inflection point of the curve, where it switches from being convex to concave. In the thermodynamic limit $L\to\infty$, this point will approach the actual critical value $\alpha_{c}$, where a non-analiticity develops, signalling the presence of a second order phase transition. Close to this point, in the so-called critical scaling region, observables such as the stationary density of immunes $N_{\rm{SS}}$ will acquire a scaling form, i.e., under a rescaling of all distances by a factor $b$

[TABLE]

where $\Delta_{N}$ and $\Delta_{\alpha}$ are the scaling dimensions of $N_{\rm{SS}}$ and the control parameter $\alpha-\alpha_{c}$, respectively, and we recall that, the system being two dimensional, $L$ scales like a surface (i.e., $\propto b^{2}$). By fixing $b=1/\sqrt{L}$ we can then write

[TABLE]

with $f$ a universal (up to a multiplicative constant) scaling function. Using standard conventions, we relate the scaling dimensions to the static critical exponents $\beta$ (do not confuse this $\beta$ with the one of Eq. (6) in the main text) and $\nu$ (see e.g. Ref. Grassberger (1983)) according to

[TABLE]

which therefore implies that

[TABLE]

i.e., once rescaled by $L^{\frac{\beta}{2\nu}}$ and plotted as functions of $L^{\frac{\Delta_{\alpha}}{2}}(\alpha-\alpha_{c})$, all datasets obtained at different system sizes should collapse onto the same master curve $f$.

In order to verify this, we have used the critical exponents characterizing the general epidemic process, which are those of the dynamic isotropic percolation universality class Grassberger (1983) ($\beta=5/36$ and $\nu=4/3$). The numerical results are displayed in Figs. S3 (GEP constraint $\overline{R}_{k}$, i.e. facilitation rate proportional to the number of infected neighbors) and S4 (stricter constraint $\Pi_{k}$, i.e. facilitation rate only $>0$ when a single infected neighbor is present) for different system sizes, where one can observe a collapse for both types of dynamics, strongly indicating that they belong to the same universality class and that the latter is isotropic percolation for the stationary properties.

S3.1 Irrelevance of higher-order terms in the infected-site density

We exploit here the results known for the effective field-theoretical description of the GEP Cardy (1983); Janssen (1985); Cardy and Grassberger (1985) to verify that the new terms generated by (A) the classical facilitation down (S12) and by (B) the more stringent constraint $\overline{R}_{k}\to\Pi_{k}$ are irrelevant in a renormalization-group sense and thus are not expected to alter the universal properties of the transition. First, let us rewrite the classical uniform mean-field equations (S20)-(S22) as

[TABLE]

where now we distinguish

[TABLE]

Following Ref. Janssen (1985), we start from the deterministic equations and formally integrate the last two to find

[TABLE]

The implicit expression for $G(t)$ as a function of $R(t)$ can now be substituted in the first equation. The subsequent step would be to expand the resulting r.h.s. in powers of the infected density $R$, keeping only the lowest orders. We thus immediately notice that

[TABLE]

is equivalent to $\mathcal{F}^{(GEP)}$ up to higher orders, which would at most rescale the coupling constants by some amount. In particular, at lowest order

[TABLE]

which shows that, since $z>1$ and $G\geq 0$, the correction never makes the linear term vanish (i.e., it amounts to an irrelevant rescaling which can be reabsorbed in the coupling itself). For this reason, in the following we only consider $\mathcal{F}^{(GEP)}=1$, to avoid purposeless complications. This will also make it relatively easy to distinguish GEP terms (surviving when $c=0$) from new ones (vanishing for $c=0$). For later convenience, we also introduce the notation

[TABLE]

By Taylor expanding in $R$ we thus find

[TABLE]

One can now re-establish the dependence on spatial coordinates by switching the variables with continuous classical fields $R(t)\to R(\vec{x},t)$, $\mathcal{R}_{i}(t)\to\mathcal{R}_{i}(\vec{x},t)$ and furthermore introduce stochastic fluctuations in the form of a random Gaussian source $\xi{\vec{x},t}$ which obeys

[TABLE]

with $K$ a positive constant and where the proportionality of the variance to the amplitude of $R$ accounts for the fact that all configurations with $R\equiv 0$ are absorbing, and thus fluctuationless — i.e., whenever the system becomes depleted of infected individuals the stochastic dynamics halts.

The equation of motion for $R$ then takes the form of a Langevin equation

[TABLE]

where the space and time dependence is intended and only the first (diffusive) term in a derivative expansion is kept with an effective coupling $D$. We have introduced shorthands for the coupling constants, where the $u_{i}$s belong to the GEP, whereas the $c_{i}$s are “new”. For the reader’s convenience, the “non-GEP” terms have also been highlighted in red color.

Langevin equations can be generically mapped onto equivalent effective field theories via the so-called Martin-Siggia-Rose-de Dominicis-Janssen (MSRDJ) formalism Täuber (2014). Here we do not review it, but simply exploit the results found in Ref. Janssen (1985) to determine the engineering scaling dimension of the new terms. To this end, we should mention that the construction of the effective action in the MSRDJ formalism comes at the cost of introducing an additional field $\widetilde{R}$ (referred to as response field, as it encodes properties of response functions). In the MSRDJ action, the deterministic part of the Langevin equation (S73) (i.e., everything but $\xi$) appears multiplied by the response field. Therefore, the effective action in our case reads

[TABLE]

where the higher orders in $\widetilde{R}$ are generated by the properties of the noise $\xi$ and are therefore “GEP-like”. One typically fixes the scaling dimension of the coordinates to be $\left[\vec{x}\right]=-1$, which implies, due to the diffusive nature of the problem, that $\left[t\right]=-2$ for time (in other words, the dynamic critical exponent is $2$). Taking from Janssen (1985) the scaling dimensions of the fields

[TABLE]

we thus find

[TABLE]

we are now able to calculate the dimensions of all coupling simply by requiring each term of the action to have net dimension zero. As an example, we work out here $\left[u_{3}\right]$; we impose

[TABLE]

i.e., $\left[u_{3}\right]=3-d/2$, which is $>0$ for $d<6$. Therefore, this term is only relevant in spatial dimensions smaller than $d_{c}=6$, which represents the upper critical dimension for dynamic isotropic percolation. The other couplings’ dimensions can be calculated in an analogous fashion, yielding

[TABLE]

These values are all $<0$ at the upper critical dimension $d=d_{c}=6$, showing that the corresponding terms in the action are strictly irrelevant and only need to be accounted for if the coupling $u_{3}$ is fine-tuned to [math]. We have shown above, however, that this never happens in our equations, even if the more stringent constraint (B) is enforced.

We finally discuss the effect of an additional process which introduces recovery without immunization, i.e., a decay $\left|r\right\rangle\to\left|g\right\rangle$. Assuming it occurs uniformly throughout the system at a rate $\lambda$ the equations of motion get modified according to

[TABLE]

By defining $\mathcal{G}=G-\lambda/(z\alpha)$, however, we can recast them in the same form as before

[TABLE]

Hence, one can follow exactly the same steps and find a Langevin equation like (S73) up to the substitution $G_{0}\to\mathcal{G}_{0}$. The introduction of this process thereby amounts to an overall rescaling of the couplings of the effective action (S74):

[TABLE]

The vanishing of the quadratic term occurs at $\mathcal{G}_{0}=\kappa/(z\alpha)$, at which $u_{3}=z\alpha\kappa>0$, implying that the new process does not affect the fundamental properties of the theory, i.e., the leading non-linearity at the critical point is still the cubic one. Of course, the actual position of the critical point is shifted to larger values of $\alpha$, as one could expect considering that this process hinders the growth of a self-sustaining front of infection (an outbreak).

S4 Repeating pattern of $N_{\rm{SS}}$ in the quantum case

We give here the details on how the theoretical (blue) curve in Fig. 4(a) in the main text can be determined. We start from the inhomogeneous mean-field equations for the quantum case, which we reproduce below for the reader’s convenience:

[TABLE]

Firstly, by measuring energies in units of $\kappa$ and time in units of $\kappa^{-1}$ we can set, without loss of generality, $\kappa=1$. Secondly, we neglect the third equation, since the remaining three close among themselves. Thirdly. we define the rescaled variables

[TABLE]

and consistently

[TABLE]

With this simple rescaling, the equations read

[TABLE]

showing no residual dependence on the driving parameter $\Omega$. A generic solution for $R_{m}$ of Eqs. (S86)-(S89) will be a function of time, $\Omega$ and the initial conditions (denoted in the following with an additional subscript “[math]”), i.e.

[TABLE]

Correspondingly, for $Q_{m}$ and Eqs. (S92)-(S94) we can write

[TABLE]

Out of the rescaling transformation (S90) we can thus connect the two expressions above according to

[TABLE]

Via a completely analogous procedure, one also finds

[TABLE]

Since the left-hand sides are independent of $\Omega$, the right-hand sides are as well. This implies that dynamic trajectories at different values of $\Omega$ can be identified up to an appropriate rescaling of the initial condition and of the overall amplitude, e.g.

[TABLE]

which, re-expressing the initial condition in terms of the original variables, can be recast as

[TABLE]

Before proceeding, it is important to mention a consequence of the physical constraint $G_{m}+R_{m}+N_{m}=1$ $\forall m$: if a process at frequency $\Omega^{\prime}$ starts from a vanishing density of immunes $N_{0k}\equiv 0$, the corresponding rescaled one at frequency $\Omega\neq\Omega^{\prime}$ will instead feature a non-vanishing one

[TABLE]

In order to apply the considerations above to the problem under study, we repeat here the main assumptions we make. First of all, from the numerical dynamics we have noticed that, to a reasonable approximation, every outbreak seems to produce a spatially uniform density of immune sites and, moreover, it appears that outbreaks always propagate from a very localized patch of infected sites in the very neighborhood of the center of the lattice. Hence, we impose that

• (I)

The density of immunes after the $(j-1)$-th outbreak is $N_{m}\equiv N^{(j)}$ for all sites $m$ in an appropriate large neighborhood of the center.

• (II)

Denoting for brevity with an index $\mathbf{k}$ the central site of the lattice, at some time $t_{j}$ after the $(j-1)$-th outbreak, the aforementioned large neighborhood in (I) reaches a configuration in which the center is the only site featuring a non-vanishing infected density $R_{\bf{k}}>0$ and vanishing healthy density and coherence $G_{\bf{k}}=\Sigma_{\bf{k}}=0$, whereas for all the others $R_{m}=\Sigma_{m}=0$ and $G_{m}>0$.

Furthermore, from the stationary data we have verified that the decreasing branch of $N_{\rm{SS}}$ after the first jump is reasonably well fitted by a power-law behavior $\propto 1/\Omega$. Therefore, we further assume that

• (III)

The density $N^{(2)}$ left behind by the first outbreak obeys $N^{(2)}=A_{1}/\Omega$ for some $A_{1}>0$ for all $\Omega$ such that a first outbreak is produced ($\Omega>\widetilde{\Omega}_{1}$ in the main text notation).

We now consider a process at frequency $\Omega$ which starts from the same initial conditions we employ throughout the main text, i.e.,

[TABLE]

and which produces at least $(j-1)$ outbreaks. In the following, for the sake of brevity we shall refer to these initial conditions as the “standard” ones. By (I) and (II), at some time $t_{j}$, in an approximately circular neighborhood of the center which reaches up to the $(j-1)$-th wavefront, the process satisfies

[TABLE]

Since the equations of motion do not explicitly depend on time, we can reset our clock at $t_{j}$ and consider the subsequent dynamics as a new process which starts at $t=0$ from initial conditions given by (S103), for which we thus introduce the notation $R_{k}(t_{j})\to R_{0k}^{(j)}$, $G_{k}(t_{j})\to G_{0k}^{(j)}$, $\Sigma_{k}(t_{j})\to\Sigma_{0k}^{(j)}$ $\forall k$. We note now that

[TABLE]

i.e. the initial conditions of the process at $t_{j}$ are rescaled by $1-N^{(j)}<1$ with respect to the standard ones. Hence, we can now use Eq. (S100) in order to identify a rescaled process which shows the same dynamical behavior up to multiplication by an overall amplitude and starts from standard initial conditions, evolving under a frequency which we rename $\Omega^{\prime}\to\Omega^{(j)}$ (as in the main text) to keep in mind that the process we are interested in is evolving after the production of $(j-1)$ outbreaks. Comparing the arguments of the l.h.s. of Eq. (S100) with (S103), we find

[TABLE]

implying

[TABLE]

as reported in the main text. Recalling the notation $\widetilde{\Omega}_{j}$ for the position of the $j$-th jump, we can draw now some predictions: if $\Omega^{(j)}<\widetilde{\Omega}_{1}$ the equivalent process evolving under $\Omega^{(j)}$ will not produce any outbreak. Consequently, the process we were originally interested in will neither, meaning that the process at fixed $\Omega$ will stop after $(j-1)$ outbreaks and produce a stationary density $N_{\rm{SS}}=N^{(j)}$ up to subextensive contributions generated by the finite propagation of the central infection, which becomes negligible in the thermodynamic limit. We thereby see that, by definition, $\Omega<\widetilde{\Omega}_{j}$. If instead $\widetilde{\Omega}_{1}<\Omega^{(j)}<\widetilde{\Omega}_{2}$, the equivalent process at $\Omega^{(j)}$ will produce a single outbreak, after which it will stop. Hence, the original process will produce, after $t_{j}$, a $j$-th outbreak, but a $(j+1)$-th will not take place. This implies that $\widetilde{\Omega}_{j}<\Omega<\widetilde{\Omega}_{j+1}$. This identification can be further extended to any number of outbreaks. By recalling that $N^{(j)}$ is the density of immunes left behind by the $(j-1)$-th outbreak taking place in a process of frequency $\Omega$, and is thus a function of $\Omega$ itself ($N^{(j)}=N^{(j)}(\Omega)$), we can also write the relation

[TABLE]

between the positions of the jumps. Since the l.h.s. does not depend on $j$, this is easily generalized to

[TABLE]

We specialise this to $i=0$

[TABLE]

and to $i=1$

[TABLE]

Now, we recall that — up to subextensive contributions which vanish in the thermodynamic limit — $N^{(j)}(\Omega)=N_{\rm{SS}}(\Omega)$ $\forall\,\,\widetilde{\Omega}_{j-1}<\Omega<\widetilde{\Omega}_{j}$. Calling $\widetilde{N}_{j}^{(\pm)}=N_{\mathrm{SS}}(\Omega\to\widetilde{\Omega}_{j}^{\pm})$ the top ($+$) and bottom ($-$) of the $j$-th jump, we find from Eqs. (S108) and (S109) that

[TABLE]

These relations, extended to generic $\Omega$s, identify the enveloping black curves in Fig. 4(a) in the main text. More precisely, using the first jump as a reference, and recalling that $N_{1}^{(-)}=0$, we find the expressions used in the plot

[TABLE]

where we numerically estimated $\widetilde{\Omega}_{1}\approx 0.8571\kappa$ (where we have reinstated the units $\kappa$) and $N_{1}^{(+)}\approx 0.898$.

We now make use of assumption (III) as well. First of all, this tells us that $N_{\rm{SS}}(\Omega)=A_{1}/\Omega$ $\forall\,\,\widetilde{\Omega}_{1}<\Omega<\widetilde{\Omega}_{2}$, which allows us to fix the constant $A_{1}$ in such a way that the fitting curve passes through the tip of the first jump: $A_{1}=\widetilde{\Omega}_{1}\widetilde{N}_{1}^{(+)}$. We now proceed by induction: we assume that $N^{(j)}(\Omega)=A_{j-1}/\Omega$ $\forall\,\,\Omega>\widetilde{\Omega}_{j-1}$ — which corresponds to (III) for $j=2$ — and we wish to show the same is true for $N^{(j+1)}$. To this end, we consider a process at frequency $\Omega$ after its $(j-1)$-th outbreak started, such that the density of immunes in the proximity of the center is $N^{(j)}=A_{j-1}/\Omega$. According to Eq. (S106), its equivalent process will take place at a rescaled frequency

[TABLE]

We now rewrite Eq. (S100) as

[TABLE]

Since $R_{m}+G_{m}+N_{m}=1$ at all times, we can also write

[TABLE]

Now, the first outbreak produced by the process in the r.h.s. corresponds to the $j$-th one for the process in the l.h.s., leading to the identity

[TABLE]

for the densities of immunes left behind by either. Since by (III) we know that

[TABLE]

this implies

[TABLE]

By the inductive step in the formulation (S113), we then arrive at

[TABLE]

or equivalently

[TABLE]

which, with the identification $A_{j}=A_{1}+A_{j-1}$, concludes our inductive proof. In particular, this also tells us that $A_{j}=jA_{1}=j\widetilde{\Omega}_{1}\widetilde{N}^{(+)}_{1}$.

Now, exploiting the fact that $N^{(j+1)}\left(\Omega\right)=N_{\rm{SS}}\left(\Omega\right)$ $\forall\,\,\widetilde{\Omega}_{j}<\Omega<\widetilde{\Omega}_{j+1}$, we find

[TABLE]

where $\theta$ denotes the Heaviside step function ($\theta(\omega>0)=1$, $\theta(\omega<0)=0$). The expression above corresponds to the blue line in Fig. 4(a) in the main text.

The considerations above also allow us to estimate the threshold frequencies $\widetilde{\Omega}_{j}$ for $j>1$, given the properties of the first jump, i.e. the values $\widetilde{\Omega}_{1}$ and $\widetilde{N}_{1}^{(+)}$. These values correspond to the points where the (blue) theoretical curve intersects the (black) enveloping ones (S112). Setting $i=1$ in Eq. (S107), we find

[TABLE]

i.e.

[TABLE]