Critical behavior of a two-step contagion model with multiple seeds

TL;DR

This paper analyzes how multiple infection seeds influence the critical behavior of contagion models, revealing a transition from hybrid to continuous epidemic spread as seed density increases, with distinct critical exponents.

Contribution

It provides analytical and numerical insights into the phase transition nature of multi-seed contagion models, extending understanding beyond single-seed scenarios.

Findings

Hybrid epidemic transition at low seed density

Transition becomes continuous at high seed density

Distinct critical exponents from single-seed models

Abstract

A two-step contagion model with a single seed serves as a cornerstone for understanding the critical behaviors and underlying mechanism of discontinuous percolation transitions induced by cascade dynamics. When the contagion spreads from a single seed, a cluster of infected and recovered nodes grows without any cluster merging process. However, when the contagion starts from multiple seeds of $O(N)$ where $N$ is the system size, a node weakened by a seed can be infected more easily when it is in contact with another node infected by a different pathogen seed. This contagion process can be viewed as a cluster merging process in a percolation model. Here, we show analytically and numerically that when the density of infectious seeds is relatively small but $O(1)$, the epidemic transition is hybrid, exhibiting both continuous and discontinuous behavior, whereas when it is sufficiently large and reaches a critical point, the transition becomes continuous. We determine the full set of critical exponents describing the hybrid and the continuous transitions. Their critical behaviors differ from those in the single-seed case.

Figures

Tables

Table 1. Table 1: List of the critical exponents of the SWIR models with a single seed and with multiple seeds.

	type	$\beta$	${\beta }^{\prime }$	$\gamma$	${\gamma }^{\prime }$	$\overline{\nu }$	${\overline{\nu }}^{\prime }$
${\rho }_0=1/N$	single seed	0	-	-	-	3	-
	multiple seeds	1/2	-	$1$	-	$2$	-
${\rho }_0={\rho }_c$	multiple seeds	1/3	1/3	$1.47\pm 0.05$	$1.28\pm 0.05$	$2.13\pm 0.1$	$2.13\pm 0.1$

Equations

$$\mathrm{S}+\mathrm{I}$$

$$\mathrm{W}+\mathrm{I}$$

$$\mathrm{I}$$

$$m(\kappa )={\rho }_0+(1-{\rho }_0)\sum _{k=1}^{\infty }{P}_d(k)\sum _{\ell =1}^k\left(\genfrac{}{}{0ex}{}{k}{\ell }\right)q^{\ell }(1-q{)}^{k-\ell }{P}_R(\ell ).$$

$$P_S(\ell )=(1-\kappa -\mu {)}^{\ell }.$$

$$P_W(\ell )=\sum _{j=0}^{\ell -1}(1-\kappa -\mu {)}^j\mu (1-\eta {)}^{\ell -j-1},$$

$$q_{n+1}={\rho }_0+(1-{\rho }_0)\sum _{k=1}^{\infty }\frac{k{P}_d(k)}{z}\sum _{\ell =1}^{k-1}\left(\genfrac{}{}{0ex}{}{k-1}{\ell }\right)q_n^{\ell }(1-q_n{)}^{k-1-\ell }{P}_R(\ell )\equiv {\rho }_0+(1-{\rho }_0)f(q_n),$$

$$f(q_{n})=1-\Big{(}1+\dfrac{\mu}{\eta-\kappa-\mu}\Big{)}e^{-(\kappa+\mu)q_{n}z}+\dfrac{\mu}{\eta-\kappa-\mu}e^{-\eta q_{n}z}.$$

$$F(m,{\rho }_0)\equiv f(m)-\frac{m}{1-{\rho }_0}+\frac{{\rho }_0}{1-{\rho }_0}.$$

$$F(m,{\rho }_0)=\frac{{\rho }_0}{1-{\rho }_0}+am+b{m}^2+c{m}^3+O(m^4),$$

$$a$$

$$b$$

$$c$$

$$F(m,{\rho }_0)\equiv \frac{{\rho }_0}{1-{\rho }_0}+am+b{m}^2+c{m}^3.$$

$$m^{\pm }=-\frac{b}{3c}\pm \sqrt{\frac{b^2}{9c^2}-\frac{a}{3c}},$$

$$m_0^{-}(\infty )-⟨m_0^{-}(N)⟩\sim N^{-\beta /\bar{\nu }}$$

$$F(\kappa_{c},m_{0})=\dfrac{dF}{dm}\bigg{|}_{\kappa_{c},m_{0}}=\dfrac{d^{2}F}{dm^{2}}\bigg{|}_{\kappa_{c},m_{0}}=0.$$

$$F(\kappa_{c}+\delta\kappa,m_{0}+\delta m)\simeq\dfrac{\partial F}{\partial\kappa}\bigg{|}_{\kappa_{c},m_{0}}(\delta\kappa)+\dfrac{1}{6}\dfrac{\partial^{3}F}{\partial m^{3}}\bigg{|}_{\kappa_{c},m_{0}}(\delta m)^{3}+\dots=0$$

Full text

I Introduction

Nonequilibrium dynamic transitions driven by cascade dynamics on complex networks have attracted considerable attention recently review1 ; review2 ; review3 . The spreading of epidemic disease on complex networks review_epidemics ; watts ; dodds ; krapivsky ; grassberger_pre ; althouse_prl ; chen ; porter ; grassberger_nphy ; althouse_pnas ; janssen ; janssen_spinodal ; chung ; hasegawa1 ; choi_2016 is an instance, in which a pathogen is transmitted from an infected node (e.g., a person) to a susceptible neighbor, who then becomes infected with a certain probability. If the transmission probability is sufficiently large (small), the pathogen spreads out to a macroscopic scale (remains local). An epidemic transition occurs between these two limits. The extent of spreading also depends on the structure of an underlying network review1 ; barabasi . When degree distribution of a network is highly heterogeneous, diseases can spread out massively even for a small transmission probability, so that an epidemic transition point can be zero sf . Information spreading in social media from one page to others may be modeled in a similar manner watts ; dodds .

Among the several epidemic models, one of the simple contagion models is the so-called susceptible-infected-recovered (SIR) model sir ; sir_newman , in which each node has one of three states, susceptible (denoted as $S$), infected ($I$), or recovered ($R$). Initially, all the nodes are in state $S$ except for one seed node in state $I$. The contagion process starts from a single node in state $I$. Each node in state $I$ transmits pathogens to its neighbors in state $S$ and infects each of them with probability $\kappa$; then, it changes its state to $R$ with unit probability. This contact process is repeated until the system reaches an absorbing state in which no infected node is left in the system. When the probability $\kappa$ is sufficiently small (large), the order parameter defined as the density of nodes in state $R$ after the system falls into the absorbing state, becomes $o(N)$ [$O(N)$]; i.e., the system falls into a subcritical (supercritical) state. In between, an epidemic transition occurs at $\kappa_{c}$, and the system exhibits critical behavior. It is known that when the dynamics starts from a single seed on Erdős-Rényi (ER) random networks ER , the SIR model undergoes a continuous percolation transition following the universal behavior of ordinary percolation.

The SIR model with multiple seeds has been considered hasegawa2 , in which two percolation transitions occur successively at $\kappa_{c1}$ and $\kappa_{c2}$ as $\kappa$ is increased. The density of nodes in state $R$ is finite for $\kappa>\kappa_{c1}$, whereas the density of nodes in state $S$ disappears for $\kappa>\kappa_{c2}$. Thus, there exists a state of coexisting nodes in states $R$ and $S$ between $\kappa_{c1}$ and $\kappa_{c2}$.

The SIR model was extended to a two-step contagion model, in which a weakened state ($W$) can exist between the $S$ and $I$ states. Accordingly, this model is called the SWIR model janssen ; janssen_spinodal . Nodes in state $W$ are involved in the reactions $S+I\to W+I$ and $W+I\to 2I$, which occurs in addition to the reactions $S+I\to 2I$ and $I\to R$ in the SIR model. The properties of the epidemic transition in the SWIR model were extensively investigated for the single-seed case grassberger_pre ; janssen ; janssen_spinodal ; chung ; hasegawa1 ; choi_2016 . The order parameter defined as the density of nodes in state $R$ after an absorbing state is reached, displays a discontinuous transition, whereas other physical quantities such as the outbreak size distribution exhibit critical behaviors. Thus, the phase transition occurring in the SWIR model with a single seed is regarded as a mixed-order phase transition choi_2016 . The dynamic rule of the SWIR model is rather so simple that its underlying mechanism for the discontinuous behavior of the order parameter was disclosed universal . Moreover, the mechanism turned out to be universal in other models such as $k$-core percolation kcore1 ; kcore2 ; kcore3 ; baxter_prx , the cascading failure model on interdependent networks buldyrev ; son ; zhou ; bk ; review_mcc , and the epidemic-related models watts ; dodds ; krapivsky ; grassberger_pre ; althouse_prl ; chen ; porter ; grassberger_nphy ; althouse_pnas .

Here, we investigate the phase transitions of the SWIR model with multiple seeds. The model with multiple seeds has been investigated in Refs. hasegawa1 ; janssen_spinodal ; hasegawa3 : The authors of Refs. hasegawa1 ; hasegawa3 used the mean-field approach and performed numerical simulations, obtaining the phase diagram as a function of the reaction rates. The order parameter exhibits either a discontinuous or continuous transition depending on the density of the infectious seeds and mean degree of a given network hasegawa1 ; hasegawa3 . In Ref. janssen_spinodal , the discontinuous transition is regarded as a spinodal transition, because there is no co-existence phase in the system while the order parameter jumps. Even though such results were obtained, the properties of the phase transitions and critical behaviors were not deeply investigated yet.

Here, we reveal that the spread of contagion in the SWIR model with multiple seeds proceeds differently from that in the SWIR model with a single seed: in the multiple-seed case, the reactions $W+I\to 2I$ often occur even in early time steps, because nodes in states $W$ and $I$ involved in that reaction can originate from different seeds (see Fig. 1). We note that the number of multiple seeds was taken as $O(N)$. On the contrary, in the single-seed case, such reactions rarely occur until the system reaches a characteristic dynamic step $n_{c}(N)\sim N^{1/3}$: When dynamic step $n$ is less than $n_{c}(N)$, the reactions $S+I\to 2I$ and $I\to R$ are dominant but the number of nodes in $R$ still remains as $o(N)$. The contagion spreads in the form of a branching tree. When the dynamics reaches $n_{c}(N)$, the branching tree forms long-range loops due to finite-size effect. Once such loops form, the reaction $W+I\to 2I$ occurs massively, in which the nodes in state $W$ were generated in early time steps. Thus, the density of nodes in state $R$ increases drastically as many as $O(N)$ in short time steps. Due to these different contagion mechanisms, the properties of epidemic transitions in the multiple seed case become different from those in the single seed case. We will determine the full set of critical exponents describing the phase transitions in the multiple seed case, and compare them with those obtained in the single seed case choi_2016 .

This paper is organized as follows: In Sec. II, we present the rules of the SWIR model in detail. In Sec. III, we set up the self-consistency equation to derive the mean-field solution using the local tree approximation of the order parameter for the epidemic transition on the ER networks. We show that, depending on the initial density of infectious nodes, different types of phase transition can occur. In Sec. IV, we report numerical results for the epidemic transitions. In the final section, a summary and discussion are presented.

II The SWIR model

The SWIR model with multiple seeds is simulated on ER networks with $N$ nodes. Initially, $N\rho_{0}$ nodes are selected randomly from among those $N$ nodes and assigned to state $I$; the other nodes are assigned to state $S$. At each time step $n$, the following processes are performed. (i) All the nodes in state $I$ are listed in random order. (ii) The states of the neighbors of each node in the list are updated sequentially as follows: If a neighbor is in state $S$, it changes its state in one of the two ways: either to $I$ with probability $\kappa$ or to $W$ with probability $\mu$. If a neighbor is in the state $W$, it changes to $I$ with probability $\eta$, where $\kappa$, $\mu$, and $\eta$ are the contagion probabilities for the respective reactions. (iii) All nodes in the list change their states to $R$. This completes the time step, and we repeat the above processes until the system reaches an absorbing state in which no infectious node is left in the system. The reactions are summarized as follows:

[TABLE]

The order parameter exhibits a discontinuous transition at a transition point $\kappa_{c}$ when $\rho_{0}$ is less than a critical value $\rho_{c}$, and it shows a continuous transition at $\kappa_{c}$ when $\rho_{0}=\rho_{c}$ for given parameter values $z$, $\mu$, and $\eta$, where $z$ is the mean degree of a given ER network.

III Self-consistency equation and physical solutions

In an absorbing state, each node is in one of three states: the susceptible $S$, weakened $W$, and recovered $R$ states. The order parameter $m(\kappa)$, the density of nodes in state $R$ in an absorbing state, is written using the local tree approximation as

[TABLE]

The first term in Eq. (5), $\rho_{0}$, is the initial density of infected nodes. In the second term, the factor $(1-\rho_{0})$ represents the probability that a node is originally in state $S$. $P_{d}(k)$ is the probability that a randomly selected node has degree $k$; $q$ is the probability that an arbitrarily chosen edge leads to a node that is in state $R$ but not infected through the chosen edge in the absorbing state. Thus, $P_{d}(k)\binom{k}{\ell}q^{\ell}(1-q)^{k-\ell}$ is the probability that a node has degree $k$ and $\ell$ of them are in state $R$ in the absorbing state. $P_{R}(\ell)$ is the conditional probability that a node is finally in state $R$, provided that it was originally in state $S$ and its $\ell$ neighbors are in state $R$ in the absorbing state.

Similarly to $P_{R}(\ell)$, we define $P_{S}(\ell)$ as the conditional probability that a node remains in state $S$ in the absorbing state, provided that it has $\ell$ neighbors in state $R$ and was originally in state $S$. $P_{W}(\ell)$ is defined similarly. We note that for a certain node to have $\ell$ neighbors in state $R$ in the absorbing state means that the node receives $\ell$ attempts to infect it when the recovered neighbors are in state $I$. Thus, a node still remaining in state $S$ with $\ell$ neighbors in state $R$ has to be unchanged from $\ell$ infection attempts through the entire process. Thus, we obtain

[TABLE]

Next, the probability $P_{W}(\ell)$ is given as

[TABLE]

where $j$ denotes the number of attacks that a node sustains before it changes to state $W$. Using the relation $P_{S}(\ell)+P_{W}(\ell)+P_{R}(\ell)=1$, one can determine $P_{R}(\ell)$ in terms of $P_{S}$ and $P_{W}$.

The local tree approximation allows us to define $q_{n}$ similarly to $q$ but now at the tree level $n$. The probability $q_{n+1}$ can be derived from $q_{n}$ as follows:

[TABLE]

where the factor $kP_{d}(k)/z$ is the probability that a node connected to a randomly chosen edge has degree $k$. As a particular case, when the network is an ER network with $P_{d}(k)=z^{k}e^{z}/k!$, $f(q_{n})$ becomes

[TABLE]

Eq. (8) reduces to a self-consistency equation for $q$ for given epidemic parameter values in the limit $n\to\infty$. Once we obtain the solution of $q$, we can obtain the outbreak size $m(\kappa)$ using Eq. (5). For ER networks, however, $m(\kappa)$ becomes equivalent to $q$, thus the solution of the self-consistency equation Eq. (8) yields the order parameter. We remark that the methodology we used here is similar to those used in previous studies of epidemic spreading on complex networks dodds ; chung ; janssen_spinodal ; grassberger_pre ; hasegawa1 .

Hereafter, we set $\mu=\kappa$ for convenience and define a function

[TABLE]

Using formula (9), we approximate $F(m,\rho_{0})$ in the limit $m\to 0$ as

[TABLE]

where

[TABLE]

For convenience, we neglect the higher-order terms and redefine $F(m,\rho_{0})$ as

[TABLE]

Depending on the relative magnitudes of $a$ and $b$, various solutions of the self-consistency equation $F(m,\rho_{0})=0$ can be obtained. However, we need to check whether these solutions are indeed physically relevant in the steady state when we start epidemic dynamics from a certain initial condition. The stability criterion was established in a previous work choi_2016 : The solution $F(m_{0},\rho_{0})=0$ is stable if and only if $\partial F(m,\rho_{0})/\partial m|_{m=m_{0}}<0$.

The equation of state in the steady state can be obtained using $F(m,\rho_{0})=0$. We notice that $F(m=0,\rho_{0})=\rho_{0}/(1-\rho_{0})>0$ and $F(m=\infty,\rho_{0})=-\infty$ because $c<0$, as shown in Fig. 2. We examine the solutions of $dF(m,\rho_{0})/dm=0$, which are obtained as

[TABLE]

where $a$, $b$, and $c$ are given in formulas (12)–(14). Note that $a$ depends on $\rho_{0}$. At these extreme points $m^{\pm}$, $F(m,\rho_{0})$ has either a local maximum or a local minimum. For a given $\rho_{0}$, $z$, and $\eta$, both $m^{\pm}$ values exist, and they are positive in the range of $\kappa_{d}<\kappa<\kappa_{a}$, where $b^{2}/9c^{2}-a/3c=0$ at $\kappa=\kappa_{d}$, and $a=0$ at $\kappa=\kappa_{a}$. For a given $z$ and $\eta$, diverse types of phase transitions occur depending on $\rho_{0}$. When $\rho_{0}$ is less than a certain value $\rho_{c}$, the order parameter jumps at a transition point. On the other hand, when $\rho_{0}\geq\rho_{c}$, the order parameter increases continuously with $\kappa$. At $\rho_{0}=\rho_{c}$, $m^{+}=m^{-}=m_{0}$ and $F(m_{0},\rho_{c})=0$ at $\kappa=\kappa_{d}=\kappa_{c}$, as schematically shown in the blue (lower) curve in Fig. 3.

III.1 When $\rho_{0}<\rho_{c}$

When $\rho_{0}<\rho_{c},$ there exists a range of $\kappa$ in which $F(m,\rho_{0})=0$ has more than one solution, as shown in Fig. 2. The order parameter $m$ versus $\kappa$ is shown in Fig. 4(a) and (b). In particular, when $\kappa$ has a certain value $\kappa_{c}^{-}$, $m^{-}$ obtained using Eq. (16) satisfies $F(m^{-},\rho_{0})=0$. The $m^{-}$ value at $\kappa_{c}^{-}$ is denoted as $m_{0}^{-}$. We also define $\kappa_{c}^{+}$ and $m_{0}^{+}$ similarly to $m^{+}$ in Eq. (16). We note that $\kappa_{c}^{+}<\kappa_{c}^{-}$. Depending on the magnitude of the reaction probability $\kappa$ relative to $\kappa_{c}^{+}$ and $\kappa_{c}^{-}$, the order parameter behaves differently, as follows:

i) For $\kappa<\kappa_{c}^{+}$, there exists one stable solution $m=m_{d}(\kappa)$, which increases slowly with $\kappa$. It is obtained that ${m_{d}}\approx\rho_{0}/(1-\kappa z)+O(\rho_{0}^{2})$.

ii) At $\kappa=\kappa_{c}^{+}$, there exist two solutions, $m_{d}$ and $m_{0}^{+}$ ($m_{d}<m_{0}^{+}$). However, $m_{0}^{+}$ is not accessible because $m_{d}$ is stable.

iii) When $\kappa_{c}^{+}<\kappa<\kappa_{c}^{-}$, there exist three solutions, $m_{d}$, $m_{m}$, and $m_{u}$, with relative magnitudes $m_{d}<m_{m}<m_{u}$; however, the solution $m_{m}$ is unstable. Thus, only $m_{d}$ is accessible from the initial density $\rho_{0}<m_{d}$. The order parameter behaves as $m_{0}^{-}-m_{d}(\kappa)\sim(\kappa_{c}^{-}-\kappa)^{1/2}$ for $\kappa<\kappa_{c}^{-}$. Thus, the critical exponent of the order parameter is obtained as $\beta=1/2$.

iv) At $\kappa=\kappa_{c}^{-}$, there exist two stable solutions, $m_{0}^{-}$ and $m_{u}$. Thus, the order parameter jumps between the two values, exhibiting discontinuous behavior. Hence, a hybrid phase transition occurs at the point $(\kappa_{c}^{-},m_{0}^{-})$.

v) For $\kappa>\kappa_{c}^{-}$, there exists one solution, denoted as $m_{u}(\kappa)$, which increases with $\kappa$ as $m_{u}(\kappa)-m_{u}(\kappa_{c}^{-})\sim(\kappa-\kappa_{c}^{-})$.

III.2 When $\rho_{0}=\rho_{c}$

When $\rho_{0}=\rho_{c}$, there exists a reaction probability $\kappa_{c}$ that satisfies the relation $F(m_{0},\rho_{c})=dF(m,\rho_{c})/dm|_{m=m_{0}}=0$, and $b^{2}-3ac=0$. Thus, the two solutions, $m_{0}^{-}$ and $m_{0}^{+}$, reduce to the same one, which is denoted as $m_{0}$. The function $F(m,\rho_{c})$ versus $m$ is shown in Fig. 3, and the order parameter $m$ versus $\kappa$ is shown with the analytic solution and simulation data in Fig. 5. At $\kappa_{c}$, singular behavior occurs, and the order parameter $m$ behaves as $m-m_{0}\sim|\kappa-\kappa_{c}|^{1/3}$ on both sides. The derivation of this exponent $1/3$ is presented in the Appendix.

IV Numerical results

To estimate various critical exponents, we perform extensive numerical simulations on ER networks with mean degree $z=8$. For simplicity, the reaction probability $\mu$ is set equal to $\kappa$, and $\eta=1/2$. With these parameter values, we determine $\rho_{c}$ as precisely as $0.00747762$, which we will use in numerical analysis later. For $\rho_{0}<\rho_{c}$, we take $\rho_{0}=2\times 10^{-3}$ in the simulations. We take the average over 50 different dynamics samples for each of 1,600–4,000 network configurations. Thus, 80,000–200,000 configuration averages were taken to obtain each data point.

IV.1 When $\rho_{0}<\rho_{c}$

We found analytically that the order parameter behaves as $m_{0}^{-}-m_{d}(\kappa)\sim(\Delta\kappa)^{\beta}$ with $\beta=1/2$ in the thermodynamic limit, where $\Delta\kappa\equiv\kappa_{c}^{-}-\kappa$. The main panel of Fig. 6 shows $m_{0}^{-}-\langle m_{d}(\kappa)\rangle$ versus $\Delta\kappa\equiv\kappa_{c}^{-}-\kappa$ in a double logarithmic scale. Data points in the figure are obtained numerically from systems of several selected sizes $N$, and the dashed line is obtained from the analytic solution of Eq. (8) by taking the limit $n\to\infty$, which is valid in the thermodynamic limit. We find that the data points saturate to constant values asymptotically as $\Delta\kappa\to 0$, whereas they overlap with the dashed curve as $\Delta\kappa$ is increased. As shown in the inset, also in a double logarithmic scale, the dashed curve follows the line with a slope of $0.5$ in the region $\Delta\kappa<\Delta\kappa^{*}\approx 10^{-4.2}$; however, it deviates from the line in the opposite region beyond $\Delta\kappa^{*}$. This fact implies that conventional finite-size scaling analysis is valid for systems with size larger than $O(10^{8})$. However, it would be impractical to perform simulations with such huge system sizes.

Following the conventional finite-size scaling theory,

[TABLE]

at $\kappa_{c}^{-}$. We check this relation in Fig. 7. For small system sizes $N$, $\beta/\bar{\nu}$ seems to be about 0.2, whereas it is estimated to be $\approx 0.24$ for large $N$. Again the crossover occurs between the system sizes $N=10^{7}$ and $10^{8}$. We could obtain a more reliable value for the exponent ratio $\beta/\bar{\nu}$ for somewhat larger system sizes, but that is impractical.

The fluctuation of the order parameter $\chi(\kappa)\equiv N(\langle m_{d}^{2}\rangle-\langle m_{d}\rangle^{2})$ diverges as $\sim(\kappa_{c}^{-}-\kappa)^{-\gamma}$. For finite systems of size $N$, it is expected that $\chi\sim N^{\gamma/\bar{\nu}}$ at $\kappa=\kappa_{c}^{-}$. From the simulation data, we obtain $\gamma/\bar{\nu}\approx 0.5$, as shown in Fig. 8.

With the measured values $\beta/\bar{\nu}\approx 0.24$ and $\gamma/\bar{\nu}\approx 0.5$ and the analytic result $\beta=1/2$, we guess $\bar{\nu}=2$ and then $\gamma=1$. If we use those values, then the hyperscaling relation $2\beta+\gamma=\bar{\nu}$ would hold.

IV.2 When $\rho_{0}=\rho_{c}$

At $\rho_{0}=\rho_{c}$, the jump in the order parameter does not appear, and $m_{0}^{+}=m_{0}^{-}$ at $\kappa=\kappa_{c}$ in the thermodynamic limit. In finite systems, however, the order parameter can still exhibit a jump in some samples. Thus, the order parameter distribution $p(m)$ accumulated over different samples exhibits two separate peaks, as shown in Fig. 9. We regard the data points of $p(m)$ in the region $m<m_{0}$ ($m>m_{0}$), where $m_{0}$ has the theoretical value $0.171405\dots$, as those obtained from $m_{0}^{-}(N)$ ($m_{0}^{+}(N)$) for different samples. At $\kappa=\kappa_{c}$, in finite systems, we obtain the power-law behaviors $m_{0}-\langle m_{0}^{-}(N)\rangle\sim N^{-\beta/\bar{\nu}}$ with $\beta/\bar{\nu}\approx 0.153$ (Fig. 10) and $\langle m_{0}^{+}(N)\rangle-m_{0}\sim N^{-\beta^{\prime}/\bar{\nu}^{\prime}}$ with $\beta^{\prime}/\bar{\nu}^{\prime}\approx 0.164$ (Fig. 11).

For $\kappa<\kappa_{c}$, the fluctuation of the order parameter $\chi\equiv N(\langle m_{d}^{2}\rangle-\langle m_{d}\rangle^{2})$ behaves as $\sim N^{\gamma/\bar{\nu}}G[(\kappa_{c}-\kappa)N^{1/\bar{\nu}}]$ with a certain scaling function $G$. On the other hand, for $\kappa>\kappa_{c}$, we obtain that $\chi^{\prime}\equiv N(\langle m_{u}^{2}\rangle-\langle m_{u}\rangle^{2})$ behaves as $\sim N^{\gamma^{\prime}/\bar{\nu}^{\prime}}G^{\prime}[(\kappa-\kappa_{c})N^{1/\bar{\nu}^{\prime}}]$ with a certain scaling function $G^{\prime}$. We numerically obtain $\gamma/\bar{\nu}\approx 0.69$ (Fig. 12) and $\gamma^{\prime}/\bar{\nu}^{\prime}\approx 0.6$ (Fig. 13).

On the basis of the numerically obtained values $\beta/\bar{\nu}\approx 0.53$ and $\gamma/\bar{\nu}\approx 0.69$, and the theoretical value $\beta=1/3$, we estimate $\bar{\nu}\approx 2.179$ and $\gamma\approx 1.5$. Those values are confirmed in Fig. 14 for $\chi N^{-\gamma/\bar{\nu}}$ versus $(\kappa_{c}-\kappa)N^{1/\bar{\nu}}$, in which the data collapse well with the choices of $\bar{\nu}\approx 2.13\pm 0.1$ and $\gamma/\bar{\nu}\approx 0.69$. The measured values of the exponents satisfy the hyperscaling relation $(2\beta+\gamma)/\bar{\nu}\approx 0.996$ well. Similarly, for $\kappa>\kappa_{c}$, on the basis of the numerical values $\beta^{\prime}/\bar{\nu}^{\prime}\approx 0.164$ and $\gamma^{\prime}/\bar{\nu}^{\prime}\approx 0.6$, and the theoretical value $\beta^{\prime}=1/3$, we obtain $\bar{\nu}^{\prime}\approx 2.03$ and $\gamma^{\prime}\approx 1.218$. Data for $\chi^{\prime}$ for different system sizes collapse well into a single curve with the choices of $\bar{\nu}^{\prime}=2.13\pm 0.1$ and $\gamma^{\prime}/\bar{\nu}^{\prime}=0.6$ (Fig. 15). These values yield $(2\beta^{\prime}+\gamma^{\prime})/\bar{\nu}^{\prime}\approx 0.91-0.93$, which deviates slightly from the expected value of unity that would satisfy the hyperscaling relation. To obtain those results, we used the numerical values $\rho_{c}\approx 0.00747762$ and $\kappa_{c}\approx 0.108021$. We remark that $\beta=\beta^{\prime}=1/3$ is obtained analytically.

V Summary and discussion

We investigated the properties of phase transitions in the SWIR model with a finite density $\rho_{0}$ of initially infected seeds janssen . A node in the state $S$ can change its state to weakened ($W$) or infected ($I$) when it comes in contact with an infected node from the same or a different root. A weakened node can also change its state to infected ($I$) when it contacts an infected node from the same or a different root. The reaction probabilities $\kappa$ and $\mu$ in Eqs. (1) and (2), respectively, serve as control parameters. For convenience, we take $\kappa=\mu$. We found that for a given network, there exists a critical density of seeds $\rho_{c}$ such that for $\rho_{0}<\rho_{c}$, the order parameter, the density of nodes in state $R$ in the absorbing state, increases continuously with the critical exponent $\beta=1/2$ as $\kappa$ is increased up to a transition point $\kappa_{c}^{-}$ and then jumps to a finite value, followed by a continuous increase. Accordingly, the order parameter behaves as $m(\kappa)=m_{0}^{-}-b(\kappa_{c}^{-}-\kappa)^{1/2}$ for $\kappa<\kappa_{c}^{-}$, where $b$ is a positive constant. At $\kappa_{c}^{-}$, the order parameter is discontinuous by $\Delta m=m_{u}(\kappa_{c}^{-})-m_{0}^{-}$. Thus, the order parameter itself exhibits a hybrid phase transition. This pattern is different from that for the single-seed case, in which the order parameter jumps from $m=0$ to a finite value, and thus $\beta=0$. The fluctuation of the order parameter diverges at the transition point $\kappa_{c}^{-}$ according to a power-law with the exponent $\gamma$. For the correlation size exponent $\bar{\nu}$ measured in finite systems, we find that the hyperscaling relation $2\beta+\gamma=\bar{\nu}$ holds reasonably well.

As $\rho_{0}$ is increased, the jump shrinks and becomes zero at $\rho_{c}$. For $\rho_{0}=\rho_{c}$, the transition becomes continuous. We determined a complete set of critical exponents describing the phase transition at $\kappa_{c}$. The critical exponents are listed in Table I.

Acknowledgements.

This work was supported by the National Research Foundation of Korea by Grant No. NRF-2014R1A3A2069005.

Appendix A Derivation of the critical exponent $\beta$ at $\rho_{c}$

Here we introduce an analytical method to determine the critical exponents $\beta=1/3$ at $\rho_{0}=\rho_{c}$. It is already noted in Sec. III-B that for $\rho_{0}=\rho_{c}$,

[TABLE]

We consider a line of the solution $F(\kappa,m)=0$ near $(\kappa_{c},m_{0})$ by expanding $F(\kappa_{c}+\delta\kappa,m_{0}+\delta m)$ as

[TABLE]

where only nonzero terms are considered. Since $\delta\kappa$ and $(\delta m)^{3}$ are two lowest terms in Eq. (19) and their coefficients have the opposite sign to each other, $\delta m\sim(\delta\kappa)^{1/3}$ when $\delta\kappa\ll 1$. Thus for both cases of $\kappa<(>)\kappa_{c}$, the critical exponents $\beta=\beta^{\prime}=1/3$.