Directed Cascades Generate New Critical Universality Classes

Spectral stability determines whether a branching cascade ultimately dies out, but not how large it can become before extinction. This distinction is central to cascade dynamics: a system can be asymptotically stable and yet generate a large transient population capable of reaching rare states. Branching processes provide a natural framework for this problem, with applications ranging from population growth and epidemics to neutron multiplication, neuronal and fracture avalanches, self-exciting point processes, and earthquake triggering [1, 2, 3, 4, 5, 6]. In the scalar case, the mean fertility, or branching ratio, $R$ determines asymptotic behavior, while at criticality the total cluster-size distribution obeys $p(S) \sim1/S^\tau$ with the universal mean-field exponent $\tau= 3/2$ for finite offspring variance [1]. Strong event-to-event fertility fluctuations can change this local universality class. If the expected fertility $\mu$ has a power-law tail $p(\mu) \sim1/\mu^{1+\gamma}$ with $1 < \gamma< 2$, the scalar critical law becomes $p(S) \sim1/S^{1+1/\gamma}$ [7]. In ETAS seismicity, for example, the combination of exponentially increasing productivity $\sigma\propto e^{\alpha m}$ with the Gutenberg–Richter magnitude distribution $p(m) \propto e^{-\beta m}$ produces precisely such a broad fertility distribution, with $\gamma= \beta/\alpha$ [6, 7, 8, 9].

Most natural cascades, however, are multichannel rather than statistically homogeneous. Events may belong to different species, epidemic stages, spatial regions, failure modes, neuronal populations, market sectors, or faulting mechanisms, with type-dependent reproduction and conversion rules. A multitype branching process is described by a non-negative matrix $\mathbf{R}$, where $R_{ij}$ is the mean number of type-$i$ descendants produced directly by one type-$j$ event. Multivariate Hawkes processes are the continuous-time counterpart: the integrated excitation kernels form the same matrix $\mathbf{R}$, and asymptotic stability is set by the spectral condition $\rho(\mathbf{R}) < 1$ [5, 10, 11, 12, 13, 14, 15, 16].

The spectral radius, however, is insensitive to how triggering pathways are organized.

Here we show that this directed architecture is an additional physical control parameter. The main surprise appears when critical components are connected successively: their number along a directed path changes the cascade-size exponent itself. For finite-variance branching, the standard $3/2$ exponent becomes $5/4$ for two successive critical components and $9/8$ for three, generating a hierarchy that approaches $1$ as the number of critical stages increases. The physical mechanism is simple: the entire fluctuating population produced by one critical stage becomes the random number of seeds injected into the next. For broad fertility distributions, the same nesting mechanism generates the corresponding $\gamma$-dependent hierarchy.

Directed architecture also controls amplification away from criticality. Feed-forward pathways can produce polynomial transient growth while leaving the eigenvalues, and hence the asymptotic stability threshold, unchanged. Their depth therefore controls how strongly activity can be amplified before it eventually decays. This separation between spectral stability, transient amplification, and critical scaling applies directly to multitype branching and Hawkes processes, and it can strongly enhance the probability of rare terminal events. In earthquake triggering, we show below that it provides a stationary mechanism for the excess abundance of foreshock–mainshock sequences relative to one-channel ETAS predictions.

The origin of these feed-forward pathways lies in the non-reciprocity of cross-channel excitation: type $j$ may trigger type $i$ much more efficiently than type $i$ triggers type $j$. Grouping mutually interacting types into strongly connected components (SCCs) then separates the recurrent dynamics within components from the directed transfer of activity between them. The SCCs form a directed acyclic graph, yielding the standard Frobenius representation of decomposable multitype branching processes [17, 18, 19]. In Frobenius form, $\mathbf{R} = \mathbf{A} + \mathbf{B}$, where $\mathbf{A}$ is block diagonal and contains the SCCs, while $\mathbf{B}$ contains the strictly feed-forward couplings between them. The spectrum of $\mathbf{R}$ is entirely determined by the SCC blocks in $\mathbf{A}$: $\mathbf{A}$ fixes asymptotic stability, whereas the nilpotent matrix $\mathbf{B}$ leaves the spectrum unchanged and transfers activity downstream.

Two distinct depths then govern the effects of the directed architecture. The feed-forward depth $L$, defined by $\mathbf{B}^{L} \ne\mathbf{0}$, $\mathbf{B}^{L+1} = \mathbf{0}$, is a purely structural property of the feed-forward backbone and sets the highest polynomial order of transient amplification. Critical scaling is governed instead by the critical depth $d(j) \le L + 1$, the maximum number of critical SCCs encountered along a directed path accessible from a seed of type $j$. Thus, $L$ controls the extent of transient amplification, whereas $d(j)$ determines the critical scaling exponent. If $p(S,j)$ denotes the total cascade-size distribution from a seed of type $j$, and $\langle S\rangle_j$ its mean, then $d(j)$ fixes the tail exponent $p(S,j) \sim S^{-\tau_j}$, while the feed-forward structure controls the magnitude of the response.

This role of the critical depth $d(j)$ is quantified by the hierarchy of cascade-size exponents

For finite-variance branching, the standard mean-field exponent $3/2$ is recovered for $d(j) = 1$, but becomes $5/4$ for two successive critical components and $9/8$ for three, approaching $1$ as the critical depth increases. For broad fertility distributions, the corresponding hierarchy is $1/\gamma^{d(j)}$. These results follow from a direct statistical argument. Consider a finite-variance critical branching class, for which the total progeny $S$ has the asymptotic tail

Now consider two consecutive critical classes. Let $S_1$ be the total progeny generated in the upstream class. Since each event in the first class has a finite probability of triggering the second class, the number of downstream seeds is, in first approximation, proportional to $S_1$. Each seed generates an independent critical progeny with the same tail as (2). The sum of $S_1$ such progenies has the characteristic scale

a standard consequence of summing variables with a $1/2$ survival-tail exponent [20, 21]. From (2), the total cascade size $S \sim S_2$ therefore satisfies $\Pr(S > s) \sim\Pr(S_1 > s^{1/2}) \sim s^{-1/4}$, and hence $p(S) \sim S^{-5/4}$. More generally, each additional critical class changes the characteristic scale according to $S_{\ell+1} \sim S_\ell^2$, yielding $\Pr(S > s) \sim s^{-2^{-(d-1)}}$, $p(S) \sim s^{-(1+2^{-(d-1)})}$.

The same argument applies to broad fertilities without modification of its logic. A scalar critical cascade now satisfies $\Pr(X > x) \sim1/x^{1/\gamma}$. The sum of $N$ independent copies has scale $N^\gamma$ [20, 21], so each additional critical channel divides the survival exponent by $\gamma$. This gives $\Pr(S > s) \sim1/s^{1/\gamma^d}$, $p(S) \sim1/s^{1+1/\gamma^d}$, which establishes the physical origin of (1): the complete fluctuations of one branching cascade become the random number of seeds of the next. The generating-function formulation makes the assumptions underlying (1) explicit (see END MATTER).

Numerical simulations (details can be found in Supplemental Material [22]) confirm the predicted scaling for both finite-variance offspring distributions and fertility distributions with tail index $1 < \gamma< 2$. We consider the representative architecture

where columns label parent types and the off-diagonal elements describe the directed transitions $3 \to2$ and $2 \to1$, both with coupling strength $q$. As detailed in Supplemental Material [22], cascades initiated from different seed types probe different critical depths, $d(j)=1,2,3$. Figure 1 shows excellent agreement between the numerical simulations and the theoretical predictions for all seed types and for both fertility regimes.

We now turn to the mean cascade size $\langle S\rangle_j$ and examine how the feed-forward structure enhances it near criticality. For $\rho(\mathbf{R})<1$, the expected total cascade size is $\langle S(j)\rangle=\mathbf{1}^{\mathsf{T}}(\mathbf{I}-\mathbf{R})^{-1}\mathbf{e}_j$. The resolvent admits an exact expansion over directed paths. A path $P:i_0=j\to i_1\to\cdots\to i_r$ contributes

where $B_{i_{a+1}i_a}$ is the mean offspring number associated with the directed transition $i_a\to i_{a+1}$. This expression has a direct physical interpretation: the factors $B_{i_{a+1}i_a}$ quantify how efficiently activity is transmitted from one component to the next, while each factor $1/(1-\lambda_{i_a})$ measures the mean persistence within component $i_a$. Thus, near-critical components act as long-lived amplifying stages, and the directed couplings successively transfer their enlarged populations downstream. Components that remain uniformly subcritical contribute only finite factors and therefore affect the amplitude, but not the critical singularity.

The combinatorics simplify considerably when all cross-channel couplings of the feed-forward matrix $\mathbf{B}$ in the decomposition $\mathbf{R}=\mathbf{A}+\mathbf{B}$ share a common strength $q>0$. To characterize the contribution of the feed-forward structure to the mean cascade size, consider a directed path $P$ containing $\ell(P)$ feed-forward steps, and let $d(P)$ denote the number of SCCs along the path whose distance from criticality is $\epsilon\to0^{+}$. The contribution of such a path has the asymptotic form $W(P)\lesssim C_P\frac{q^{\ell(P)}}{\epsilon^{d(P)}}$, where $C_P$ remains nonsingular as $\epsilon\to0^{+}$. Summing over all paths accessible from a seed $j$ gives $\langle S(j)\rangle=\sum_P W(P)\lesssim\frac{C_j}{\epsilon^{d(j)}}$, where $d(j)=\max_P d(P)$. Thus, the number of near-critical SCCs encountered along a path determines the order of the critical divergence, whereas the nonsingular prefactor $C_j$ collects the contributions from the feed-forward structure and from all paths attaining the maximum critical depth. The structural extent of this amplification is characterized by the feed-forward depth $L=\max_P\ell(P)$, which sets the maximum number of feed-forward steps that can be traversed. Hence, increasing $L$ allows activity to propagate through longer feed-forward pathways and can strongly enhance the amplitude of $\langle S(j)\rangle$, without changing its critical exponent. In this sense, $L$ controls the structural capacity for transient amplification, whereas $d(j)$ determines the singularity of the mean cascade size.

To elucidate the mechanism of transient amplification, we first consider the homogeneous case $\mathbf{A}=\lambda\mathbf{I}$. For a feed-forward chain of $d$ classes, the contribution reaching

the terminal class after $k$ generations is

Thus a chain with $d-1$ feed-forward steps produces a polynomial of degree $d-1$, and the maximal degree is $L$. The eigenvalues determine the eventual exponential decay. This is the characteristic separation between spectral stability and transient response in non-normal systems [23, 24, 25]. In the branching context, the polynomial factors enumerate ordered sequences of type conversions.

Figure 1 confirms the two universality hierarchies. At strict subcriticality, the critical laws persist over an intermediate range whose extent increases upon approaching the transition. For finite-variance offspring, the distribution eventually develops a cutoff. For power law fertility, it crosses over to a faster, fertility-dominated algebraic tail [7]. Directionality nevertheless produces large amplification throughout the subcritical crossover regime.

We now apply this mechanism to earthquake triggering and precursory sequences. The physical basis and empirical construction of a non-normal, mechanism-resolved Hawkes–ETAS model are developed in Ref. [26]. In this description, earthquake types represent faulting mechanisms or other seismologically meaningful triggering channels. Normal, strike-slip, and thrust earthquakes provide three natural classes of triggering events, since both their magnitude statistics and their triggering properties depend systematically on faulting style [27, 28, 29].

Let $R_{ij}$ denote the baseline productivity from channel $j$ (say thrust) to channel $i$ (say strike-slip): a type-$j$ earthquake of magnitude $m$ produces, on average, $R_{ij}^{\mathrm{ETAS}}(m)=R_{ij}e^{\alpha(m-m_{\min})}$ direct descendants of type $i$. Assuming the GR magnitude density $p(m)=\beta e^{-\beta(m-m_{\min})}$, averaging this productivity over magnitudes gives, for $\beta>\alpha$,

At the level of first moments, the existence of a broad distribution of fertilities therefore preserves the entire multitype dynamics after the single replacement $\mathbf{R}\mapsto\overline{\mathbf{R}}^{\mathrm{ETAS}}=y\mathbf{R}$: generation profiles, resolvents, and transient polynomial amplification retain the form of homogeneous fertilities (no marks) This closure uses only the mean mark factor $y$. Critical fluctuations are more subtle and cannot in general be obtained by the same scalar renormalization. The complete fertility distribution fixes the leading nonlinear singularity of the offspring generating function. Here the productivity law and the GR distribution generate a power law fertility tail with $\gamma=\beta/\alpha$. Consequently, $\rho(y\mathbf{R})$ determines asymptotic stability, the directed architecture fixes the nesting depth, and $\gamma$, rather than $y$ alone, changes each nested square root into a $\gamma$th root. This separates the structural amplification mechanism from the fertility-dependent critical behavior.

In the conventional one-channel ETAS model [6, 9], earthquake magnitudes are independent GR draws within each triggering tree, so foreshocks arise when a small event initiates a cascade that later reaches a larger magnitude. Although this mechanism generates foreshock sequences [30, 31], standard one-channel ETAS simulations substantially underpredict their observed abundance [32, 33, 34]. Directed multichannel branching provides a stationary mechanism for this excess: it amplifies the population of descendants, and therefore the number of opportunities for a cascade to reach the upper magnitude tail, without requiring time-dependent preparation or evolving model parameters. Such multichannel triggering is physically natural in heterogeneous fault systems. Distinct faulting mechanisms can coexist and interact within the same tectonic environment because deformation generally cannot be accommodated by a single mode alone [26, 35]. A population generated in one faulting channel can therefore seed another, whose fluctuating population may in turn activate a further channel. Crucially, what is transmitted downstream is not a representative event but an entire random population. The cascade-size fluctuations produced in one channel thus become the random number of triggering opportunities supplied to the next, providing the physical origin of the amplification.

Figure 2(a) shows one realization of this process with the matrix in (4) with $\lambda_1=\lambda-q$. In this case, increasing $q$ is compensated by reducing part of the diagonal reproduction, so that the total column sum is independent of $q$. The cascade begins in type 3, whose events reproduce over several generations and eventually seed type 2. The resulting type 2 population subsequently seeds type 1, which continues to generate events over later generations. The three populations overlap in time: downstream activation does not require the upstream population to have become extinct. Moreover, type 1 activity can reappear at late generations because different branches of the genealogy reach the downstream channel after different delays. Panel (a) thus makes explicit the mechanism underlying the amplification: successive random populations are transmitted through the directed sequence $3 \to2 \to1$.

This transmission affects not only the mean number of descendants but the complete cascade-size distribution. It consequently increases the number of opportunities to sample the upper GR magnitude tail. To quantify this effect, we initiate a cascade with a type 3 source earthquake of magnitude $m_0=m_{\min}$ and define

where $m_{\max}$ is the largest magnitude reached anywhere in the resulting multichannel cascade.

Figure 2(b) shows that $Q_{m_{\mathrm{target}}}(q)/Q_{m_{\mathrm{target}}}(0)$ increases strongly with directed coupling $q$, even though the spectral radius and column count are kept fixed. The relative enhancement is larger for the more extreme target magnitudes. Panels (a) and (b) therefore connect the microscopic mechanism to its statistical consequence: successive transmission through distinct triggering channels generates a larger and more broadly fluctuating descendant population, which raises the probability of reaching the remote magnitude tail. This result does not constitute a prediction of individual earthquakes. It provides a population-level mechanism that increases the frequency of sequences classified retrospecttively as foreshock–mainshock sequences and gives an explanation for the high foreshock rates observed in empirical catalogs relative to standard one-channel ETAS predictions [32, 33, 34].

The same logic applies beyond seismicity. Whenever rare extreme events are reached through several asymmetric reproduction or conversion stages, a one-channel model underestimates the probability of large terminal events. Examples include progression between epidemic stages or host populations, mutation and phenotypic differentiation, successive failure modes in heterogeneous materials, activation across neuronal populations or cortical layers, reaction and particle-multiplication chains, and contagion across economic sectors. In each case, precursor abundance may reflect the directed organization of otherwise stationary subcritical dynamics rather than a temporal drift toward instability.

We acknowledge CINECA for computational resources through the ISCRA C project IsCc7. DS acknowledges support from the Center for Computational Science and Engineering at Southern University of Science and Technology. EL acknowledges funding from the Italian Ministero dell’Università e della Ricerca under PRIN 2022 (“re-ranking of the final lists”), project 2022KWTEB7, CUP B53C24006470006. GP acknowledges support from the DyDéCo project and from the CPJ Chair at ENS de Lyon and CNRS.

End Matter

Consider a directed chain of critical strongly connected components (SCC) labelled $a=1,\ldots,d$, with triggering from class $a+1$ to class $a$. Class $a=1$ is the terminal downstream class. Let $p_a(S)$ be the probability that a cascade initiated by one event in class $a$ contains $S$ events in total, including the initial event and all descendants in class $a$ and in the downstream classes. Its probability generating function is

The large-$S$ behavior of $p_a(S)$ is determined by the singularity of $H_a(z)$ as $z\to1^-$. We therefore introduce

so that $x\to0^+$ probes the large-cascade limit.

We must distinguish $H_a$, which describes the complete cascade, from the generating function of the direct offspring of one event. For a parent in class $a+1$, let $K_{a+1}$ be the number of direct offspring remaining in class $a+1$, and let $L_{a+1,a}$ be the number transmitted to the next downstream class $a$. Their joint offspring generating function is

The corresponding mean offspring numbers are

Critical reproduction within class $a+1$ means $\langle K_{a+1}\rangle=1$, while directed transmission requires $q_a>0$.

For finite, nonzero offspring variance, the expansion of the direct offspring generating function around $(s,t)=(1,1)$ can be written as

where

The last equality follows from $\langle K_{a+1}\rangle=1$. The coefficients $C_{a+1}$ and $D_{a+1}$ describe mixed and downstream offspring fluctuations, respectively.

The exact tree equation for a cascade initiated in class $a+1$ is

Substituting $z=1-x$, $H_{a+1}=1-\epsilon_{a+1}$, and $H_a=1-\epsilon_a$ into Eq. (14), and then using the expansion (12), gives

The terms $-\epsilon_{a+1}$ on the two sides cancel precisely because the mean same-class reproduction is critical. Retaining the dominant terms then yields

For the terminal critical class, $A_1\epsilon_1^2\sim x$, and hence $\epsilon_1\asymp x^{1/2}$. Since $\epsilon_1\gg x$ as $x\to0^+$, (16) implies recursively

Thus every additional critical SCC takes one further square root of the downstream singularity, giving

Standard singularity analysis [36] then gives $p_d(S)\sim S^{-1-2^{-d}}$.

If the fertility distribution instead has a power-law tail with $1<\gamma<2$, its variance diverges and the leading nonlinear term in Eq. (12) is $A_{a+1}u^\gamma$ rather than $A_{a+1}u^2$. The same substitution gives

so that each additional critical SCC takes a $\gamma$th root: $\epsilon_d\asymp x^{\gamma^{-d}}$. Consequently, $p_d(S)\sim S^{-1-\gamma^{-d}}$. For a general seed of type $j$, $d$ is replaced by the largest accessible critical depth $d(j)$. The coupling coefficients determine the prefactors, whereas the critical depth and the local fertility distribution determine the tail exponent. The complete calculation is presented in the Supplemental Material [22].

Expected total progeny

Consider a multitype branching process with $n$ types. Let $R_{ij}$ be the mean number of offspring of type $i$ generated by one individual of type $j$. The process is subcritical whenever $\rho(\mathbf{R}) < 1$, where $\rho(\mathbf{R})$ denotes the spectral radius. If the cascade starts from one individual of type $j$, the expected total progeny (including the ancestor) is exactly

where $\mathbf{e}_{j}$ is the $j$th canonical basis vector and $\mathbf{1}$ denotes the vector whose entries are all equal to one. Equation (S1) is exact for every finite-type Galton-Watson branching process with finite first moments.

Frobenius decomposition and finite resolvent expansion

As in the Letter, we first group mutually reachable types into strongly connected components (SCCs) and order the SCCs topologically. The offspring matrix then has the Frobenius decomposition

where $\mathbf{A}$ is block diagonal and contains the internal dynamics of the SCCs, whereas $\mathbf{B}$ contains only feed-forward couplings between distinct SCCs. It is strictly block triangular and therefore nilpotent. We use exactly the convention of the Letter: the feed-forward depth $L$ is defined by

Starting from $\mathbf{I} - \mathbf{R} = \mathbf{I} - \mathbf{A} - \mathbf{B}$, we factorize

Indeed,

Taking the inverse gives

Using the Neumann identity

with $\mathbf{X} = (\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$, we obtain

Multiplication by the block-diagonal matrix $(\mathbf{I} - \mathbf{A})^{-1}$ preserves the strictly block-triangular structure. Hence the expansion terminates exactly at $r = L$. The SCC blocks in $\mathbf{A}$ determine the spectrum and asymptotic stability, whereas successive factors of $\mathbf{B}$ enumerate feed-forward transfers between SCCs.

Expansion over directed paths

For explicit scalar path weights, specialize each SCC to a single type, so that

Equation (S6) admits a simple genealogical interpretation. The diagonal resolvent satisfies

which represents the amplification produced by repeated reproduction inside type $i$. Each factor of $\mathbf{B}$ corresponds instead to one directed conversion between different types.

Consider an arbitrary directed path

Its contribution to the expected cascade size is

The total expected progeny is therefore obtained by summing over every directed path accessible from the initial type,

Equation (S10) is exact and provides a complete path representation of the resolvent.

Homogeneous conversion chains

For homogeneous conversion weights, we write

where $\mathbf{G}$ is the binary adjacency matrix of the directed feed-forward graph. Suppose now that every allowed conversion has the same weight,

For a path containing $\ell$ directed conversions,

(129) becomes

The numerator depends only on the number of directed conversions, whereas the denominator depends only on the distance from criticality of the classes visited by the path.

Homogeneous chain of critical SCCs

The homogeneous multicritical limit used in the Letter corresponds to

Every class visited by the longest directed path approaches criticality simultaneously. If this path contains $d$ critical SCCs, the exact path weight gives

where

counts the total weight of maximal critical paths accessible from the seed. This is the homogeneous specialization of the path contribution stated in the Letter. Here $d$ is the number of critical SCCs on this path, while the number of feed-forward links is $d-1$.

Range of validity

(133) does not hold for an arbitrary non-normal branching matrix. Rather, it is the homogeneous specialization of the exact path expansion (S10). The correct general statement is therefore:

For a decomposable multitype branching process, the expected total progeny admits an exact expansion over directed paths. Each path contributes the product of its conversion weights divided by the product of the distances from criticality of the classes visited. Consequently, if the initial type has access to a chain of $d$ successive classes that become critical simultaneously, the expected cascade size diverges as $(1-\lambda)^{-d}$. For homogeneous conversion chains this reduces to $S \sim q^{d-1}(1-\lambda)^{-d}$.

Noncritical segments and critical depth

Assume that a path $P$ visits $d(P)$ classes approaching criticality, $1-\lambda_i=\epsilon$, while the remaining $m$ classes remain uniformly subcritical, $1-\lambda_i=\epsilon+q$. Since the path contains $n=d(P)+m$ classes, its length is $\ell=n-1$. Its contribution is therefore

(134) is exact.

Critical scaling

Suppose first that $q$ remains finite while $\epsilon\to0$. Then

and (134) becomes

Since $d(P)=n-m$, this simplifies to

Hence, every non-critical class contributes only a finite resolvent factor, while each critical class contributes one singular factor $(1-\lambda)^{-1}$. Consequently, the critical divergence depends only on the number $d(P)$ of simultaneously critical SCCs visited by the path.

Maximizing over all paths accessible from a seed of type $j$ defines the critical depth used in the Letter,

This must not be confused with $L$: the feed-forward depth depth counts all feed-forward links, whereas $d_*(j)$ counts only critical SCCs on an accessible path. When $d(P)=1$, all noncritical segments contribute only finite factors and the path diverges as $(1-\lambda)^{-1}$.

Matrix powers and transient amplification

For a homogeneous scalar chain, write $\mathbf{R}=\lambda\mathbf{I}+q\mathbf{G}$, where $\mathbf{G}^{L}\ne0$ and $\mathbf{G}^{L+1}=0$. Since $\lambda\mathbf{I}$ and $\mathbf{G}$ commute,

Along the longest feed-forward chain, the dominant contribution is

For the three-type chain $3 \to2 \to1$ used in the Letter, $L=2$ and the homogeneous case gives

The off-diagonal terms contain polynomial prefactors multiplying the exponential decay. In particular,

whereas

Although every matrix element eventually decays because $\lambda<1$, these polynomial corrections produce a finite-generation maximum before the asymptotic exponential decay sets in.

In particular, there exists a value $k_{\max}$ such as off-diagonal elements present a maximum value. The calculation of $k_{\max}$ and the maximum value $y_{\max}$ is immediate in the $\lambda=\lambda_{1}$ case. Taking for instance $a_{13}(k)=\frac{1}{2}q^{2}k(k-1)\lambda^{k-2}$, we obtain

and

For $\lambda\ne\lambda_{1}$, the analytical expressions become more cumbersome.

The matrix exhibits the same qualitative behavior: the off-diagonal elements initially increase, reach a maximum, and finally decay exponentially. Therefore, transient amplification is a generic consequence of directed propagation pathways rather than a peculiarity of the degenerate case.

The calculation of $k_{\max}$ in this case $\lambda_{1}\ne\lambda$ is more complicated. Considering $a_{13}(k)=\frac{q^{2}}{(\lambda-\lambda_{1})^{2}}\left[k-k\lambda\lambda_{1}^{k-1}+(k-1)\lambda_{1}^{k}\right]$, we find that $k_{\max}$ is the solution of the transcendental equation

The distribution of cluster size $p(S,j)$

Generating functions

Let $H_{j}(z)$ denote the generating function of the total progeny produced by a cascade initiated in critical class $j$,

where $p(S,j)$ is the probability that the total cascade size equals $S$. For a directed critical chain, we label the classes by $a=1,\ldots,d$, starting from the terminal downstream class $a=1$. The index $j$ is instead reserved for a generic seed type. Near the singular point $z=1$, we introduce

The asymptotic behavior of $\epsilon_{a}(x)$ determines the tail of the cascade-size distribution through standard singularity analysis.

Joint offspring generating function and its expansion

The total-cascade generating function must be distinguished from the generating function of the direct offspring of one event. For an event in class $a+1$, let $K_{a+1}$ be the number of direct offspring that remain in class $a+1$, and let $L_{a+1,a}$ be the number transmitted to the next downstream class $a$. Their joint offspring generating function is

Its first derivatives at $(1,1)$ are the corresponding mean offspring numbers. Critical same-class reproduction and nonzero feed-forward coupling therefore mean

Assuming finite second moments, we set $s=1-u$ and $t=1-v$ and expand around $(1,1)$:

The curvature in the same-class direction is

where the last equality uses $\left\langle K_{a+1}\right\rangle=1$. The coefficients $C_{a+1}$ and $D_{a+1}$ encode mixed and downstream offspring fluctuations. No Poisson assumption is required.

Terminal critical class

Consider first the terminal critical class. Its generating function satisfies

Substituting $H_1 = 1 - \epsilon_1$, $z = 1 - x$, and using the one-variable specialization of (125) gives $1 - \epsilon_1 = (1 - x)(1 - \epsilon_1 + A_1\epsilon_1^2 + o(\epsilon_1^2))$. Expanding to leading order yields $x \sim A_1\epsilon_1^2$, or equivalently

This is the classical square-root singularity responsible for the well-known $S^{-3/2}$ critical branching law.

Recursive critical classes

Suppose now that critical class $a+1$ feeds class $a$ through a non-zero coupling. The exact tree equation is

The first argument accounts for offspring remaining in class $a+1$; the second accounts for offspring transmitted to class $a$, each of which generates its own complete downstream cascade. Substitution of $z = 1 - x$ and $H_b = 1 - \epsilon_b$ into (129), followed by use of (125), cancels the term linear in $\epsilon_{a+1}$ because same-class reproduction is critical. The leading balance is

The omitted terms $\epsilon_{a+1}\epsilon_a$, $\epsilon_a^2$, and $x\epsilon_{a+1}$ are asymptotically smaller. Since (128) implies $\epsilon_1 \gg x$, induction gives $\epsilon_a \gg x$ at every subsequent level. Hence,

Each additional critical feed-forward class therefore takes one further square root of the singularity generated downstream.

Induction

Starting from $\epsilon_1 \sim x^{1/2}$, successive application of the previous relation gives $\epsilon_2 \sim x^{1/4}$, $\epsilon_3 \sim x^{1/8}$, and by induction,

Therefore each additional critical class halves the singularity exponent.

Singularity analysis

By the Transfer Theorem of Flajolet and Sedgewick [36], if $1 - H(z) \sim(1 - z)^\alpha$, then, using $\Gamma(-\alpha) = -\frac{\Gamma(1-\alpha)}{\alpha}$, the coefficients satisfy

In our specific case

where $C_d$ is a non-universal constant. Consequently,

with

This establishes the finite-variance hierarchy stated in the Letter.

Conditions of validity

The above derivation relies only on the following assumptions:

1. the process is decomposable, or equivalently the critical classes can be arranged in feed-forward order;

2. each critical class possesses finite, strictly positive offspring variance;

3. successive critical classes are connected by non-zero coupling;

4. the entire chain is accessible from the initial type;

5. all classes along the chain become critical simultaneously;

6. the offspring generating functions are analytic around the critical point.

Remarkably, Poisson statistics are never required.

Several accessible paths

Suppose now that several directed critical chains are accessible from the initial seed. Let $d(P)$ denote the number of critical SCCs visited by path $P$; uniformly subcritical SCCs on the same path are not counted. The dominant asymptotic behavior is governed by the deepest accessible critical chain,

provided at least one maximal-depth path has non-zero weight. Therefore

All shorter chains contribute only subleading corrections.

Universality statement

The previous derivation establishes the following result.

Universality theorem.

For a decomposable critical multitype Galton–Watson branching process with finite, non-zero offspring variances, the total progeny distribution generated by a seed having access to a chain of $d$ successive critical classes obeys

The exponent depends exclusively on the number of successive critical classes along the longest accessible feed-forward chain. Offspring distributions, coupling strengths, and microscopic details modify only the non-universal prefactor.

(136) therefore represents a universality law for decomposable finite-variance multitype branching processes. It should not be interpreted as applying to arbitrary critical multitype branching matrices. In particular, irreducible critical multitype branching processes generally recover the ordinary mean-field exponent $\tau= \frac{3}{2}$, whereas the hierarchy $\tau_d = 1+2^{-d}$ originates from the successive nesting of critical feed-forward classes.

(136) is therefore universal only within the class of finite-variance decomposable critical multitype branching processes. The distinction is structural. In a decomposable process, the critical classes can be ordered so that activity passes only downstream from one class to the next. Each additional critical class in the longest accessible chain adds another level of critical nesting, giving

For example, two critical classes coupled only in one direction may have a mean offspring matrix of the form

A seed that can access both critical classes has $d=2$ and hence $\tau_2 = 5/4$. By contrast, an irreducible critical matrix such as

has mutually communicating types and constitutes a single critical class; it therefore exhibits the ordinary finite-variance mean-field exponent $\tau= 3/2$. Thus the exponents $\tau_d < 3/2$ for $d > 1$ are not a generic property of critical multitype branching, but arise specifically from a feed-forward chain of successive critical classes.

Crossover from a decomposable to an irreducible critical process

Let us introduce a weak backward coupling $\eta> 0$ while keeping the matrix critical. For this, we retune the diagonal reproduction rates and write

The two eigenvalues are

so that $\rho(\mathbf{R}_{\eta}) = 1$ for every $\eta> 0$. At the same time, $\mathbf{R}_{\eta}$ is irreducible for every $\eta> 0$, because the two classes communicate in both directions. Thus $\eta= 0$ is a singular limit: exactly at $\eta= 0$ the system contains two successive critical classes, whereas for any $\eta> 0$ it contains a single irreducible critical class.

To estimate the crossover scale, let $\epsilon= 1-z$, and denote by $x$ and $y$ the deficits of the total-progeny generating functions from unity. Expanding the offspring generating functions about the critical fixed point, finite non-zero offspring variances give, to leading order,

where $c_1,c_2 > 0$ are non-universal constants determined by the second moments of the offspring distributions. Only the powers of $\epsilon$ are needed for the following scaling argument.

In the strictly decomposable limit $\eta= 0$, hence $a=0$, (143) gives

Substituting this result into (142) yields

and therefore

The singularity $1-H(z) \sim(1-z)^{1/4}$ implies

as expected for a chain of two critical classes.

For $0 < \eta\ll1$, the decomposable scaling remains valid as long as the terms introduced by the critical retuning are asymptotically negligible. In (143), the first such correction is the linear term $ay$. Using the decomposable estimate $y \sim\varepsilon^{1/2}$, this term becomes comparable to the leading term $\varepsilon$ when

Thus, the reducible fixed-point scaling is observed for $\varepsilon\gg a^2 = q\eta$. In contrast, sufficiently close to $z = 1$, the weak backward coupling is resolved and the process crosses over to the irreducible critical behavior.

Coefficient asymptotics probe the singularity on the scale $1-z \sim S^{-1}$. The cross-over occurs for $\varepsilon_\times= a^2 = q\eta$, which corresponds to a cluster-size crossover

up to a non-universal multiplicative constant. For fixed $q > 0$, the extent of the intermediate decomposable regime consequently diverges as

The resulting crossover for the distribution of cluster sizes has the scaling form

Thus an arbitrarily weak backward coupling does not immediately erase the $5/4$ law. Instead, after retuning the process to criticality, the feed-forward exponent survives as an intermediate asymptotic over a range of cluster sizes that grows as an inverse power of the feedback strength. The true $S \to\infty$ asymptotic exponent is nevertheless $3/2$ for every fixed $\eta> 0$, because the two classes then form a single irreducible critical class.

Appendix D: The total progeny in the ETAS model

In this section we derive the marked counterpart of the mean-resolvent formula used in the Letter. The derivation closely parallels that of the unmarked branching process, the only additional ingredient being the magnitude-dependent fertility of each event.

Marked multitype ETAS process

Consider a multitype ETAS model with $n$ earthquake types. The mean number of direct aftershocks generated by an earthquake of magnitude $m$ and type $j$ is assumed to factorize into a magnitude-dependent productivity and a type-transition matrix,

where

is the standard ETAS productivity law, $\mathbf{R}$ is the multitype triggering matrix, and $\mathbf{e}_j$ denotes the unit vector corresponding to the initiating type.

Throughout this section, we assume that descendant magnitudes are independently drawn from the Gutenberg–Richter distribution,

so that the average productivity is finite.

Average productivity

The mean productivity of an earthquake is

Substituting Eqs. (S2) and (S3) gives

Consequently, the average offspring matrix for every generation after the initiating earthquake is simply

At the level of first moments, the marked ETAS process therefore differs from the corresponding unmarked branching process only through the multiplicative factor $y$. This reduction does not determine higher moments or the singular part of the total-progeny distribution, which retain information about the full mark-dependent fertility law.

Expected cascade size

Consider an initiating earthquake of fixed magnitude $m_0$ and type $j$. Generation zero consists only of the initiating event, $N^{(0)} = 1$. The expected number of first-generation aftershocks is

After the first generation, magnitudes are independently sampled from the Gutenberg–Richter distribution. Consequently, each subsequent generation acquires an average productivity factor $y$.

The expected number of second-generation events is therefore

Proceeding recursively,

(130) is the marked analogue of the generation expansion of an ordinary multitype branching process.

Exact resolvent formula

The expected total number of earthquakes generated by the initiating event, including the root earthquake itself, is

Substituting (130) gives

Extracting one factor of $\mathbf{R}$,

Whenever $\rho(y\mathbf{R})<1$, the geometric series converges,

leading to the exact expression

(135) is the exact marked counterpart of the resolvent formula

derived for the unmarked process.

Homogeneous nilpotent chain

Consider now the homogeneous decomposition

where, consistently with (123), $\mathbf{B}^{L}\ne\mathbf{0}$, $\mathbf{B}^{L+1}=\mathbf{0}$. We write the nilpotent feed-forward matrix as

where $\mathbf{G}$ is a strictly triangular adjacency matrix satisfying $\mathbf{G}^{L}\ne\mathbf{0}$, $\mathbf{G}^{L+1}=\mathbf{0}$. Therefore,

Factoring out the diagonal term,

Using again the Neumann expansion,

and exploiting the nilpotency of $\mathbf{G}$,

Substituting into (135) yields

Three-type chain

For the elementary chain $3\to2\to1$, the structural depth is $L=2$. Since $\mathbf{G}^{3}=\mathbf{0}$, (141) reduces to

If the cascade starts from the upstream type ($j=3$) and its initiating magnitude is also averaged over the Gutenberg–Richter law, one obtains

This expression is exact. For a fixed initiating magnitude $m_0$, the corresponding formula contains the root-dependent prefactor derived below in (108).

Equivalent resolvent representation

The expected total number of earthquakes generated by an initiating event of type $j$ and magnitude $m_0$, including the initiating event itself, is

Using (130) we obtain

Extracting one factor of $\mathbf{R}$ gives

Since $\rho(y\mathbf{R})<1$,

and therefore

Introducing the magnitude-averaged ETAS offspring matrix

we use the matrix identity

We immediately obtain

Hence, the exact resolvent formula can be written as

Equation (108) is the exact marked analogue of

for ordinary multitype branching processes. The only modification is the replacement of the branching matrix $\mathbf{R}$ by its effective marked counterpart $y\mathbf{R}$, together with the prefactor associated with the initiating earthquake.

Critical scaling

We now investigate the behavior of (108) close to the ETAS critical point,

For the homogeneous decomposition

(141) gives

Substituting this expression into (108) yields

In this homogeneous critical limit, let $d_*(j)$ be the number of critical SCCs on the deepest path accessible from the initiating type $j$. Such a path contains $d_*(j)-1$ feed-forward links, and

The leading singular contribution accessible from $j$ is therefore

Introducing the geometric factor

we obtain

Equation (115) is the fixed-root, marked analogue of the homogeneous path scaling in the Letter. The only modification with respect to the unmarked process is the replacement

together with the multiplicative factor $a(m_0)y^{d_*(j)-2}$, which depends on the productivity of the initiating earthquake but leaves the critical exponent unchanged.

Averaging over the initiating magnitude

The previous expression refers to an initiating earthquake with fixed magnitude $m_0$. If the initiating magnitude itself is assumed to be drawn from the Gutenberg–Richter distribution, its average productivity is

Therefore

Taking the average of (108) over the initiating magnitude immediately gives

Hence

Remarkably, after averaging over the initiating magnitude, the ETAS expression becomes formally identical to the resolvent formula of the unmarked multitype branching process, the only difference being the replacement of the branching matrix $\mathbf{R}$ by its effective counterpart $y\mathbf{R}$.

General ETAS scaling theorem

The previous derivation establishes the following result.

ETAS mean-scaling theorem. Consider a decomposable multitype ETAS process with Gutenberg–Richter distributed magnitudes and productivity law $a(m)=e^{\alpha(m-m_{\min})}$. If an initiating event of type $j$ has critical depth $d_*(j)$, the expected total number of events generated by an earthquake of magnitude $m_0$ obeys

where

and $C_j$ depends only on the topology of the accessible directed chains. After averaging over the initiating magnitude, it holds that

whose leading critical divergence is

These statements concern the mean cascade and its critical divergence. For all first-moment dynamical quantities, magnitude averaging renormalizes the interaction matrix according to $\mathbf{R}\mapsto y\mathbf{R}$, and hence moves the critical point from $\lambda=1$ to $y\lambda=1$, together with a non-universal amplitude factor. The distributional critical behavior is not fixed by $y$: it depends on the leading nonlinear singularity generated by the complete fertility distribution, as derived in the next section.

Cascade-size distribution in the ETAS model

We now extend the previous generating-function derivation to the magnitude-dependent ETAS process. This step exposes the separation that is invisible at the level of the mean: the graph and its critical depth remain structural properties, whereas the mark distribution selects the local branching singularity that is recursively nested along the graph.

As in the finite-variance case, we label the classes along a directed critical chain by $a=1,\ldots,d$, with $a=1$ denoting the terminal downstream class. We define

and introduce

For the marked ETAS process, the offspring generating function has the non-analytic expansion

For the directed chain $3\to2\to1$, the offspring generating functions become

The total-tree generating functions satisfy

For type 1,

which immediately gives

Hence

recovering the scalar marked branching law

For type 2,

Near $z=1$,

so that $\epsilon_2 \gg\epsilon_1$. Consequently,

yielding

with

Therefore,

Proceeding identically for type 3,

which gives

where

Applying standard singularity analysis,

The scalar marked-process exponent,

is therefore replaced by the non-normal exponent

Since $1<\gamma<2$, one has

showing that directed non-normal amplification systematically produces substantially heavier critical tails than the scalar marked process. More generally, a seed with critical depth $d_\ast(j)$ has tail exponent $\tau_j=1+\gamma^{-d_\ast(j)}$. Thus the mean dynamics depends on marks only through $y$, while the critical cascade-size law depends on the full fertility tail through $\gamma$; the feed-forward architecture controls how many times that local singularity is composed.

Numerical simulations

We test the analytical predictions by direct Monte Carlo simulations of both the unmarked multitype branching process and the marked ETAS process. Each realization starts from a single ancestor of prescribed type $j_0$. For the unmarked process, a parent of type $j$ generates a Poisson-distributed number of offspring of type $i$ with mean $R_{ij}$. For the marked ETAS process, a parent of type $j$ and magnitude $m$ generates a Poisson-distributed number of offspring of type $i$ with mean $R_{ij}e^{\alpha(m-m_{\min})}$. Each offspring magnitude is independently drawn from the Gutenberg–Richter distribution

The branching procedure is iterated generation by generation until the cascade becomes extinct. For each realization, we record the total cascade size $S$, including the initial ancestor.

We use the three-type matrix

where the off-diagonal elements describe a preferential propagation pathway $3\to2\to1$.

We investigate two regimes. In the compensated case, we set

so that each row has the same total weight $\lambda$. In the homogeneous multicritical case, we set

In both cases, $q\ge0$ controls the strength of the feed-forward coupling. In all cases

We always assume that the magnitude of the initial earthquake is set to $m_0=m_{\min}$ and indicate with $j_0$ its type. To test these theoretical predictions, we fix the Gutenberg–Richter exponent to $b=1$ and consider three different productivity exponents, $\alpha=0.77,\ 0.81,\ 0.83$, corresponding to $\lambda=1.19,\ 0.98,\ 0.90$, chosen so that the effective branching ratio remains close to criticality ($B_r\simeq1$).

For each value of $\alpha$, we consider four values of the coupling parameter, $q=0,\ 0.1,\ 0.2,\ 0.3$, keeping $\lambda_1=\lambda$. For the largest coupling ($q=0.3$), we additionally investigate two reference cases:

• $\lambda_1=\lambda-q$ with $j_0=3$;

• $\lambda_1=\lambda$ with $j_0=2$.

Figure S1–Figure S3 compare the Monte Carlo simulations with the analytical predictions. For all three values of $\alpha$, the numerical distributions are found to be in excellent agreement with the predicted power law behavior. In particular, the simulations clearly reproduce the three different scaling exponents predicted by the theory,

At the same time, we recover the standard isotropic exponent when $\lambda_1=\lambda-q$ with $p_3(S)\simeq3p_{\mathrm{scalar}}(S)$, where $p_{\mathrm{scalar}}(S)$ denotes the corresponding scalar-process distribution.

Figure S4 explores the case $\gamma>2$ and we observe that the theoretical scaling behaviors

are recovered in numerical simulations.