Multicritical Infection Spreading

TL;DR

This paper investigates the multicritical behavior of the disordered contact process, revealing universal properties and activated dynamical scaling, and establishing its connection to the multicritical quantum Ising model.

Contribution

It uncovers the universal multicritical behavior of the disordered contact process and links it to the quantum Ising model universality class.

Findings

Universal multicritical behavior observed in 2D and 3D.

Activated dynamical scaling with exponents matching strong disorder RG predictions.

Connection established between classical multicritical contact process and quantum Ising model.

Abstract

The contact process is a simple infection spreading model showcasing an out-of-equilibrium phase transition between a macroscopically active and an inactive phase. Such absorbing state phase transitions are often sensitive to the presence of quenched disorder. Traditionally, a phase transition in the disordered contact process is either triggered by dilution or by locally varying the infection rate. However, when both factors play an important role, a multicritical point emerges that remains poorly understood. Here, we study the multicritical contact process by large-scale Monte Carlo simulations in two and three dimensions. The multicritical behavior is found to be universal and exhibits ultra-slow, activated dynamical scaling, with exponents consistent with those predicted by the strong disorder renormalization group method. This finding indicates that the multicritical contact process belongs to the same universality class as the multicritical quantum Ising model, opening future directions to measure quantum entanglement properties via classical simulations.

Tables

Table 1. Table 1 : Multicritical exponents. Comparison of the critical exponents from (i) Monte Carlo simulations of the multicritical CP (MCP), (ii) quantum Monte Carlo study of the quantum Ising model (QMC), and (iii) strong disorder renormalization group results for the quantum Ising model (SDRG). Our top four exponents are measured directly through power-law scaling, the next three exponents are directly measured but are impacted by a non-universal microscopic time scale, while the rest are calculated from the values above. The number in the parentheses gives an estimate of the uncertainty on the last digits. (“ − - ” is used for exponents of unknown value.)

	$d=2$					$d=3$
	MCP [32]	QMC [33]	SDRG [34]	MCP (this work)		SDRG [34]	MCP (this work)
${\left(\psi \overline{\delta }\right)}^{-1}$	1.08(14)	1.50(1)	1.258(5)	1.27(5)		0.685(9)	0.65(7)
$\overline{\Theta }/\overline{\delta }$	0.15(28)	0.99(2)	0.516(9)	0.55(7)		0.05(3)	0.03(12)
$\psi \overline{\Theta }$	0.14(24)	0.66(1)	0.410(6)	0.43(4)		0.08(4)	0.03(19)
$\beta$	0.81(7)	0.865(9)	1.099(7)	1.12(4)		1.64(4)	1.44(22)
$\overline{\Theta }$	0.25(45)	$-$	0.579(25)	0.61(10)		0.09(4)	-0.03(17)
$\psi$	0.57(4)	$-$	0.708(20)	0.73(7)		0.93(2)	1.05(20)
$\overline{\delta }$	1.63(10)	$-$	1.12(4)	1.11(7)		1.57(6)	1.6(4)
$d_f$	1.07(12)	1.332(4)	1.205(3)	1.21(3)		1.54(2)	1.49(9)
${\nu }_{\Delta }$	0.88(10)	1.295(5)	1.382(7)	1.42(11)		1.123(10)	0.94(24)
${\overline{\nu }}_{\parallel }$	0.50(9)	$-$	0.98(3)	1.01(9)		1.04(3)	0.9(4)
${\nu }_p$	$-$	$-$	1.168(10)	$-$		0.86(1)	$-$

Equations

$$\rho (t)\sim [\ln (t/t_0){]}^{-\bar{\delta }},$$

$$N(t)$$

$$R(t)$$

$$P(t)$$

$$R$$

$$N$$

$$\frac{1}{\psi \bar{\delta }}$$

$$\frac{\bar{\Theta }}{\bar{\delta }}$$

$$\psi \bar{\Theta }$$

$${\frac{\partial \ln N}{\partial \lambda }|}_{\lambda ={\lambda }_c}\sim P^{-1/\beta }.$$

$$\rho (\Delta )\sim {\Delta }^{\beta },$$

$${\xi }_{\Delta }\sim |\Delta {|}^{-{\nu }_{\Delta }},$$

$$\ln {\xi }_{\parallel }\sim {\xi }_{\Delta }^{\psi },$$

$$\rho (L,t,\Delta )=b^{-x}\tilde{\rho }(L/b,\ln t/b^{\psi },\Delta b^{1/{\nu }_{\Delta }}),$$

$$P(L,t,\Delta )=b^{-x}\tilde{P}(L/b,\ln t/b^{\psi },\Delta b^{1/{\nu }_{\Delta }}).$$

$$N(L,t,\Delta )=b^{d-2x}\tilde{N}(L/b,\ln t/b^{\psi },\Delta b^{1/{\nu }_{\Delta }}).$$

$$R(L,t,\Delta )=b\tilde{R}(L/b,\ln t/b^{\psi },\Delta b^{1/{\nu }_{\Delta }}).$$

$$\bar{\delta }\equiv x/\psi$$

$$\bar{\Theta }\equiv (d-2x)/\psi .$$

$${\xi }_{p-p_c}\sim |p-p_c{|}^{-{\nu }_p}.$$

Full text

Introduction

Even with recent events, our understanding of infection spreading remains limited. As one of the simplest models, the contact process (CP) [1, 2] captures the ability of each agent to infect their neighbors at some infection rate and, once infected, the ability to become healthy, dictated by a recovery rate. On a lattice, the CP is equivalent to the susceptible-infected-susceptible (SIS) model, motivated by scenarios where the agents can get infected multiple times, as is the case for the common cold and COVID. Most theoretical and computational studies of infection spreading focus on the simplified case in which all agents share the same infection and recovery rates — a far cry from real-life scenarios. In the absence of heterogeneities, the CP exhibits an out-of-equilibrium absorbing state phase transition between a macroscopically active and a fully inactive phase that falls in to the universality class of directed percolation (DP) [3, 4, 5, 6].

A perturbative argument, known as the Harris criterion [7], predicts that a phase transition is sensitive to weak time-independent quenched disorder in $d$ dimensions whenever $d\nu<2$, where $\nu$ is the correlation length critical exponent [8]. As a result, in the presence of sufficiently strong disorder, a qualitatively new emergent behavior is expected, characterized by a new exponent $\nu^{\prime}\geq 2/d$ [8]. For the DP universality class, the Harris criterion implies that the critical behavior is modified by weak quenched disorder below four dimensions.

Disorder is typically introduced to the CP either by randomly drawing individual infection (and/or healing) rates from a given distribution [9, 10] or by diluting of the underlying lattice [11, 12]. Depending on the introduced heterogeneities, two distinct non-equilibrium transitions can be observed in the disordered CP. A generic transition takes place in the case of random infection rates or when the dilution of the lattice is below the percolation threshold; a percolation transition is observed when the lattice is critically diluted. As shown in Fig. 1, the two transitions meet at a multicritical point, which is the focus of this work.

At both transitions, the DP critical behavior is replaced by ultra-slow, activated scaling—manifestation of an infinitely disordered fixed point (IDFP) [14, 15]—where the role of time is played instead by its logarithm. While the percolation IDFP is characterized by the critical exponents of the standard percolation [16], the generic IDFP has a new set of “quantum percolation” critical exponents that are known through the so-called strong disorder renormalization group (SDRG) method [14, 15]. The SDRG approach was first developed for disordered quantum systems [17, 18] and was shown by Daniel Fisher to be asymptotically exact in the vicinity of the critical point in one dimension [19, 20]. The disordered CP allows for a similar, but not identical, SDRG treatment [10, 21] that is expected to lead to the same universality class for strong enough disorder below four dimensions [22]. In one dimension, the disordered CP follows the SDRG predictions [23]. In both two and three dimensions, the simulation results of the CP at the generic transition [11, 12, 13] indicate an IDFP with exponents that are close to those obtained by numerical implementations of the SDRG [24, 25, 26, 27, 28]. However, it is an open problem whether the IDFP is attractive for arbitrarily weak disorder [29], or if instead there exists a weak disorder regime, characterized by a line of disorder-dependent fixed points [10, 21, 30, 31].

In this work, we aim to characterize the emergent behavior of the disordered CP at the multicritical point in two and three dimensions (see Fig. 1). While it is not a priori clear whether the SDRG provides accurate predictions for the CP above one dimension, recent implementations of the SDRG method for the multicritical point of the quantum Ising model [34, 35, 36] are a natural point of comparison for our results. In $d=2$, the only simulations of the multicritical CP [32] led to critical exponents that are inconsistent with the SDRG predictions [34]. Moreover, the exponents obtained in Ref. [32], which we report in Table 1, violate the Harris criterion [8]. While systems that violate the Harris criterion are known to exist [37, 38], careful testing would be needed to ascertain whether it is the case for the multicritical CP. To the best of our knowledge, there are currently no simulation results for the $d=3$ multicritical CP, and SDRG predictions [34, 35] are the only available results.

Alternative approaches to the study of the disordered quantum Ising model—such as quantum Monte Carlo or tensor network methods—remain extremely challenging and are restricted to relatively small system sizes. A recent quantum Monte Carlo study of the diluted quantum Ising model [33] considered $d=2$ systems up to linear size $L=64$, obtaining estimates for a subset of the critical exponents at the multicritical point, reported in Table 1. However, as acknowledged by the authors, one must be careful in interpreting these results, which are subject to biases that go beyond the reported statistical error. As evident from Table 1, the three independent sets of estimates are typically not consistent within the reported errors, calling for further studies to characterize the quantitative aspects of this universality class.

The rest of this paper is organized as follows. In Sec. I we define the contact process, summarize the relevant scaling expectations, and describe the Monte Carlo simulations. In Sec. II we present our Monte Carlo results and the extrapolated multicritical exponents in two and three dimensions, followed by a Discussion. Detailed scaling expectations at an IDFP are summarized in the Appendix A.

I Model and Methods

Let us consider a $d$-dimensional hypercubic lattice of sites, each of which is in one of two states: active (infected) or inactive (healthy). Time-evolution is modeled as a continuous-time Markov chain, where every infected node heals at a rate $\mu$ and every healthy node is infected at a rate $\lambda n/2d$ where $n$ is the number of infected nearest-neighbors. Without loss of generality, the healing rate $\mu$ can be set to one, hence the statistical behavior of the system is determined solely by the infection rate $\lambda$. For small infection rates, the system is unable to sustain the presence of active sites indefinitely. The healing rate eventually pushes the dynamics toward a fully inactive state, from which it cannot escape. On the other hand, when $\lambda$ is large, the system is able to maintain a stationary density of active sites [3]. The two regimes are separated by a non-equilibrium phase transition, occurring at a critical infection rate $\lambda_{c}$.

When studying the CP, the initial conditions matter, even informing which quantities can be conveniently studied. One natural option is that of a fully-active initial state. In this case, the fraction of active sites $\rho(t)$ will start at one and decrease over time, eventually settling at some finite value when $\lambda>\lambda_{c}$ or decaying to zero when $\lambda<\lambda_{c}$ [5]. At the critical infection rate, $\rho(t)$ is usually expected to decay over time as a power-law. However, in the presence of strong enough quenched disorder—such as dilution of the lattice—this behavior is replaced by the so-called activated scaling [14] of an IDFP, in which the fraction of infected sites decays as a power-law of the logarithm of time:

[TABLE]

where $t_{0}$ is a microscopic time scale. For a summary of the relevant scaling expectations at an IDFP, see the Appendix A.

An alternative (and often more informative) choice of initial conditions is that of a single infected site. In this “spreading” case, above the critical infection rate, there is a finite probability that the active cluster will spread from the initial active site to span the entire system. When $\lambda<\lambda_{c}$, the system will always eventually reach an inactive state. Note that, due to the stochastic nature of the CP, the system can become trapped in a fully inactive state even when the infection rate is above $\lambda_{c}$. In fact, when the infection rate is only slightly larger than $\lambda_{c}$, most realizations will meet this fate, but a finite fraction, however small, will sustain indefinitely. With these initial conditions, there are three quantities of interest: the average number of active site $N(t)$, the root-mean-square distance of the infected sites from the initial infection site $R(t)$, and the probability $P(t)$ that the system contains at least one active site after time $t$. $N(t)$ and $R(t)$ should both be taken to be averages over an ensemble of realizations. When $\lambda<\lambda_{c}$, $N$, $R$, and $P$ all eventually decay to zero; when $\lambda>\lambda_{c}$, $N$ and $R$ grow to infinity while $P$ asymptotically approaches a finite value. At an IDFP, we once again expect the relevant quantities to scale asymptotically as power-laws of the logarithm of time (see Appendix A for more details):

[TABLE]

A well-known duality can be established between the spreading and fully active initial conditions [39, 40], from which it follows that they share the same critical infection rate $\lambda_{c}$ and that $\rho(t)$ and $P(t)$ share the same critical exponent $\bar{\delta}$. Hence, by expressing these observables as a function of each other, we can cancel the dependence on the non-universal microscopic time scale $t_{0}$, obtaining the following power-law scaling forms:

[TABLE]

Note that these exponent combinations are not independent, as they are all related to $\beta/\nu_{\Delta}$ via the relations (Appendix A):

[TABLE]

In this work, we perform Monte Carlo simulations of the two and three-dimensional CP on a critically diluted lattice according to the following, well-established procedure [12, 13]. First, a bond-diluted square (cubic) lattice is generated, in which bonds are independently removed with probability $p_{c}=\frac{1}{2}$ ($p_{c}=0.751188$), i.e. the bond percolation threshold on a square (cubic) lattice [41]. Two-dimensional simulations are run on critically-diluted finite lattices of $2000\times 2000$ sites. However, the edge of the lattice has not been reached by the infected cluster in any of our simulations, so the finite system size plays no role. Three-dimensional simulations are run on a lattice whose size is increased whenever the infected cluster reaches its edge, so that its size is effectively infinite. At time $t=0$, the system is initialized with a single active site, then time-evolution is simulated through repeated random updates. During each update, one active site, $i$, is randomly selected, and heals (becomes inactive) with probability $p_{\mathrm{heal}}=\frac{1}{1+\lambda}$. If site $i$ does not heal, one of its neighbors, $j$, is randomly selected. If site $j$ is inactive and the bond connecting $i$ and $j$ has not been removed, site $j$ is infected (becomes active). Each update corresponds to a time increment $\Delta t=\frac{1}{n_{a}}$, where $n_{a}$ is the current number of active sites in the sample.

II Results

II.1 Critical scaling in two dimensions

We run $10^{6}$ independent simulations for a set of infection rates $3.4\leq\lambda\leq 3.57$ up to a final simulation time $T=10^{5}$, each using an independently generated critically diluted lattice. Close to the estimated location of the multicritical point, we increase the number of simulations to $10^{7}$ and the final time to $T=10^{6}$.

The number of active sites $N$, the radius of the active cluster $R$, and the fraction of surviving realizations $P$ are measured from the simulations. Analyzing these observables directly as functions of time is difficult. The value of the microscopic time scale $t_{0}$—which is a priori unknown—can impact any conclusions drawn from finite-time simulations. Even locating the critical infection rate can pose a challenge: while at criticality the observables are expected to exhibit power-law behavior asymptotically, in practice one finds that finite-time effects make identifying such behavior highly non-trivial. To address this problem, $N$, $R$, and $P$ are often plotted against each other [12, 13]. Since the same $t_{0}$ is expected to apply to all three observables, the value of $t_{0}$ becomes irrelevant, allowing for an easier identification of $\lambda_{c}$. The resulting plots are shown in Fig. 2–2. Each line in these plots traces out the average trajectory of an ensemble of simulations through the space of observables $N$, $R$, and $P$. Note that $P(t)$ is a monotonically decreasing function of time (the fraction of surviving simulations), so along each line, time flows in the negative $P$ (positive $N$, $R$) direction. As expected, we observe trajectories diverging from each other over time: simulations below the multicritical infection rate $\lambda_{c}$ die out ($N,R,P\rightarrow 0$), while a finite fraction of simulations above $\lambda_{c}$ approach a stationary density of active sites ($N,R\rightarrow\infty$, $P\rightarrow\mathrm{const.}$).

To identify the multicritical infection rate, we study the slope of each line as a function of $P$. Near $\lambda_{c}$, we should find lines whose slope remains finite as $P\rightarrow 0$ (i.e., $t\rightarrow\infty$ for $\lambda\leq\lambda_{c}$). In this limit, the slope of each line corresponds to a combination of critical exponents, as given in Eq. (3), whose value should be related to the other two exponent combinations as in Eq. (4). The resulting slopes are shown in Fig. 2–2. The $R$ versus $P$ curves prove to be especially informative. As $P\rightarrow 0$, the corresponding slopes appear to curve upward or downward for larger and smaller infection rates, respectively. Curves with $3.5575\leq\lambda\leq 3.56$ appear approximately linear as functions of $P$, suggesting that this is the location of the critical region. Furthermore, extrapolating as $P\rightarrow 0$, we find a value of the corresponding exponent combination $(\psi\bar{\delta})^{-1}$ in agreement with the SDRG prediction within the uncertainty. The slopes of the $N$ versus $P$ and $N$ versus $R$ curves in this proposed critical region do not appear linear when plotted against $P$. Nevertheless, the $P\rightarrow 0$ estimates of the respective exponent combinations are consistent with the SDRG predictions.

As our best proxy, we estimate the location of the multicritical infection rate of the two-dimensional CP to be at $\lambda_{c}^{(2)}=3.5588(13)$. As an additional test, we check whether the extrapolated values of the exponent combinations are related to one another as in Eq. (4). We find that the only infection rates that yield exponent combinations in agreement with Eq. (4) within the uncertainty are $\lambda=3.5575$ and $3.55875$, further supporting our estimate of the location of the critical region. The numerical results for the three exponent combinations can be found in Table 1, showing good agreement with the SDRG predictions.

We turn now to the task of estimating the microscopic time scale $t_{0}$ and measuring the critical exponents $\bar{\Theta}$, $\psi$, and $\bar{\delta}$ individually. The primary feature of the critical point is the asymptotic activated scaling summarized in Eq. (2). In practice, when plotting $N$, $R$, and $P$ against $\ln(t/t_{0})$, the slope of the critical curve should remain finite in the $t\rightarrow\infty$ limit, with the corresponding critical exponent as its limiting value. Having previously located the critical region, we can in turn use this information to estimate $t_{0}$ by requiring that the slopes of the critical lines allow for good extrapolations to a finite value in the infinite-time limit. Since in the disordered CP all observables are most naturally understood as functions of $\ln(t/t_{0})$, we choose to implement the $t\rightarrow\infty$ limit as $[\ln(t/t_{0})]^{-1}\rightarrow 0$. Through this procedure, we arrive at an estimated range for the microscopic time scale of $5.5\lesssim\ln(t_{0})\lesssim 7$, similar to previous estimates of the microscopic time scale in the bond-diluted CP [12].

Fig. 3 shows plots of $N(t)$, $R(t)$, and $P(t)$ and their corresponding slopes using $\ln(t_{0})=6.3$. Building on this estimate of $t_{0}$, we can perform infinite-time extrapolations of the slopes of the near-critical lines to arrive at numerical estimates of the critical exponents. Two sources of uncertainty need to be taken into account when performing these extrapolations: on the value of the critical infection rate and on the value of the microscopic time scale. The resulting $\bar{\Theta}$, $\psi$, and $\bar{\delta}$ are reported in Table 1; they are all compatible with the SDRG predictions within their respective uncertainties.

II.2 Critical scaling in three dimensions

To study the critical behavior of the three-dimensional multicritical CP, we run at least $10^{4}$ independent simulations for a set of infection rates $5\leq\lambda\leq 6.2$ up to a final simulation time of at least $T=10^{4}$, each using an independently generated critically diluted lattice. As we approach the estimated location of the multicritical point, we progressively increase the number of simulations and the final time, reaching $10^{6}$ simulation up to $T=10^{6}$ for near-critical lines. Compared to our two-dimensional simulations, the resulting curves are noisier, owing to both smaller sample sizes and fast-decreasing fractions of surviving simulations. The analysis of the three-dimensional simulations is analogous to the two-dimensional case, as detailed in Section II.1.

The multicritical point is located by plotting the observables $N$, $R$, and $P$ against each other and studying the slope of each curve as a function of $P$ (Fig. 4). As in $d=2$, the $R$ versus $P$ plot proves to be the most informative. Curves with $5.58\leq\lambda\leq 5.6$ appear approximately linear as functions of $P$, and linear extrapolations as $P\rightarrow 0$ are consistent with the SDRG estimate of the corresponding exponent combination, $(\psi\bar{\delta})^{-1}$. Likewise, the slopes of the two remaining sets of curves in the critical region, though not flat, extrapolate to values consistent with the SDRG predictions in the $P\rightarrow 0$ limit, as reported in Table 1. As our best estimate, we identify the location of the multicritical point of the three-dimensional CP to be at $\lambda_{c}^{(3)}=5.59(1)$. Checking for consistency, we find that the only set of exponent combinations in agreement with Eq. (4) corresponds to $\lambda=5.6$.

Next, using the critical lines as a guide, we estimate the microscopic time scale $t_{0}$. We require that the slopes of the near-critical lines appear approximately linear when plotted against $[\ln(t/t_{0})]^{-1}$ to allow for a good infinite-time ($[\ln(t/t_{0})]^{-1}\rightarrow 0$) extrapolation, since they are expected to remain finite in this limit. Through this procedure, we arrive at a range of estimates for the microscopic time scale of $6\lesssim\ln(t_{0})\lesssim 8$. We note that, due to the noisier curves, estimating $t_{0}$ for the three-dimensional model proved to be more challenging than in the two-dimensional case, yielding a more approximate range. Taking into account the uncertainty on $t_{0}$ as well as on $\lambda_{c}$, we extrapolate the slopes of $N(t)$, $R(t)$, and $P(t)$ to find numerical estimates of each of their corresponding critical exponents (Fig. 5). The resulting values, reported in Table 1, are all in agreement with the SDRG estimates within the error bars.

II.3 Off-critical scaling

To access more critical exponents, we estimate $\beta$ by studying the scaling of the observables away from the critical infection rate. With a fully-active initial state, the stationary density of active sites is the relevant order parameter, scaling with the control parameter $\lambda$ as $\rho(\lambda)\sim|\lambda-\lambda_{c}|^{\beta}$. Duality implies [39, 40] that in the case of spreading initial conditions, the survival probability $P(\lambda)$ follows the same scaling. However, going after the dependency of $P$ on the infection rate directly requires measuring the infinite-time, stationary survival probability close to the critical point, which is made challenging by the ultra-slow dynamical scaling. To get around this problem, we follow Ref. [13] and instead study the derivative of the logarithm of the observables with respect to $\lambda$. The relevant scaling form for this analysis is:

[TABLE]

Analogous formulas hold for $\partial_{\lambda}\ln P$ and $\partial_{\lambda}\ln R$. These formulas are especially useful since they do not rely on our estimate of $t_{0}$ (with its uncertainty).

Our analysis of the off-critical scaling behavior proceeds as follows. First, we consider a set of near-critical $N(\lambda,t)$ curves and use their distance from a candidate critical curve to estimate $\partial_{\lambda}\ln N$ as a function of time. To increase the accuracy of our estimate, we average the result over several near-critical curves, but exclude the two curves closest to the critical one, since the corresponding distances are greatly affected by statistical noise. $\partial_{\lambda}\ln N$ is then plotted against $P$ and the slope of the resulting curve is used to estimate $\beta$ by extrapolating to $P\rightarrow 0$. The full procedure is summarized in Fig. 6.

In the $d=2$ model, we perform this analysis separately for all three curves within the critical region ($3.5575\leq\lambda\leq 3.56$). We find that the slope of $\partial\ln N/\partial\lambda$ as a function of $P$ is initially close to the CP Monte Carlo results of Ref. [32] and to the quantum Monte Carlo results of Ref. [33]. However, as $P\rightarrow 0$, we observe a crossover to a different scaling behavior, occurring approximately when the $N(\lambda,t)$ curves cross the microscopic time scale $t_{0}$ (Fig. 6). Performing a linear extrapolation of the slope of $\partial_{\lambda}\ln N$ as $P\rightarrow 0$ yields a value of $\beta$ consistent with the SDRG prediction [34]. The resulting numerical estimate of $\beta$ is reported in Table 1, where the primary source of the reported statistical error is the uncertainty on the value of $\lambda_{c}^{(2)}$.

The same analysis is carried out for each curve in the proposed critical region of the $d=3$ model ($5.58\leq\lambda\leq 5.6$). The larger uncertainty on $\lambda_{c}$, as well as the noisier $N(\lambda,t)$ curves, lead to a larger uncertainty on the extrapolated value of $\beta$. The resulting numerical estimate, reported in Table 1, is nonetheless consistent with the value predicted by the SDRG [34]. Carrying out the same analysis using $\partial_{\lambda}\ln P$ or $\partial_{\lambda}\ln R$ in place of $\partial_{\lambda}\ln N$ yields similar, if noisier, results.

Discussion

In this work, we have determined the critical exponents of the multicritical CP through large-scale Monte Carlo simulations in both two and three dimensions. In two dimensions, there have been three independent estimates of the critical behavior, by the strong disorder renormalization group (SDRG) [34], a quantum Monte Carlo simulation [33] and a CP simulation [32]. In three dimensions—to the best of our knowledge—there have been no prior estimates for the multicritical exponents.

In both $d=2$ and $d=3$, our exponents are within the error consistent with those obtained by the SDRG method for the multicritical quantum Ising model [34]. This serves as the first confirmation of the multicritical SDRG results, indicating that the approximations applied in the numerical SDRG scheme are indeed valid. This opens the possibility that our methodology could also show an improved agreement between the critical CP simulations and the SDRG results.

In two dimensions, our exponents satisfy the Harris criterion, a qualitative difference from those obtained by Dahmen et al. [32]. Note that in Ref. [32], the Monte Carlo simulations only used the largest connected component of each critically diluted lattice. However, as discussed in Ref. [34], this methodological difference only leads to a small shift in a subset of the exponents. The main factor explaining the difference between our results and those of Ref. [32] is the choice of a fully-infected initial condition. In this type of simulations, the density of active sites is the only accessible observable and any interpretation of the results, including locating the critical region, is greatly impacted by the microscopic time scale.

As for the quantum Monte Carlo results of Ref. [33], future work would ideally improve on three aspects: (i) obtain larger sizes for the estimates to converge, (ii) characterize the missing exponents — most importantly quantifying $\psi$ to check the presence of activated scaling — and (iii) improve the reported error bars by incorporating systematic effects in addition to the currently included statistical error.

The observation that the disordered (critical and multicritical) CP and the quantum Ising model likely belong to the same universality class has widespread consequences. Prominently, it opens a direct way to measure in equilibrium quantum entanglement effects in the classical out-of-equilibrium CP, as the dynamics is governed by the same underlying activity clusters [42]. Specifically, subsystem entanglement has been characterized recently in both $d=2$ and $d=3$ for the critical [27] and multicritical [35] quantum Ising model, and was found to show a universal dependence on the shape of the subsystem [36]. However, there is no direct quantum measurement that could validate these results. Therefore, a simple classical model that could test universal quantum information results is a very attractive possibility. Characterizing quantum entanglement between multiple subsystems is another rich direction that has been recently demonstrated to provide universal results for the disordered quantum Ising chain [43].

Acknowledgements.

We thank W. T. Engedal and R. Juhász for helpful comments and discussion. This work was supported by the National Science Foundation under Grant No. PHY-2310706 of the QIS program in the Division of Physics, and Northwestern’s Office of Undergraduate Research with the generous support of the Hung-Farinelli Family under the 2023 Summer Undergraduate Research Grant.

Appendix A Scaling expectations at an IDFP

Here we provide a summary of the relevant scaling relations at an IDFP, following Refs. [10, 21]. The control parameter is traditionally chosen as $\Delta_{0}\equiv\overline{\ln(\lambda/\mu)}$ [14], with the overbar standing for a disorder average. Note that here, we use $\mu=1$. The deviation from the critical point $\Delta_{0}^{*}$ is then given by $\Delta\equiv\Delta_{0}-\Delta_{0}^{*}$. The $\rho$ activity density of the system can be chosen as the order parameter, vanishing as

[TABLE]

close to the transition in the active phase. As usual, the spatial correlation length $\xi_{\Delta}$ diverges close to the critical point as

[TABLE]

with the correlation length exponent $\nu_{\Delta}$. At an IDFP, the time and length scales are related to each other as

[TABLE]

with $\psi$ being the tunneling exponent, leading to $\overline{\nu}_{\parallel}=\nu_{\Delta}\psi$. The order parameter for finite size $L$ and time $t$ is expected to transform as

[TABLE]

when the length is rescaled by a factor $b$. Hence, we obtain the relation $x\equiv\beta/\nu_{\Delta}$. Then, the fractal dimension is given by $d_{f}=d-x$.

The survival probability is expected to scale as the order parameter [5, 6]

[TABLE]

With the binary variables $n_{i}$ coding active ($n_{i}=1$) and inactive ($n_{i}=0$) sites, the spatio-temporal correlation function $C[n_{0}(t=0),n_{i}(t)]$ summed over the position $i$ yields for the scaling of the average number of active nodes:

[TABLE]

The spread is defined as the root-mean-square of the distance of active nodes from the origin, averaged for the surviving samples up to time $t$. Using the same scaling form as the spatio-temporal correlation function we obtain

[TABLE]

At the critical point ($\Delta=0$) of the infinite system ($L=\infty$), it follows from the above relations that the observables depend on time asymptotically as given by Eq. (2), with the exponents given in terms of the earlier ones as

[TABLE]

and the hyperscaling relation

[TABLE]

At the multicritical IDFP, there are two relevant directions, corresponding to two control parameters. In addition to moving away from the multicritical point by changing the infection rates (controlled by $\Delta$), we can also move away by changing the dilution parameter, $p$, leading to an additional correlation length exponent $\nu_{p}$.

[TABLE]