Epidemics greatly affect Homo sapiens [1] [2] [3] , as well as all other animals. Epidemic modeling goes back at least to Bernoulli who modeled the spread of smallpox [4, 5] . The methods and concepts developed in this field have been applied to modeling recent human epidemics like HIV, H1N1 Swine Flu, Zika Virus, and COVID-19, and also to modeling non-biological processes like the spread of rumors or cultural fads, diffusion of innovations, computer viruses, etc. (see [6] [7] [8] [9] and references therein).
One general lesson that we learn about epidemics is that they often operate close to a critical regime. Subcritical epidemics are actually the most numerous, but they quickly die out. Supercritical epidemics can kill a finite fraction of the population thereby destroying the environment that the virus needs to thrive and multiply. Epidemics close to the critical regime are kind of optimal in a view of never-ending competition between the hosts and the viruses. Furthermore, humans are now sufficiently advanced to devise and employ the containment measures to suppress supercritical epidemics to critical.
Mathematical models of epidemics are over-simplified. In physics, for instance, we know that electrons are identical. In contrast, organisms are different, even twins are different. This heterogeneity is rarely taken into account, and it is far from clear how to model it in a reasonable way. In physics, we usually rely on binary and symmetric interactions. Both these features are questionable in the realm of epidemics. Other realistic features are also mostly ignored. However, the populations where the epidemics spread are usually very large and the lore from statistical physics tells us that in large systems qualitative behaviors can be predicted even if one greatly simplifies the model. This is especially true in critical regimes.
The very concept of critical regimes comes from epidemic modeling. This concept clearly emerges from the well-known susceptible-infected-recovered (SIR) process [10] [11] [12] [13] [14] , a toy model that mimics the spread of infection. According to the rules of the SIR process, infected individuals recover (become immune or die) with equal rates and every infected individual transmits a disease to every susceptible individual with the rate R 0 /N , where N is the population size. Thus each infected individual on average spreads the infection to R 0 individuals before recovery. Therefore the behavior of the SIR process greatly depends on whether the reproduction number R 0 is smaller or larger than the recovery rate which we set to unity. When R 0 < 1, i.e., for subcritical SIR processes, outbreaks quickly end, namely just a few individuals catch the disease. For supercritical SIR processes (R 0 > 1), the outbreak may affect only a few individuals, e.g. starting from a single infected individual the size of the outbreak is finite with probability R −1 0 . With complimentary probability, rN + O( √ N ) individuals catch the disease before the outbreaks dies out; the fraction r = r(R 0 ) is implicitly determined by r + e −R0r = 1 (1) Huge outbreaks killing finite fractions of the population continue to devastate animal species. They also used to decimate human societies [1] [2] [3] , e.g., the Black Death killed about 50% of the European population [3] . Preventive and containment measures such as quarantine, improved hygiene, etc. suppress supercritical infectious diseases and often drive them to a critical situation. This critical state is effectively self-organized. Indeed, suppressing the disease to the subcritical regime may be possible but costly and psychologically difficult to maintain when the number of newly infected starts to decrease exponentially. Therefore, if the outbreak is not quickly and completely eradicated, the containment measures are relaxed and the system may return to the supercritical stage, the disease gets again out of control, so the containment measures are tightened driving the system back to the subcritical state. It would be interesting to devise a self-organized process of the spread of infection with dynamics similar to the critical SIR process. In this paper, however, we merely consider the critical SIR process with many initially infected individuals.
The SIR processes are often treated using a deterministic framework [10] [11] [12] [13] [14] . This framework can be applied to the supercritical regime where it gives e.g. the simplest derivation of Eq. (1). Stochastic effects are unavoidable, however, for the critical and subcritical SIR processes (see [15] [16] [17] [18] ). When the population of susceptible is finite, finite-size corrections become important, particularly for the critical SIR process [19] [20] [21] [22] [23] [24] [25] .
In this paper, we study the critical SIR process and hence we employ stochastic methods. We consider finite populations. The population size is assumed to be large, N 1. We focus on the situation when the initial number of infected individuals is also large: k 1. Most epidemics start with a single infected individual. We want to describe the SIR process that has begun at the supercritical regime and exhibited an exponential growth regime. Preventive and containment measures subsequently suppressed the reproduction number to unity. Ignoring the earlier regime yields the critical SIR process with a certain large number k of initially infected individuals.
The same critical SIR process with a large number of initially infected individuals has been recently studied, using large-scale simulations and scaling arguments, by Radicchi and Bianconi [26] . Our analysis relies on exact calculations and asymptotic methods. Our analytical and asymptotic predictions qualitatively agree with simulation results [26] , explain the major features, and perhaps suggest slightly different scaling fits which are a little simpler than the fits used in Ref. [26] . The chief reason for subtle behaviors are algebraic tails, the average size and duration of outbreaks are especially sensitive to these tails.
Considerable insight into the behavior of finite systems can be gained from the analysis of the infinite-population limit where the critical SIR process is equivalent to the critical branching process with continuous time. This is a classical subject, but we are studying an unusual setting with an arbitrary number k of initially infected individuals. We derive several exact results in this infinitepopulation limit. Our analysis of finite systems employs asymptotic and scaling methods.
The outline of this paper is as follows. In Sec. II, we consider the infinite-population limit, present exact results for the outbreak size distribution and show that exact results approach a simple scaling form in the most interesting situation when k 1. We then analyze the critical SIR process in a population with N 1 individuals. The size of an outbreak is a random quantity, its average and variance are studied using scaling and heuristic arguments (Sec. III) and asymptotically exact analysis (Sec. IV). In Sec. V, we investigate the duration of outbreaks. Several technical calculations are relegated to the Appendices A-C.
In the infinite-population limit, the SIR process reduces to the branching process [27] [28] [29] [30] [31] . Branching processes involve duplication and death. Symbolically,
For the critical branching process, the rates of duplication and death are equal.
Branching processes have numerous applications, e.g., they mimic cell division and death; for adult organisms, the critical branching process is appropriate as the number of cells remains (approximately) constant.
We begin with a classical setting when one individual was initially infected. Let A n be the probability that exactly n individuals catch the infection before the epidemic is over. With probability 1 2 , the initially infected individual joins the population of recovered before infecting anyone else, so A 1 = 1 2 . Further, A 2 = 1 2 A 2 1 since at the first step a new individual must get infected, and then both must recover without spreading infection. Proceeding along these lines we arrive at the recurrence
reflecting that the first infection event creates two independent infection processes [28] . A solution to (2) is found by introducing the generating function
converting the recurrence (2) into a quadratic equation
whose solution reads
Expanding A(z) in powers of z we find
In particular, the probabilities A n are given by for n = 1, . . . , 11. Generally when the critical branching process begins with k initially infected individuals, infection processes originated with each individual are independent. Hence the probability A (k) n that exactly n individuals catch the infection before the epidemic is over can be expressed via the probabilities A m ≡ A (1) m corresponding to the classical situation with one initially infected individual:
The generating function
is therefore
This generating function encapsulates all A (k)
n , but it is desirable to obtain more explicit representation.
Needless to say, A (k) n = 0 when n < k. When n ≥ k, one can express A are a first few explicit formulas:
The general formula is
where x denotes the largest integer ≤ x. Plugging (6) into the sum in Eq. (10), one reduces the sum to a hypergeometric series [32] . Using identities involving hypergeometric series [32], as well as identities involving the gamma function (particularly, the duplication formula), we have computed the sum in (10) and arrived at a neat final formula
Note that substituting (6) into the recurrence (7) one can directly compute These results are recovered from the general solution (11) thereby providing a consistency check. As another consistency check we note that (11) agrees with normalization
The sequence A (k) n has a single peak located at n = k when k ≤ 4, while for k ≥ 5 the peak is at n = ν(k) > k, see Fig. 1 
at n = ν(k), then decays and eventually approaches
This asymptotic behavior is straightforwardly deduced from (11) . Using the general solution (11) together with the Stirling formula one can establish the behaviors of b(k) and ν(k) in the k → ∞ limit. One gets
with
The general solution (11) approaches the scaling form
in the scaling limit
The scaled distribution Φ(µ) has a single peak (Fig. 2 ) and it vanishes faster than any power of µ when µ → 0.
We also note that from Eqs. (18)- (19) one easily recovers (16)- (17) . Due to the algebraic tail (15) , the moments defined by n a = n≥k n a A (k) n diverge when a ≥ 1 2 . To perform the summation in the a < 1 2 range where the moments converge it is convenient to consider modified moments with n a replaced by Γ(n+1)/Γ(n+1−a). Such moments admit an analytical expression n≥k Γ(n + 1)
Identities which we have used in computing the sum in (20) are standard, they can be found e.g. in Ref.
[32]. When k 1, the modified moments are asymptotically equal to the regular moments since Γ(n+1) Γ(n+1−a) → n a for n ≥ k 1, so Eq. (20) simplifies to
The moments n a remain finite for finite populations, but diverge as N → ∞ if a ≥ 1 2 . This divergent leading behavior of the moments in finite populations is computed in the next section; in the simplest case of k = 1, it is given by Eq. (24) .
Consider the critical SIR process with k initially infected individuals, but in a finite population. Denote by N the size of the population and by A n (k; N ) the probability that the size of the outbreak is n. In the infinitepopulation limit, A n (k; ∞) ≡ A (k) n is given by (11) ; it approaches the scaling form (18)-(19) when k 1. The probability distribution A n (N ) ≡ A n (1; N ) for the critical SIR process starting with a single initially infected individual has been previously investigated in Refs. [21] [22] [23] [24] [25] . The probability disrtibution A n (N ) acquires a scaling form
in the scaling limit
This scaling was proposed and numerically supported in [21] . Kessler and Shnerb [22] derived the scaling function F (ν). The moments n a which diverge in the infinitepopulation limit when a ≥ 1 2 are now finite. In this range, the scaling (22)-(23) implies the following leading asymptotic behavior
Two most important cumulants are the average E k (N ) = n and the variance V k (N ) = n 2 − n 2 . When k = 1, these cumulants exhibit the following growth with the population size
The amplitudes can be (numerically) computed using the exact expression [22] for the scaling function to give
In the general case, the distribution A n (k, N ) depends on three variables. The interesting range is
The first scaling in (27) follows from (23), the second is an outcome of (19) and (23). In the scaling region (27), the distribution A n (k, N ) is expected to acquire a scaling form
with ν = n/N 2/3 , see Eq. (23) , and the scaled initial number of infected individuals κ = k/N 1/3 . More precisely, (28) should hold in the scaling limit (23) and
The normalization condition, n≥k A n (k, N ) = 1, gives ∞ 0 dν G(κ, ν) = 1 and explains the pre-factor in (28) . The two-variable distribution G(κ, ν) is unknown, so let us discuss the average outbreak size which is the basic quantity with simpler scaling behavior. The average outbreak size E k (N ) depends on two variables and its conjectural scaling behavior is
The two scaling behaviors, (28) and (30), are compatible if Ψ(κ) = ∞ 0 dν νG(κ, ν). We now show that the scaled distribution Ψ(κ) has simple extremal behaviors
To establish the small κ asymptotic we note that in the k N 1/3 regime, the infectious processes generated by each initially infected individual are mutually independent. Thus
Comparing (30) and (32) we obtain Ψ(κ) C 1 κ as κ → 0 as asserted in (31).
To establish the large κ behavior we first mention that E N (N ) = N . This suggests that Ψ(N 2/3 ) ∼ N 1/3 , from which Ψ(κ) ∼ √ κ as κ → ∞. This argument is heuristic. In Sec. IV, we derive the large κ asymptotic asserted in (31). In Appendix A, we present an elementary derivation based on the observation that the behavior in the κ → ∞ limit is essentially deterministic. The deterministic analysis is significantly simpler than the exact approach, but it cannot be used to study fluctuations while using the exact approach we also compute the variance (Sec. IV).
Simulation results [26] are in a reasonably good agreement with the (conjectural) scaling behavior (29)- (30) . The numerical data [26] well agree with (31) in the small κ limit: Ψ(κ) = C 1 κ with C 1 ≈ 1.5 in simulations, the analytical prediction is C 1 = 1.4528 . . .. In the large κ limit, numerical data in [26] were fitted to √ κ up to a logarithmic correction.
We now turn to the variance. For sufficiently small k, we rely again on the mutual independence of k infectious processes to deduce
The scaling region is given by (29) . The natural scaling behavior of the variance compatible with (33) is
where Ψ 2 (κ) C 2 κ when κ → 0. The asymptotically exact behaviors in the complimentary κ → ∞ limit is established in Sec. IV. Summarizing,
Simulation results [26] support an algebraic decay in the κ → ∞ limit: Ψ 2 ∼ κ −γ . The numerical uncertainty [26] in the magnitude of the exponent γ is rather significant, γ = 0.75 ± 0.15.
The critical SIR process admits an exact treatment. Denote by s, i and r the population sizes in the susceptible, infected and recovered individuals. The entire population is constant:
Due to the constraint (36), the state of the process can be described by any pair of variables s, i, r. We choose (i, x) with x = N − s. For the critical SIR process, in the interesting regime s is close to N , viz. N − s N , and hence x is a more convenient variable than s. The constraint (36) shows that x = i + r, so x ≥ i.
Infection and recovery events are symbolically
Denote by t(i, x) the number of transitions from the state (i, x) to termination. We are mostly interested in t(k, k), i.e., starting with k infected and no recovered. The process terminates at some state (0, n), where n is the size of the outbreak. The rules (37a) and (37b) show that the quantity i − 2x decreases by 1 in each transition. Thus starting at (i, x) = (k, k) gives i−2x = −k−T after T transitions, and in particular
The rates of the processes (37a) and (37b) imply that they occur with probabilities
The stochastic transition time t(i, x) evolves according to the rules
A. Average number of transitions Averaging (40) we find that T 1 (i, x) = t(i, x) satisfies
To avoid cluttering of formulae, we write p ± ≡ p ± (x) and T 1 ≡ T 1 (i, x) when there is no confusion. The recurrence (41) should be solved subject to the boundary condition
The boundary-value problem (41)-(42) admits an exact solution
One can verify by direct substitution that Eq. (43) with amplitudes determined from the recurrent relations (44) satisfies (41)-(42). In Appendix B, we show that the guess form (43) is rather natural. Specializing (43) to i = x = k we obtain
from which the average size of an outbreak is
Thus, we have found an exact formula (46) for the average size of the outbreak. We have not succeeded in extracting an asymptotic behavior of the sum in (46). There are two technical challenges. First, the amplitudes, which are in principle known from the recurrence relations (44), are unwieldy. Second, the sum involves N − k terms. One can find compact exact results when N − k = O(1), see Appendix B. In the interesting range k = O(N 1/3 ), and hence the number of terms in the sum is huge.
Another line of attack on discrete problems relies on continuum methods. We assume that T 1 (i, x) is a smooth function of i and x, and we expand T 1 (i + 1, x + 1) and T 1 (i − 1, x) appearing in (41) up to the second order
Here ∂ i = ∂/∂i, ∂ x = ∂/∂x, etc. denote partial derivatives. Inserting (47) into (41) and keeping only dominant terms we obtain
Suppose that the scaling in the interesting range is
Plugging (49) into (48) we find that the terms in are comparable only when α = 1 3 and β = γ = 2 3 . Thus we re-scale the variables
and the average number of transitions
One can verify that terms not included in (48) are subdominant. For instance, computing the second derivates gives ∂ i ∂ x T 1 = O(N −1/3 ) and ∂ 2 x T 1 = O(N −2/3 ), so these derivatives can indeed be dropped.
The transformation (50)-(51) turns (48) into a partial differential equation (PDE)
for the re-scaled transition time T(I, X). We must solve (52) in the quadrant I ≥ 0 and X ≥ 0. The boundary condition (42) yields T(0, X) = 0 (53)
Solving the boundary-value problem (52)-(53) is an intriguing challenge that we leave for the future. Here we limit ourselves by a simpler problem of computing the asymptotic behavior of T 1 (k, k) when k N 1/3 . In the realm of the framework (50)-(53), we should learn how to extract the large I behavior of T(I, X). This can be done by noting that when I 1, the diffusion term can be dropped from Eq. (52). Thus we arrive at the first order PDE
Introducing new variables
The solution is
with an arbitrary function f (u). The boundary condition (53) gives T = 0 when v = −u. This fixes f (u) = 2 √ 2u.
We want to determine T 1 (k, k) = N 2/3 T(κ, N −1/3 κ). Thus I = κ 1 and X = 0 as we always consider large populations, N 1. More precisely, setting X = 0 amounts for a tacit assumption κ N 1/3 . Summarizing, our asymptotic results are valid in the range
The upper and lower bounds are well separated when N 1/3 1. This is satisfied for large populations, yet the convergence may be slow as the effective small parameter is N −1/3 . Thus T 1 (k, k) = N 2/3 T(κ, 0) when the bounds (58) are obeyed. Using Eq. (57) we get T 1 (k, k) = √ 8kN which we insert into E k (N ) = [k + T 1 (k, k)]/2 obtained after averaging Eq. (38). Keeping only the leading term gives
2kN . This completes the derivation of the large κ behavior announced in Eq. (31).
Taking the square of Eq. (40) and averaging we find that T 2 (i, x) = t 2 (i, x) satisfies
Again we shortly write p ± ≡ p ± (x) and T 2 ≡ T 2 (i, x). We now subtract the square of Eq. (41) from Eq. (59) and find that the variance
Similar to (47) we expand the variance
Plugging the expansions (47) and (62) into (61) and keeping only dominant terms we obtain
We use the same rescaled variables (50) as before, and seek the variance in the scaling form
The transformation (50) and (64) turns (63) into
When I 1, we can drop again the diffusion term from (65). We also use the asymptotic expression (57) for T(I, X) and arrive at the first order PDE
Using the variables (55) we recast (66) into
We set again I = κ and X = 0 in (68) and arrive at the asymptotic V(κ, 0) = 4 2/κ when κ 1. Using Eq. (38) we find V k (N ) = N 4/3 1 4 V(κ, 0) which is indeed the large κ behavior announced in Eq. (35).
Some basic features of the duration of the outbreaks in the critical SIR process in a finite system can be extracted from the temporal behaviors in the infinite-population limit. We first recall these infinite-population results in the simplest case with one initially infected individual [16, 25] . The probability P i (t) to have i infected individuals at time t satisfieṡ
where dot denotes the derivative with respect to time.
A solution of an infinite set of equations (69a)-(69b) subject to the initial condition P i (0) = δ i,1 reads
This soultion can be verified by a direct substitution, or derived using e.g. generating function techniques [see Appendix C for details]. The probability that the outbreak is still alive at time t, is
The average number of infected individuals in outbreaks which are still alive at time t is therefore
A general solution of Eqs. (69a)-(69b) describing an infinite-population limit subject to an arbitrary number of initially infected individuals is somewhat cumbersome, it is presented in Appendix C.
In a finite population, the infection process eventually comes to an end. To estimate heuristically this final time t f one uses (72) to express [21] the final size of the outbreak through the final time: n f ∼ t f dt i ∼ t 2 f . The maximal outbreak size scales as n * ∼ N 2/3 (see e.g. [21, 22, 25] ) and hence the maximal duration is
The average duration of the outbreak is formally
Recalling that P 1 (t) = (1+t) −2 in the infinite-population limit (equivalently, for the critical branching process), we notice that the integral in (74) diverges. We should use, however, the finite upper limit given by (73). This leads to an estimate
where the subscript reminds that the process begins with a single infected individual. The logarithmic growth of the average duration time was predicted by Ridler-Rowe many years ago [19] , albeit with incorrect amplitude; the correct amplitude 1/3 is easy to appreciate [33], it was argued and numerically supported in [21, 25] . The above argument also suggests the more precise asymptotic
Since the logarithm is a slowly growing function, the subleading constant term c 1 significantly contributes to the average duration. The computation of the sub-leading term requires much more comprehensive analysis than what we have used so far. If the number of initially infected individuals is sufficiently small, k N 1/3 , one can generalize the prediction (75) for the average duration of an outbreak without using a complete solution of the infinite-population limit (Appendix C), it suffices to rely on the independence of infection processes generated by each initially infected individual. The probability that the infection is over at time t is P k 0 , with P 0 given by (70b). Thus −dP k 0 /dt = kP k−1 0 /(1 + t) 2 is the probability density that the infection is eradicated at time t, from which
implying that the average duration of the outbreak exhibits a simple logarithmic scaling with amplitude proportional to the initial number of infected individuals.
To guess the scaling form we notice that (77) can be rewritten as E k [t] k ln(N 1/3 k −1 ) with the same leading accuracy. Thus one plausible scaling form is
with Θ(κ) −κ ln κ when κ → 0. The variance can be probed similarly to the average, see (75). One establishes the scaling law
but not an amplitude. Independence gives V 1 [t] ∼ kN 1/3 when k N 1/3 , and hence the hypothetical scaling form of the variance is
Simulations [26] support the scaling law (80) and the linear small κ behavior. Simulations also suggest [26] that the scaled distribution Θ 2 (κ) is inversely proportional to κ in the large κ limit. Thus
The scaling behavior of the average outbreak duration seems rather tricky. A scaling form with pre-factor being a product of an algebraic and a logarithmic in N terms has been used in [26] to fit simulation data. This scaling form is notably different from three other scaling forms that all have a standard purely algebraic in N pre-factor [cf. Eqs. (30) , (34), and (80)]. It may be worthwhile to try to fit numerical data for the average outbreak duration with the standard scaling form (78), but allowing the possibility of an anomalous small κ behavior such as Θ(κ) −κ ln κ argued above.
We have studied the critical SIR process starting with a large number of initially infected individuals, k 1. Particularly interesting behaviors emerge when k scales as a cubic root of the population size, k ∼ N 1/3 . The critical SIR process exhibits large fluctuations, so we have relied on the stochastic formulation. We have treated the problem using a combination of exact calculations and asymptotic methods. We have focused on the size and duration of the outbreaks, more precisely on the average and variance of these quantities. The analysis of the size of outbreaks is rather detailed, Secs. III-IV.
Our analytical and asymptotic predictions qualitatively agree with simulation results [26] . Whenever there is a little discrepancy between simulation results and theoretical predictions, it seems that data may be fitted using slightly different scaling forms, simpler than the fits used in Ref. [26] . The chief reason for subtle behaviors are algebraic tails, the average size and duration of outbreaks are especially sensitive to these tails.
There are many remaining challenges. For instance, we have found the average size of outbreaks (Sec. IV). The exact expression for the average size involves the sum of N − k terms, with amplitudes determined recurrently but becoming increasingly cumbersome. We have not succeeded in extracting a clean asymptotic behavior of this sum. We have also shown that continuum methods in principle allow one to determine the scaling functions, one should just solve linear PDEs with non-constant coefficients. We have used continuum methods in extracting asymptotic behaviors of the scaling functions.
The derivation of the exact average size in Sec. IV can be probably generalized to establish the variance, and perhaps all cumulants (that is, to compute the cumulant generating function).
[31] P. Haccou, P. Jagers, and V. A. Vatutin, Branching processes: variation, growth, and extinction of populations For the SIR process starting with a single infected individual, the deterministic framework is applicable only in the super-critical regime if the process enters the runaway regime. In this situation, the epidemic ends only after a finite fraction of the population [see Eq. (1)] catches the disease, and this fraction can be determined using the deterministic framework. If the initial number of infected individuals is very large, the deterministic framework always applies and one can describe the SIR process via the system of differential equations for the densities S(t) of susceptible, I(t) of infected, and R(t) of recovered individuals.
For the critical SIR process, R 0 = 1, the deterministic rate equations readṠ = −SI (A1a)
Equations (A1a)-(A1c) are consistent with the conservation law S(t) + I(t) + R(t) = 1. The initial condition
Below we always assume that
In this range, the number of initially infected individuals, N , is large and the deterministic framework should be asymptotically correct. Setting 1 is not necessary, but this assumption implies that the initial fraction of infected individuals is small; this is natural since we assume that the containment measures making the process critical have started when still a microscopic fraction of the population caught the disease. Furthermore, assuming that 1 allows one to derive much more explicit (asymptotic) formulas than in the general case.
To solve (A1a)-(A1c) subject to (A2) it is convenient to treat S as a time variable. Dividing (A1b) by (A1a) and integrating gives
Similarly dividing (A1c) by (A1a) and integrating yields
The outbreak ends when I f = 0. Equation (A4) shows that the final fraction of susceptible is implicitly determined by
from which
1. Since R f = N −1 E k (N ) and = k/N , we can re-write (A7) as E k (N ) = √ 2kN , leading to the asymptotic Ψ(κ) = √ 2κ when κ 1. The applicability of the deterministic framework for sufficiently large k is intuitively clear. We have determined the crossover k ∼ N 1/3 when the stochastic effects become important using the stochastic framework. It would be interesting to deduce the crossover in the realm of the deterministic framework.
The applicability of the deterministic framework on the final stage with a few remaining infected individuals is questionable. We know, however, that for the critical SIR process starting from a single infected individual, large outbreaks are of the order of N 2/3 . Thus when κ 1, the outbreak size is N 2/3 √ 2κ+O(N 2/3 ). The first term describes the deterministic stage, while the second term accounts for the stochastic final stage. The deterministic first term dominates when κ 1. Therefore in this regime, one can use the deterministic framework.
In Sec. IV, we have used asymptotic methods to determine the leading behaviors. It is possible to derive exact results for the average outbreak size. Unfortunately, these results are expressed through the sum of a large number of variables with complicated individual terms.
Here we explain how one could guess these results; the verification is straightforward, so only the guess part requires a bit of effort. We have made the guess (43) by establishing a few explicit exact results. These exact results describe a non-interesting region where the initial number of susceptible individuals is small, so their virtue is that they simplify the guesswork.
We start by solving (41)-(42) in the simplest cases when N − x = O(1). First, we notice that Similarly we specialize (41) to x = N − 1 and use (B1) to obtain T 1 (i, N − 1) = 1 + N N + 1 T 1 (i − 1, N − 1) + i + 1 N + 1 which is solved subject to T (0, 1) = 0 to yield
Specializing (41) to x = N −2 and using (B2) we arrive at the recurrence
Similarly we compute 
when N 1. The asymptotic behavior of the average size of the outbreak is therefore