Directed percolation in random temporal network models with heterogeneities

The event graph representation of temporal networks suggests that the connectivity of temporal structures can be mapped to a directed percolation problem. However, similarly to percolation theory on static networks, this mapping is valid under the approximation that the structure and interaction dynamics of the temporal network are determined by its local properties, and, otherwise, it is maximally random. We challenge these conditions and demonstrate the robustness of this mapping in case of more complicated systems. We systematically analyze random and regular network topologies and heterogeneous link-activation processes driven by bursty renewal or self-exciting processes using numerical simulation and ﬁnite-size scaling methods. We ﬁnd that the critical percolation exponents characterizing the temporal network are not sensitive to many structural and dynamical network heterogeneities, while they recover known scaling exponents characterizing directed percolation on low-dimensional lattices. While it is not possible to demonstrate the validity of this mapping for all temporal network models, our results establish the ﬁrst batch of evidence supporting the robustness of the scaling relationships in the limited-time reachability of temporal networks.


I. INTRODUCTION
Connectivity is an essential characteristic of complex networks as it determines how far information or influence can spread in a network structure.Consequently, it governs the emergence and scale of any macroscopic phenomena often modelled on networks such as disease spreading, transportation, or information diffusion, to mention a few examples.Percolation theory provides a comprehensive understanding that characterizes network connectivity with various mathematical and algorithmic tools primarily developed for complex networks.For example, percolation can be mapped to late-stage results of specific epidemic processes [1][2][3][4][5], such that the size of percolating components determine the final size of the epidemic.Meanwhile, the percolation transition and its related critical behavior explain the disease outcome close to the epidemic threshold.
However, these theoretical descriptions commonly assume that the network is static, with links and nodes always present, ignoring the typical character of several complex structures where links may vary in time.Since information, disease or other effects can pass between two nodes in a network only at the time of their interactions, the temporal alternation of links may crucially influence the critical behavior and final outcome of any ongoing spreading processes [6][7][8][9][10][11][12].To characterize these processes, one needs to measure connectivity in temporal networks across time, where components are defined in terms of network nodes and links and the temporal distribution of interactions.Consequently, beyond the wellstudied structural heterogeneities of static networks, like in their node degrees, the effects of temporal correlations leading to temporal heterogeneities in the interaction dynamics, like burstiness, become important [13][14][15][16][17][18].This is especially the case for so-called limited-waiting-time processes, where an effect or information, e.g. a disease or a meme [19], arriving at a node can pass over to another node only if an interaction appears within a time window δt.Otherwise, the pathogen times out, e.g., the patient recovers or the meme becomes irrelevant, making it impossible to reach other nodes.
Similar to static networks, the connectivity of temporal networks passes through a phase-transition.However, close to this critical threshold, temporal networks exhibit different critical behavior as compared to static structures [20][21][22].For limited-waiting-time connectivity, where the control parameter is δt, this phase-transition can be theoretically understood under some simplifying assumptions about the homogeneous dynamics of connectivity [22].Since there is an embedded direction (or flow) of time, the microscopic dynamics can be fundamentally irreversible with a broken detailed balance and non-equilibrium steady-state.These results suggest that the dynamics of percolation on temporal networks are generically the same as any other system belonging to the Directed Percolation (DP) universality class, which is characterized by a one-component order parameter without additional symmetries and unconventional features such as quenched disorder [23].
The homogeneity approximations used for the derivations presented in [22], however, become less grounded when the underlying structure deviates from a random graph or if the interaction dynamics become inhomogeneous.In this paper, our goal is to build on the theory laid down in [22] to investigate further the relation between temporal networks and directed percolation.In other words, the primary objective of this manuscript is as follows: to show empirically that diverse classes of temporal networks, with various degrees of temporal and spatial heterogeneity, combined with the very general notion of limited-waiting-time reachability, will show an absorbing phase transition in connectivity that belongs to the directed percolation universality class.
In its epidemic interpretation, directed percolation can be one of the most basic non-equilibrium second-order phase transitions from fluctuating states into so-called absorbing states, which exhibit universal features, determined by symmetry properties and conservation laws.We demonstrate the precision of this mapping using extensive numerical simulations and provide further theoretical calculations to study synthetic temporal networks as directed percolation processes with a range of temporal and spatial inhomogeneities.
The remainder of Sec.I will be dedicated to laying the groundwork and presenting the context in which this manuscript is set: in Sec.I A we will discuss connectivity on temporal networks, the event graph representation and modeling spreading processes and Sec I B will introduce directed percolation and its characteristics.
The next section, Sec.II, is dedicated to an overview of our contributions.Section II A will describe our mapping of concepts of directed percolation and temporal networks.Section VI provides an overview of the theoretical results from Ref. [22], which will be extended further in Sec.VI 1 by explicitly deriving some critical exponents and scaling relations.Section II B will lay down the algorithmic techniques that make large-scale simulations of spreading processes on temporal networks possible.
Finally, in Sec.III we will describe the experimental setup and provide numerical evidence for validity of our hypothesis by application of the methods described previously, while Sec.IV provides an overview of the implications of the results and the limitations of our study.

A. Temporal networks and the event graph
A temporal network G = (V, E) provides representations of a dynamically changing complex system as a set of timed interactions known as events E between a set of entities V = {v 1 , v 2 , . . ., v n } known as nodes or vertices during an observation period T .Each event indicates a time-dependent interaction between two nodes, e.g.physical contact or communication between two people or trade between two commercial entities [24], i.e. e = (u, v, t start , t end ) such that u, v ∈ V between times of t start , t end ∈ T (t start < t end ).Note that this definition can be easily extended to directed events and to directed or undirected temporal hypergraphs.
Two events e, e ′ ∈ E are adjacent if they share at least one endpoint node in common, {u, v} ∩ {u ′ , v ′ } ̸ = ∅, and they follow each other in time such that ∆t(e, e ′ ) = t ′ start − t end > 0. Therefore, any temporal network can be represented as a higher-order static directed acyclic weighted graph known as the event graph D = (E, E D , ∆t(e, e ′ )) [20,25].Nodes of the event graph are the events of the original temporal network and the weight of a link between two connected nodes (adjacent events) is then defined as the time difference ∆t between the corresponding events.Every path on the event graph constitutes a causal chain as, by definition, a path constitutes a list of events where every two consecutive events are adjacent.Paths in event graphs are, therefore, equivalent to time-respecting paths in the corresponding temporal network representation [26].Therefore, calculating timerespecting reachability on a temporal network is equivalent to connectivity on its corresponding (static) event graph representation.The weakly connected components on an event graph determine causal domains, disjoint sets of events where there can be no causal connections whatsoever between events if they belong to two different weakly connected components.In addition, as compared to reachability, the size and distribution of weakly connected components are quantities, which are much easier to measure for temporal networks and they characterize a percolation transition if we assume an undirected network.Moreover, the sizes of these components put an upper bound on how much an effect can spread starting from one of the events in that component [20].
Temporal networks preserve the dynamic properties of the represented complex system, unlike aggregated static networks where this information is lost.Through the studies of time-varying networks, several new phenomena in human dynamics have been explored over the last decades, such as node and link burstiness [18,27,28], causal, temporal motifs [29], or the cyclic activation patterns of human interaction activities [30], to mention a few.As opposed to systems governed by homogeneous and independent processes, these correlations and the induced temporal dynamics may have significant effects on various dynamical processes evolving on temporal networks such as spreading [31,32], reachability [33,34], diffusion [9,35], and opinion formation [36].
The different dynamics of a temporal network are often straightforward to study through simulations.For example, in the case of spreading processes, transmission can be modeled by temporal network events [13][14][15][16][17].More concretely, in a physical interaction network, where nodes represent people and events represent two people coming to close proximity, each of these contact events will have a probability of transmitting the disease.The disease then spreads to all the nodes that can be reached via such infecting events from the initially infected nodes.Similarly, in a network where events represent communication of information at a specific time, such as mobile phone calls or email exchanges, it is straightforward to model the spreading of information by keeping track of the information nodes have access to at each point in time.
Many dynamics evolving on top of networks, such as some spreading processes [37][38][39], social contagion [40][41][42] ad-hoc message passing by mobile agents [43] or routing processes [44], can have a limited memory thus can only use paths constrained by limited waiting times.Limited waiting-time reachability can be modeled using the event graph, D, that contains a superposition of all temporal paths [20,26,45].In a limited waiting-time spreading process unfolding over a temporal network, ei-ther the spreading agent (e.g. the pathogen in the disease spreading) must be transmitted onward from a node within some time δt or the infection has to be renewed before that time.In other words, the node must participate in a possibly disease-carrying event in δt time, or the process stops and the node reverts to susceptible.Therefore, all the spreading paths in the network are δt-constrained time-respecting paths.Let's call two adjacent events e and e ′ as δt-adjacent if ∆t(e, e ′ ) ≤ δt.A subset of the event graph D with an upper threshold of weights no greater than δt, i.e., where directed links indicate δt-adjacency, enables us to calculate reachability for δt limited-time spreading process for the corresponding temporal network.Therefore, the event graph encapsulates a complete set of δt-constrained time-respecting paths for all values of δt simultaneously.

B. Directed Percolation
The waiting-time limit δt can be regarded as the control parameter of a continuous phase-transition, where connectivity in the event graph is determined by δtconnected paths of events.As the value of maximum waiting time decreases, more and more of the links of the event graph get removed, where each deleted link corresponding to an adjacency relationship between two events that are temporally more than δt apart.This leads to a drop in connectivity in the event graph, which is exactly equivalent to the drop in connectivity on the temporal network.In order to characterize these phase transitions, unlike characterizing the superficially similar phase transitions that take place when removing links in static (undirected) networks, we need to consider a percolation framework that can explicitly model the one-way flow of time.
Directed percolation is a paradigmatic example of dynamical phase transitions into absorbing states with a well-defined set of universal critical exponents and is often used to model phenomena with inherent directionality, such as fluids passing through porous media [23,[46][47][48].Originally introduced as a model for directed random connectivity [49], directed percolation attracted scrutiny in percolation theory in the late seventies [50].Since then, a considerable body of work has been devoted to this approach of interpretation in the literature since the critical behavior of many stochastic many-particle non-equilibrium processes can be shown to belong to the directed percolation universality class.Directed percolation has applications in various domains at multiple scales ranging from galaxies to semiconductors [51][52][53][54].
As the simplest model exhibiting a transition between active and absorbing phases [46], it is straightforward to define and implement models governed by directed percolation, e.g. in the case of lattice models [55][56][57][58][59][60][61].Directed percolation, however, does not appear to be an integrable model and its critical behavior is highly non-trivial.Moreover, it seems that the basic features of directed percolation, such as non-fluctuating states, are quite difficult to realize in nature [62].Another fundamental problem is quenched disorder due to microscopic inhomogeneities of the system [23].One of the earliest unambiguous and robust experimental realizations of a system exhibiting critical behavior in the directed percolation class was for the rather specific case of liquid crystal electrohydrodynamic convection [63].Another experimental evidence was reported in 2016 in the case of transition to turbulence [64].Due to the simplicity and robustness of directed percolation, it seems to be a good model for explaining ubiquitous phase-transitions in many real-world phenomena, especially in the so-called contact processes [4,33,[65][66][67][68][69][70][71][72] in the realm of temporal networks [24].
Before presenting the mapping between reachability in temporal networks and the concepts in directed percolation, for the remainder of this section we will review these concepts for the case of the simple infinite lattice.Let us take the example of a spreading process across time in an infinitely large d-dimensional square lattice: assume that each infected (or occupied) node can infect any of its neighbors independently with probability p at each tick of a discrete timer.Let us also assume that an infected node recovers (becomes unoccupied) in one tick of the clock after infection unless it is re-infected by a neighbor.This configuration is denoted in many sources as a d + 1-dimensional lattice, substituting the temporal axis with another discrete spatial dimension with the only difference that, unlike the other d dimensions, this one has an inherent directionality.Throughout the rest of this section, we will continue to use the space and time analogy to facilitate a better transition to modeling phenomena on temporal networks.
The dynamics of this spreading process is defined by the topology and dimensionality of the medium of percolation and competition between two processes: the probability that an infected node infects each of its neighbors in a single tick of the clock, or "reproduction" from the perspective of the spreading agent, and the time it takes for each infected node to recover, or "self-annihilation" or "death" of the spreading agent.In the many classic representations of directed percolation, the reproduction probability is often denoted by the parameter p and the "self-annihilation" is set to happen in exactly one tick of the clock.For large enough values of p, the system will forever stay in an "active state" where there is a non-vanishing density of nodes infected (occupied) at all times.Conversely, if the annihilation process has the upper hand, the system eventually transitions irreversibly into an "absorbing phase" where no occupied nodes are left in the lattice and the spreading agent is extinct.
More generally, let us say the reproduction and selfannihilation process respectively happen at rates µ p and µ r .Let us assume that at t = 0, nodes are uniformly occupied with density ρ 0 .To write a mean-field rate equation for occupation density ρ(t), we need to take into account how often more than one spreading agent (pathogens) simultaneously occupies (infects) the same node, in which case only one new node is occupied.Let us only consider the rate µ c at which two other nodes simultaneously infect a single node and assume the probabilities of three or more simultaneous infections are small.In this case, the rate equation is of the form where the control parameter τ = µ p − µ r is the manifestation of the competition between reproduction and death as described above and coupling constant g = µ c describes the events of infecting a node already infected by another neighbour [23].This equation has a steadystate at lim t→∞ ρ(t) = ρ stat (τ ) = 0 which corresponds to the aforementioned absorbing phase.Furthermore for τ > 0 the value of ρ(t) approaches a stationary occupation density of lim t→∞ ρ(t) = ρ stat (τ ) = τ /g, which is identified as the order parameter of the directed percolation process.At exactly τ = 0, occupation density decays algebraically with time ρ(t) ∼ (ρ −1 0 + gt) −1 .Naturally, for values of τ < 0 the system eventually arrives at the absorbing phase ρ(t) → 0 in finite time.
More generally, starting from a homogeneously occupied initial condition, order parameter ρ stat (τ ) of a system in the directed percolation universality scales as ρ stat (τ ) ∼ τ β , when control parameter τ is close to τ c = 0.For τ > 0, density decays algebraically as ρ(t) ∼ t −α where in the mean-field regime (i.e.d ≥ 4), β = α = 1.In the case of a spreading process controlled by a percolation probability p introduced at the beginning of this section, it can be shown that τ ∝ p − p c where critical percolation probability p c is a function of topology and dimensionality of the percolation medium [73,74].
Alternatively, we can focus on the ramifications of starting from a single seed of infection, as opposed to a homogeneous initial distribution of occupied nodes.A characteristic property of this scenario is survival probability P (t): the probability that a spreading process starting from a single seed would still be in the active phase (ρ(t) > 0) at time t.Similar to occupation density ρ(t), at criticality τ = τ c = 0 survival probability also decays algebraically with time P (t) ∼ t −δ .A second alternative for order parameter is the ultimate probability of survival P surv (τ ) = lim t→∞ P (t).When the control parameter is close to the critical threshold τ → 0 − , the ultimate probability of survival scales algebraically as P surv (τ ) ∼ τ β ′ .
Continuous phase transitions in models with time-like dimensions generally have the same system of two separate order parameters, controlled by two different critical exponents β and β ′ .For the case of directed percolation, however, "rapidity-reversal symmetry", an invariance property under time-reversal, ensures the two exponents have the same value β = β ′ [75] which implies that P (t) and ρ(t) are at least asymptotically proportional as t → ∞, and in some cases exactly equal P (t) = ρ(t) [23].Rapidity-reversal symmetry limits the number of independent critical exponents to three [76,77].

Characteristic quantities of the directed percolation
The single-source initial condition also allows us to define additional interesting characteristic quantities in the absorbing phase, which might lend themselves to experimental observation.Let's define pair-connectedness function c( ⃗ r 1 , t 1 , ⃗ r 2 , t 2 ) as the probability that a path exists from a node with spatial coordinates ⃗ r 1 at time t 1 and another in ⃗ r 2 at time t 2 .Note that the definition of spatial coordinates for nodes as a d-dimensional vector ⃗ r i implies that the percolation medium and node i is embedded in a d-dimensional space, e.g. a d-dimensional lattice.Assuming that the percolation medium is invariant with respect to translations across time and space, we can simplify the pair-connectedness function by fixing the origin on the source node and denote the pair-connectedness function as c(⃗ r, t).Mean cluster mass M is defined as the integration of the pair-connectedness function across time and space: which, with control parameter close to the critical threshold Similarly, mean spatial volume V can be defined as the number of unique nodes that will ever get infected in a single-source spreading scenario.As with the case of the cluster mass M , spatial volume scales through a power relationship V ∼ (−τ ) −υ close to the critical threshold where υ = dν ⊥ − β ′ .It is possible to think of spatial volume V as the size of the projection of the percolation cluster over the d-dimensional spatial plane, i.e., over the original d-dimensional lattice.Projection of the same cluster on the temporal dimension will define the survival time of the cluster, which is distributed according to the probability of survival P (t).The homogeneous, fully-occupied initial condition, on the other hand, allows us to study the response of a system to an external field h on the order parameter static density ρ stat .For the case of directed percolation, an external field can be implemented as the spontaneous occupation of nodes at a rate h.A positive external field deprives the system of the possibility of ever transitioning into an absorbing phase.Susceptibility χ is defined as the magnitude of the response generated by a minuscule disturbance in the external field which diverges algebraically as the control parameter τ converges to the critical threshold where γ is the same exponent as the mean cluster mass M .For the rest of this paper, when not specified, susceptibility χ is studied at minuscule values of external field (h = 0) as τ converges to the critical threshold τ c = 0.In practical terms, susceptibility is a useful tool for finding the transition point, as unlike the order parameters, we do not need to define an arbitrary threshold for what constitutes a small or large value for a quantity such as M (t) or V (t) close to the transition point in a finite system.Instead, the susceptibility will typically show a peak even in finite systems, which are discussed in more detail in Sec.I B 2.

Finite size scaling properties of the system
While the dynamics described previously explain the behavior of an infinitely large system, measuring properties of infinitely large systems is a rather involved task.Verifying that the behavior of a system at criticality is explained by a specific set of critical exponents is often easier performed by studying the finite-size scaling properties of the system.This can be carried out by measuring a set of quantities for realizations at different scales and plotting the universal scaling function of each quantity as a function of scale-invariant ratios.If the exponents used are correct, all the scaling functions of different linear system sizes for the same quantity should collapse on top of each other.
The effect of the finite size of the system manifest themselves as deviations from the scaling laws as described before and their effects are measurable after some characteristic size-dependent amount of time has elapsed since the beginning of the simulation.For example, while in an infinitely large system in active phase τ > 0 the system will forever stay in an active phase, a finite system will always have a non-vanishing probability of transitioning to the absorbing state due to fluctuation of the order parameter.These finite-size effects take place at a characteristic time t f that scales as t f ∼ l z where z = ν ∥ /ν ⊥ is the so-called dynamical exponent and l is the lateral (or linear) size of the system as opposed to system size N measured in number of nodes, where N ∝ l d .
In phenomenological scaling theory simple scaling is assumed for absorbing phase transitions.This means that large-scale properties of the system are invariant under scale transformations with the control parameter close to the critical threshold.A multiplicative transformation, or "concentration", of the control parameter τ by a factor of λ, τ → λτ would result in re-scaling of other quantities as where t and l stand for time-like and length-like quantities respectively.More specifically, scale-invariance mandates very specifically how a quantity will change under multiplica-tive scale change.As an example, let's study changes of where t is time from initial infection seed and l is the linear system size.Since this relationship is valid for all values of λ, we can remove one parameter of the function by selecting a special value where the function F (x) is referred to as the "(universal) scaling function" of its corresponding quantity, in this case, density ρ.The parameter to this function l −ν ∥ /ν ⊥ t is in itself invariant to scale transformations.This parameter and those similarly derived for other quantities are oftentimes known as "scale-invariant ratios".The function F (x) is universal, meaning that if measured to sufficient accuracy, we obtain exactly the same type of scaling function for systems with similar boundary conditions and shape for any phenomena in the directed percolation universality class [23].
The value of each exponent is only a function of a few large-scale properties of the system, such as the number of spatial dimensions of the system.There exists an upper critical dimension d c where systems with spacial dimensionality d ≥ d c all follow the same set of values for critical exponents, which are exactly equal to those derived through mean-field estimation.For the case of the directed percolation universality class the upper critical dimension has a value of d c = 4 [23].

A. Directed percolation in temporal networks
Let us now take the case of δt limited-time spreading from a single source on a temporal network.Similar to the classic directed percolation single-source spreading process, each temporal network node can participate in the spreading process by becoming infected, infecting others and recovering multiple times.Temporal networks are different from the archetypal directed percolation systems presented in Sec.I B in that they do not present a regular lattice or metric space in the spatial dimension.Furthermore, there is typically no discrete structure in the temporal dimension, which is usually modeled as a continuous axis.Nevertheless, if the various concepts such as order parameter, control parameter and cluster sizes are defined carefully, temporal networks and limited waiting-time connectivity can be mapped to directed percolation [22].
To put it in the same reference frame as with other absorbing phase transitions, changing the parameter δt, in this scenario, controls the relative occurrence of "annihilation" and "multiplication" processes.A small enough value of δt, will lead to a situation where spreading scenarios will eventually die out, at which point the system enters an absorbing phase.Similarly, as δt grows, a spreading agent will be able to avoid extinction for longer time, until after some threshold δt > δt c a random spreading scenario will not die out (in an infinitely large network).As discussed in Section I A, such spreading scenarios are closely related to various properties of the δt thresholded limited waiting-time event graph D.
As illustrated in Fig. 1, the projection of the spreading cluster over the spatial plane amounts to a subset of temporal network nodes V that has ever participated in the spreading process.This can be measured by calculating the mean number of unique temporal network nodes involved in the out-components of the event graph.The (ensemble) average number of unique nodes participating in single random source spreading processes is analogous to mean spatial volume V .The projection of the spreading cluster over the temporal axis is equal to the time window from the beginning of the spreading process to its end.The ensemble average of this time duration is analogous to mean survival time T .The sum of the duration of infectiousness for all the nodes, i.e. the integration of the pair-connectedness function, would therefore be analogous to spatial and temporal integration of the pair-connectedness function or mean component mass M .Note that the duration of the infectiousness is equal in all of the events, therefore, we use the number of reachable (i.e.possibly infection-carrying) events as a proxy for M , ignoring the overlaps.The above-defined quantities can be measured as features of the event graph.The average number of events in the out-component of a node in the event graph (equivalent to an event in the temporal network) measures the number of reachable events.The survival probability P (t) can be similarly defined over an ensemble of single-source spreading instances based on the distribution of the lifetime of each spreading scenario, accounting for the finite temporal window of the simulation of the temporal network using a Kaplan-Meier estimator [78].
Another scenario is the simulation of the spreading process from homogeneous, fully occupied, initial conditions.Translating this from classic directed percolation poses a new problem; a homogeneous initial condition cannot translate to a "full row" of occupied nodes since we are dealing with continuous-time as opposed to the typical directed percolation case of discretized time presented in Sec.I B. Rather, a better translation of the fully-occupied initial condition to continuous time is to assume all nodes to be occupied at the beginning of the observation period t = min(T ), or more accurately by assuming all nodes to be occupied for all values of t where t ≤ min(T ).Occupation density ρ(t) is defined as the fraction of infected nodes at time t.Stationary density ρ stat (τ ) is therefore defined as occupation density after the system had enough time to reach a stationary state.We can also emulate the effects of an external field h in this scenario: In continuous-time, this is equivalent to each node spontaneously becoming occupied through an independent Poisson point process with a rate of h.Susceptibility χ(τ, h) can then be measured, from Eq. ( 3), by the rate of change in stationary density as external field changes.

B. Empirical methods for estimation of characteristic quantities
In practice, we can estimate M , V , T and P (t) on the event graph by finding all the out-components, i.e., every reachable event starting from every event (see, Figure 1(c)).Calculating the exact set of out-components for every event in the event graph is time and memoryintensive.However, if we are only interested number of events or number of unique nodes that participate in those events, as opposed to the full set of events in the out-components, we can use probabilistic cardinality estimation data-structures to estimate out-component sizes with arbitrary precision in O(|E| log |E|) time, as opposed to O(|E| 2 ) time required for exact calculation [45].Minimum and maximum time of all events in the outcomponent can be exactly calculated in O(|E| log |E|) time.Calculating properties of the in-component of an event is possible through a simple reversal of direction of all links in the event graph and applying the same algorithms.
Similarly, in the homogeneous fully-occupied initial condition scenario, we do not need to directly estimate occupation density ρ(t), stationary occupation density ρ stat (τ ) and susceptibility χ(τ, h) via naive algorithms, which would explicitly compute these measures by simulating propagation.The properties of homogeneous, fully occupied, δt-constrained reachability can be estimated by marking as occupied any event that is in the outcomponent of at least one event with time −∞ < t < t 0 .This can be accomplished by running the in-component size estimation algorithm [45] once over the whole network, recording minimum observed time in in-component of each event and marking those with minimum incomponent time smaller than t 0 as occupied.In practice, temporal networks are only recorded or generated for a finite window of time t min < t < t max .As there are no adjacency relationship between events more than δt apart temporally, any event that has at least one event in its in-component with time t min < t < t min + δt can be considered occupied.Figure 1(d) shows all occupied events (dark grey) with the initial condition that assumes all nodes are occupied from −∞ < t < 0. The density of occupied events, which corresponds to particle density ρ(t), can be estimated from the event graph representation by the number of occupied nodes in a band of time divided by the area covered by the band, i.e. number of nodes multiplied by the width of the band.
Normally, calculating the effects of an external field h would require simulating a fully occupied initial condition, marking some nodes randomly selected with rate h V can be visualized as the mean size of projection of a spreading cluster on the spatial plain, i.e. the static base network on the vertical axis, whereas survival time T is equivalent to the mean size of projection of the cluster on time (horizontal) axis.While measuring a direct analogue to component mass M , integrating the pair-connectedness function across time and space equivalent to the mean sum of lengths of colored horizontal lines in (a) is not straightforward with the event graph representation.It is possible to show that the mean number of uniquely counted events involved in the spreading process, corresponding to the cardinality of the out-component of the initial event here represented by the total number of colored nodes in the event graph, show the same scaling behavior.(d) Homogeneous, fully-occupied initial condition with the occupied events shown in a darker shade than unoccupied events shows the decline and eventual stabilization of the occupation density as time grows.In this scenario, all nodes are considered occupied for time −∞ < t < 0, which translates to the occupation of all events in period 0 ≤ t < δt and all events in their out-components.
as occupied, computing their out-components, and measuring how many new events got occupied.As we are interested in the effects of a minuscule positive external field, indicated by susceptibility χ(τ, 0), we can instead calculate the effects of spontaneously marking exactly one random event in the whole network as occupied using probabilistic counting and in-components of all events (i.e., looking back in time).If the number of events in the in-component of an event e is denoted as |E in (e)| and the minimum time among all events in its in component as t in min (e) = min (u,v,t)∈E in (e) t, the expected number of spontaneously occupied events when a minuscule external field h is applied can be estimated as e∈E P occupied (e) where In this scenario, the respective value for the external field that would spontaneously occupy on average one event is proportional to h ∝ 1/|E|.We approximate ρ(t) by number of occupied events within a δt time window di-vided by spatio-temporal hyper-volume of the time window δt × |V|.The estimate for ρ(t) can in turn be used to approximate quantities like stationary density ρ stat (τ ) and susceptibility χ(τ, h).

A. Experimental setup
In this section, we focus on validating and exploring the limits to our hypothesis that δt limited-time spreading in many forms of temporal networks belongs to the directed percolation universality class.We do this by performing single-seed and homogeneous initial-condition spreading simulations following the method defined in Sec.II A and explained in detail in Sec.II B. By measuring various observables for networks of different sizes as described in Sec.I B 2, we can verify whether for each quantity the corresponding universal scaling functions collapse for systems of different finite sizes when using the same values of critical exponents β, β ′ , ν ∥ and ν ⊥ as that of DP corresponding to the dimensions of the system as a previous mean-field approximation and experimental setups for the directed percolation.
The experiments are performed on a variety of synthetic temporal networks.The generation procedure consists of generating a static base network corresponding to the aggregate network and generating events, i.e. activations or timestamps, for each link based on some temporal dynamic.In total, we analyzed 26 combinations of base networks and link-activation processes.In order to perform the finite-size scaling analysis, we computed all the statistics for ten network sizes, starting from N = 2 8 nodes and increasing the size by a factor of two until we reached N = 2 17 nodes.For the case of d-dimensional square grids where d ∈ {2, 3, 4}, however, closest powers of d to the powers of two from 2 8 to 2 17 was used with a periodic boundary condition, to provide spatial translational invariance.Each statistic was calculated as the average of at least 256 (up to 4096) realizations and each realization of the largest configuration consists of around 3.7 × 10 7 events.No sampling of spreading scenarios was required for each network's realization, as the effect of starting a spreading process from any possible combinations of nodes and times could be gathered in one pass as described in Sec.II B. See supplementary material for a more detailed overview of the experimental setup.
Static base networks are either (a) one to fourdimensional square lattice grids with periodic boundary conditions (b) random regular graphs with specified average degree [79,80] or (c) Erdős-Rényi G(n, p) random networks with specified expected average degree [81].For the random networks, we chose the average degrees 8 for the Erdős-Rényi graphs and 9 for the random regular graphs (such that both networks have the same expected excess degree).The higher degrees of random networks ensure that the probability of generating networks with large isolated components remains negligible and that even locally, the network would be of high enough dimensionality to be in the mean-field regime above the upper critical dimension d c = 4.
Temporal dynamics of the links are either governed by (a) Poisson processes, i.e. exponential inter-event times, (b) bursty processes, i.e. renewal processes with powerlaw inter-event time distributed as ∝ ∆t − γ with exponents γ ∈ {2.05, 2.2, 2.8, 5.2} and minimum interval cutoff set so that the expected inter-event would be equal to one and (c) Hawkes independent self-exciting processes with different parameter sets.The Hawkes univariate exponential self-exciting process [82], is defined by the conditional intensity function The parameters of this formulation of the Hawkes process are (1) background (or exogenous) intensity of events µ indicating the random probability of events happening without being caused through self-excitement, (2) the infectivity factor α, which can be interpreted as the ex- pected number of induced self-exciting events per each event, and (3) the rate parameter of the delay θ.Based on the properties of exponential kernel used in defining Eq. 8, 1/θ is the expected inter-event time between an event (e.g. a coincidental social interaction) and its corresponding induced self-exciting event (e.g. the follow-up social interactions) [83].
As the unit of time is arbitrary, temporal processes are scaled, without loss of generalization, so that they produce timestamps with a mean inter-event time equal to one.The processes are initialized in their stationary state, and in practice, the first timestamp for each event is generated through residual time distribution of each process, except for the case of Hawkes process where the process is allowed a burn-in time equal to the simulation time window before the first timestamp is recorded.The temporal processes of pairs of links are simulated independently of each other.Figure 2 shows a visualization of the different methods of generating event activations.Temporal networks were simulated for a time window of at least T = 64 and up to T = 8192 units of time.See supplementary material for the exact experimental setup for each system size.The difference in system sizes and time windows for the simulations were necessitated by the limitations and optimal utilization of the computational facilities.
B. Estimating the critical threshold δtc and the critical exponents β, β ′ , ν ∥ and ν ⊥ Best estimate of the critical exponents β, β ′ , ν ∥ , ν ⊥ and critical threshold δt c can be determined by finding the values of these exponents that would produce the best data collapse for the universal scaling functions corresponding to ρ(t), P (t), M (t), and V (t).The quality of collapse, in turn, can be assessed by comparing the deviation of the scaling function curves for different system sizes from the average trajectory.Here, for each of the quantities P (t), ρ(t), M (t) and V (t), we calculated one trajectory for finite-size scaling function for each system size, as defined for example for the case of ρ(t) by Eq. ( 6).As the tested value of critical exponents and δt c gets closer to the actual critical threshold, the curves for different sizes should more closely collapse on top of each other.Plotted with the correct values of critical exponents and critical threshold, we expect to see all trajectories collapse into one with the possible exception of very small values of t.To quantify the quality of a collapse, we measure the mean curve in the area where all system sizes have defined values for the scaling function and measure the root mean square difference of all points from all system sizes to the mean curve.The errors were measured after logarithmically scaling the values to account for the power-law nature of the scaling functions.Sum of errors for the collapse of P (t), ρ(t), M (t) and V (t) was used in evaluating each set of parameters.
In order to assess collapse of the universal scaling functions, we first determine a value for δt c for each network configuration.That is, the best candidate for δt c is selected based on the least total error for collapse of P (t), ρ(t), M (t) and V (t) assuming DP critical exponents.Figure 3 shows this total error of collapse for two network configurations.This shows is a clear minimum for each configuration indicating the critical value δt c , which is consistent across P (t), ρ(t), M (t), and V (t) trajectories.
The resulting estimates for δt c can be used to visually verify directed percolation critical exponents and our selected optimal value of δt c for each system by plotting the finite-size universal scaling functions of different system sizes.In total, we produce collapses for eight characteristic quantities measured in a single source or homogeneous initial conditions.Figure 5 shows these collapses measured for regular networks with bursty dynamics (renewal process with power-law inter-event times) and Erdős-Rényi networks with a Hawkes self-exciting process dynamics.The full set of plots for all of the 26 configurations are shown in Supplementary Materials.In all cases, a satisfactory collapse can be observed for at least probability of survival P (t) and density ρ(t) and in most cases, other quantities show a good collapse as well.It is important to note that quantities that de- FIG. 3. Root mean square (logarithmic) deviation of scalingcorrected functions of probability of survival P (t), density ρ(t), mass M (t) and volume V (t) for different system sizes from average trajectory shows a sharp drop at δtc due to data collapse.Each instance of the network is made through realizations of (a) Erdős-Rényi static network ⟨k⟩ = 8 and Poisson process λ = 1 activations, (b) random 9-regular networks with bursty (power-law with minimum cutoff) interevent time distribution with mean 1 and exponent γ = 2.8, (c) Erdős-Rényi static network ⟨k⟩ = 8 and Hawkes univariate exponential self-exciting process with parameters µ = 0.2, α = 0.8 and θ = 1.0 and (d) one-dimensional grid with periodic boundary conditions (a circle) and Poisson process λ = 1 link activations.Refer to Sec.III A for the definitions of the parameters.
pend on measuring values as time approaches infinity, e.g., ρ stat (δt) and χ(δt) have generally lower quality of measurement and collapse since the time to reach a stable value for these increases substantially close to criticality [46].
Table I (column "Est.δt c ") shows our best estimate of the critical threshold δt c for each configuration using the method described above.As the systems become rapidly more and more connected after the critical threshold, a lower value for the critical threshold δt c indicates higher, or more robust, spatio-temporal connectivity, meaning that the same δt limited-time spreading would result in a larger number of reachable nodes, V (δt), or larger number of reachable events, M (δt).When modeling infectious disease spreading as directed percolation on temporal networks, larger values for V (δt) and M (δt) may indicate larger epidemic sizes and the total number of human hours of infection in the population, respectively.
These results indicate that within each spatial configuration, increased burstiness (as indicated by lower value for the power-law exponent γ) generally leads to a lower TABLE I. Column "Est.δtc" shows the best candidate for critical threshold δtc selected by minimizing the collapse error of the universal scaling functions for probability of survival P (t), density ρ(t), mass M (t) and volume V (t) derived for different system sizes, assuming DP exponents.The collapse error is measured by the sum of root mean squared deviation of logarithmically scaled trajectories for all four scaling functions.For the best estimate for exponents α = β/ν ∥ and δ = β ′ /ν ∥ , reported respectively in columns "Est.α" and "Est.δ", a power-law was fitted to the head of values of probability of survival P (t) and density ρ(t) respectively for the largest system size simulated for the time period between 2 < t < 0.04 × N ν ∥ /dν ⊥ , where both functions are expected to be still mostly behaving, similar to an infinite system, according to power relations t −α and t −δ respectively.Directed percolation mean-field values for these critical exponents are δ = α = 1, which is close to the value estimated for random, high-dimensional networks.Furthermore the value of these critical exponents in a DP system are expected to be close to α = δ = 0.15946 for 1+1 dimensional, α = δ = 0.450 for 2+1 dimensional, α = δ = 0.732 for 3+1 dimensional and equal to the mean-field estimates systems α = δ = 1 for 4+1 dimensional [46], which is close to values estimated for 1 dimensional lattice and 2 to 4 dimensional square lattices.

Configuration
Est value for δt c threshold and higher connectivity.Furthermore, for the case of the self-exciting process, increasing the expected number of self-induced events, as indicated by α, generally results in a lower value for δt c (higher connectivity).While it was previously understood that a wide range of temporal inhomogeneities slows down spreading processes over temporal networks [13], these results demonstrate that certain temporal inhomogeneities, e.g. a highly bursty or self-exciting temporal dynamic, can enable a more limited spreading agent (expressed in terms of a maximum waiting time) to spread to a wider set of nodes.For example, spreading processes with maximum waiting time between 0.063 < δt < 0.084 over an Erdős-Rényi networks ⟨k⟩ = 8 will spread to a much larger set of nodes and span a longer span of time if the link activations are highly bursty (γ = 2.05) compared to a Poisson process with the same mean inter-event time, as the latter will be spreading in the sub-critical regime compared to the super-critical regime for the former.
It is also interesting to note that while the random spatial configurations, namely random 9-regular networks and the Erdős-Rényi networks ⟨k⟩ = 8, both result in networks with the same expected excess degree value, the Erdős-Rényi networks with higher levels of spatial inhomogeneity, which manifests as a wider spread degree distribution, can be observed to have a lower δt c critical threshold.While testing on a wider range of spatial (structural) inhomogeneities would be required before a conclusion is reached, these results might hint at a similar behavior as with temporal inhomogeneities, namely that introducing certain spatial inhomogeneities might result in higher connectivity in the sense that the same limited-time spreading agent can eventually spread to a wider share of the network.
Additionally, we present a method to assess the quality of a collapse for a range of different values of critical exponents (β, β ′ , ν ∥ and ν ⊥ ) and δt c .A five-dimensional grid search for optimal values for critical exponents and δt c based on the quality of collapse for P (t), ρ(t), M (t) and V (t) shows that the total error declines around critical exponent values close to that of directed percolation, i.e. β = β ′ = ν ∥ = 1 and ν ⊥ = 0.5 for mean-field regimes and their respective DP values for lower dimensionality square grid networks.Figure 4 shows for Erdős-Rényi static networks with ⟨k⟩ = 8 and Poisson process link activation, the β × β ′ plane from the five-dimensional grid search with two sandwiching parallel planes along each of the ν ∥ , ν ⊥ and δt c dimensions.This verifies that there is a minimum close to β = β ′ = ν ∥ = 1, ν ⊥ = 0.5 and δt c = 0.08421 for total error of collapse of P (t), ρ(t), M (t) and V (t).Similar plots for some other network configurations (along with a different two-dimensional slice, ν ∥ × ν ⊥ ) can be viewed in the Supplementary Material.It is important to note that while other combinations of parameters in the grid might lead to other local optima, visual inspection of the resulting collapse show that to be mainly numerical artifacts where the total error changes rapidly close to extreme values of the parameters (i.e., critical exponents and δt c ) where only a very small fraction of the trajectories for different finite sizes actually overlap.
Furthermore, for each of the critical exponents, we can measure an estimation error based on this fivedimensional parameter grid.For each exponent, we find a range of values where, assuming that all other exponents are fixed at their DP values, would produce collapses of higher or equal quality compared to the DP value of that exponent.The sizes of these ranges, which by definition includes the DP value for all exponents, provides a confidence interval for the range of possible exponent values that are able to explain the behavior of the system with at least the same quality as that of directed percolation.As shown in Table I, these errors are in most cases only a few percent, with a notable exception of the highly bursty renewal processes with γ = 2.05.Simulating power-law distributions becomes a much harder problem as the magnitude of the exponent approaches 2. Close to this exponent, it takes a larger and larger number of realizations for the properties of the population, e.g.average inter-event time for bursty temporal dynamics, to converge.It is also possible that the large estimation error is an indicator that the system is approaching a breakdown of one of the key symmetries, with the most likely candidate being rapidity-reversal symmetry based on the fact that the estimation error for β ′ is much larger than that of the other exponents.

C. Estimating critical exponents by simulating very large systems
As discussed before in Sec.I B, the effects of the finite size of the system manifest at characteristic times t f ∝ N ν ∥ /dν ⊥ in the form of fluctuations that causes the transition of the system to the absorbing phase.At times much smaller than t f the system shows approximately the scaling behavior of an infinitely large system where at criticality, ρ(t) ∼ t −α and P (t) ∼ t −δ where α = β/ν ∥ and δ = β ′ /ν ∥ .On the other hand, the power-law scaling behavior becomes visible at times comparable to the mean inter-event time of the dynamic process but not up to arbitrarily infinitesimal values of t.Given these properties, we fitted two power-law functions using the least-squares method to the results of experiment with the largest system size for the range of time 2 < t < 0.04 × N ν ∥ /dν ⊥ on ρ(t) and P (t) to derive exponent α and δ. Figure 6 shows one such fitting for a system made from Erdős-Rényi networks with ⟨k⟩ = 8 and N = 2 17 nodes and bursty (power-law with minimum cutoff) inter-event time distribution with mean 1 and exponent γ = 2.8.Table I (columns Est.α and Est.δ) shows the best estimates of these exponents, which as expected are very close to respective directed percolation critical exponents of 1 (for the mean-field regime d ≥ 4) for the case of random networks and 0.159, 0.450 and 0.732 for one-, two-and three dimensional lattices respectively [46]. FIG.
6.An example of fitting power-law functions on empirical ρ(t) and P (t) results on finite networks for deriving critical exponents α = β/ν ∥ and δ = β ′ /ν ∥ .Power-law functions were fitted on experimental results of spreading over Erdős-Rényi networks with ⟨k⟩ = 8 and N = 2 17 nodes and bursty (powerlaw with minimum cutoff) inter-event time distribution with mean 1 and exponent γ = 2.8.The fitting was performed on values in range 1 < t < 0.04 × N ν ∥ /dν ⊥ (i.e. 1 < t < 14.48) to limit the interference of finite-size effects with the scaling behavior.

IV. DISCUSSIONS
Through combining multiple methods of empirical and theoretical verification, we are able to confidently state that limited waiting-time connectivity percolation over a wide range of synthetic temporal networks incorporating a range of temporal and topological inhomogeneities show behavior compatible with the directed percolation universality class.It is of utmost importance to discuss the limitations of our method: chief among them, that our empirical finite-size simulation method, as described in Sec.II B, is not able to measure quantities which are defined at t → ∞, such as the ultimate probability of survival P surv and static density ρ stat (and therefore susceptibility χ) to the same standard of accuracy as the other quantities due to the finite size of the synthetic networks used for analysis.This is exacerbated close to the critical threshold where the equilibration time, the time required for the network to reach a stationary state, grows rapidly while the memory and computational cost of simulating a temporally larger temporal network grow linearly and log-linearly, respectively, with the increased simulated time [46].This is visible in Fig. 5e,f,k,l as a worse collapse as compared to other quantities.Also, while it is computationally much more feasible to measure susceptibility χ by inducing occupation of exactly one existing event in the temporal network (described in Sec.II B) as compared to inducing occupation of nodes at random times (as described in Sec.I B), the latter method might be more robust, especially when dealing with a temporal network with a high degree of temporal inhomogeneity.Although our experiments with this alternative method were limited to smaller system sizes, we could not observe any significant difference between the two methods for the network configurations presented in this manuscript.While a wide range of temporal dynamics and network structures with different levels of inhomogeneity are studied here, there is still a wide variety of systems that present computational and theoretical challenges.First, the effects of event-event correlations between links are not studied.It has been shown that event-event correlations, among other forms of inhomogeneity, can affect the rapidity of the spreading process on temporal networks [13].Conceptually, local event-event correlations such as temporal motifs [29], are close to temporal event graphs, which are in practice computed using isomorphisms on slightly modified temporal event graphs.Thus, incorporating temporal motifs to the framework at the level of analytical computations is an interesting future direction, as that corresponds to modifying the frequency of appearance of structural motifs in the event graphs.Second, the effect of static base networks with heavy-tail degree distributions and other more complicated network topologies are absent from this study.Here, of especial interest are the networks with heavy-tailed degree distribution with static network reachability percolation threshold at zero occupied links, e.g.p(k) ∝ k −2 .While initial results did not support the conclusion that a δt limited waiting time over this class of synthetic temporal networks would be in the directed percolation universality class, due to limitations on computational resources, we were not able to perform the analysis on the larger system sizes comparable to the other types of networks.
Depending on the physical mechanism involved in the modeled connectivity phenomenon or spreading process, alternative methods of defining the adjacency relationship might be more suitable than the one used here.For example, for the case of disease spreading over a physical contact network, the currently used definition of event graph causes a "re-infection" of the infected party, manifested as a restart of their δt duration of disease.This can be resolved by substituting each undirected event in the temporal network with two simultaneous directed events.Similarly, for a disease spreading scenario over transportation networks, such as an airplane traffic network, the time between two events (the value that is compared to the maximum duration of disease δt to determine whether two flights are adjacent) should be calculated from the departure of one flight to departure of the possibly adjacent flight and not, as it is currently presented, from the arrival of the latter to the departure of the former.This might be an important factor when dealing with scenarios in which the reasonable values for δt are comparable to the delay or the duration of the events, e.g. the time from the departure of a flight to the arrival in a spreading process over an air transport network.
For some spreading mechanisms, it might also be more suitable to replace the hard δt limited-time cutoff of adjacency requirement used in this work with a probabilistic process by measuring quantities over an ensemble of event graphs.For example, using a Poisson process instead of a δt limited-time cutoff would produce dynamics similar to simulations of SIS processes over networks while simulating results of the simulation starting at every possible starting point in one pass.Viewed this way, normal δt limited-time cutoff can be seen as a probabilistic process where the probability of adjacency is a step function at ∆t = δt.It is also possible to combine an occupation probability similar to classic directed percolation (see Sec.I B) with a δt limited-time cutoff (or a Poisson process cutoff or other forms of temporal locality constraint) to construct a two-dimensional phase diagram for each temporal network.
It would also be possible to define connectivity in the event graphs in a way that mimics the SIR process.In this case, one would need to prune some of the temporal paths in the event graph such that temporal network nodes are not repeated.This distinction is equivalent to paths and simple paths (or walks and paths, respectively) in static graphs.The algorithmic techniques employed in this work are not directly applicable to this case, and in fact, it has been recently shown that algorithmic problems in such settings can be computationally difficult.For example, in the SIR interpretation of the event graph, finding if it is possible for a node infected at a specific time to infect a given node is an NP-hard problem [84].In any case, averaging over explicit simulations of spreading scenarios is always an alternative option to the algorithms that take advantage of the redundancies in computing reachability.
Connectivity, which encapsulates several important phenomena on complex systems such as spreading pro-cesses [37][38][39] and routing dynamics [44], has not yet undergone the same level of development on temporal networks as the static networks.It has been previously suggested that connectivity on temporal networks, or other adjacent representations such as dynamic networks, might show the same properties as any other directed percolation system [20,21], a class of percolation models with built-in directionality which has enjoyed abundant attention in the past decades.In Ref. [22], we laid formal foundations by providing one-to-one analogues between concepts from directed percolation and temporal network connectivity and provided theoretical evidence supporting this hypothesis.In this work, we presented multiple accounts of empirical evidence showing that connectivity on many model temporal networks belongs to the directed percolation universality class and that this hypothesis is robust for a range of temporal and spatial heterogeneities.
This work focused mainly on establishing the vocabulary and developing the required tools in the hopes of rendering studies of connectivity in temporal networks ripe for future analysis, especially from a critical phenomena perspective.It is important to note that this work has only scratched the surface of the analytical study of connectivity on temporal networks and still, a vast body of analytical and phenomenological topics, some of which were eluded to in the previous paragraphs, remains open for future study.
The event graph representation contains many redundant adjacency relationships, e.g.triangles, or more generally feed-forward loops, that can be removed without changing reachability of nodes, producing a reduced event graph [22].Assuming the probability of two or more adjacent events happening at exactly the same time is negligible, the reduced event graph, a subset of the event graph with exactly the same reachability properties, has a maximum in-and out-degree of 2 [22,25].If we make the simplifying assumption that the reduced event graph representation of δt limited waiting-time spreading process on a specific temporal network is indistinguishable from a random directed network with the same joint inand out-degree distribution P (k in , k out ), a mean-field solution to order parameter occupation density ρ stat for a δt limited-time spreading process over temporal networks, as defined in Section II A, can be derived in the form where ⟨Q out ⟩ is the mean excess out-degree of the reduced event graph [22].This rate equation has the same form as Eq. ( 1).The solution to this equation shows a phase transition at τ c = 0 and other behavior consistent with τ = ⟨Q out ⟩ − 1 being the control parameter of directed percolation.As with Eq. ( 1) this sets two of the four critical exponents in the mean-field regime to the same values as those of mean-field DP, α = β = 1.Under the same assumption, the probabilitygenerating function representation of the out-degree distribution is G out 0 (y) = G(1, y) where G(x, y) is the joint in-and out-degree distribution probabilitygenerating function.Similarly, the excess out-degree distribution probability-generating function can be defined as where is the mean in-or out-degree on the event graph.This can be used to derive the out excess-degree distribution as Making the same assumption as above, namely that the event graph representation is indistinguishable from a random directed network with the same joint in-and out-degree distribution P (k in , k out ), we can derive the mean cluster mass, which as discussed in Sec.II A can be calculated as the number of reachable events or mean out-component size on the event graph, as which has a power-law asymptote at τ c = 0 of the form M ∼ −τ −1 which confirms the mean-field DP exponent γ = 1 [22].Deriving a closed-form solution for τ becomes prohibitively complex for many types of synthetic temporal networks that involve even the slightest traces of spatial or temporal inhomogeneities and require many simplifying approximations of the structure of networks.As the nature of the assumption are similar to the ones we used while showing the critical exponents in the mean-field regime, this alone would not be productive as a mean to validate or refute the previous theoretical claims for networks with heterogeneous structure or dynamics.Therefore, we complemented these analytical derivations of τ (from Sec.VI 1) and the critical exponents (from the mean-field approach of Ref. [22] and the current section) with measurements derived from simulations.While it would be possible to measure τ from the simulated event graphs, we elected to use δt − δt c as a stand-in for control parameter τ , similar to how p − p c was used in Sec.I B for lattices.Very close to the critical threshold τ → τ c , δt − δt c linearly approximates the control parameter τ , which would preserve the power-law relationships mentioned before at least for some neighbourhood ) finite network for random 9-regular networks in the absorbing phase δt < δtc.Also visualised is the effect of using δt − δtc as an approximation of control parameter τ , which shows similar behavior close to δtc.
of τ = τ c = 0. δt is simply a parameter of the simulation and δt c can be derived empirically for each configuration, either by trial and error or through the finite-size scaling method described in Sec.I B 2. This means that, by virtue of not relying on the methods and assumptions presented previously, we can provide a clean separation between the empirical validation and our theoretical assumptions.
It is possible to find a closed-form solution for τ for very simple systems, such as the case of random k-regular networks with Poisson process link activations.This, however, entails making simplifying assumptions about the structure of the event graph.The results of this derivation and the comparison with empirical measurements follows in Sec.VI 1.For the case of random k-regular static base networks and Poisson process activation of links with mean interevent time λ, we were able to analytically derive a closedform solution of the control parameter τ as a function of δt, k and λ.To this end, it is necessary to derive the joint degree distribution probability-generating function G(x, y) of the event graph based on the excess degree distribution of the base random k-regular network and the Poisson process [22].This leads to a formulation of out-degree and excess out-degree distribution probability-generating functions of the form This in turn, based on relation τ ) Figure 7(a) shows the relationship between the theoretically derived value of the control parameter τ from Eq. ( 13) for different random k-regular networks with a Poisson process with mean inter-event time fixed to 1.As expected, a denser network has a lower onset of criticality in terms of the maximum waiting time δt.Furthermore, a linear approximation of τ ∝ δt − δt c works quite well for these systems for the neighborhood close to τ = 0 given the lower curvature for at least the immediate surrounding of τ c .
Given that, for the event graph representation of an infinite random k-regular networks with a Poisson process activation configuration the out-degree and the excess out-degree distributions are equal, as derived in Eq. ( 12) (i.e.G out 0 (x) = G out 1 (x)), Eq. ( 11) simplifies to M = −τ −1 for τ < 0. Figure 7(b) compares this analytical solution of mean out-component size (calculated with the assumption of the randomness of the event graph) with empirical measurements of a large network.Note that, for k = 9, our best empirical estimate for δt c , δt empirical c = 0.08808, compared to the estimate from the analytical method, δt theoretical c = 0.08559, have a difference of around 3%.This is also visible when comparing empirical measurements of mean cluster mass M (δt) and the theoretical estimations for the system in Fig. 7(b).This can be attributable to the fact that the rate equation Eq. ( 9) is constructed for temporal networks under the assumption that the event graph is indistinguishable from a random directed network with the same joint inand out-degree distribution.This difference seems to suggest that certain local structures in the event graph are very slightly over-represented compared to a random directed graph with the same degree distribution.Also indicated by Fig. 7(b) is the fact that the power-law behavior of the empirical trajectory with a critical exponent of γ = −1 can quite easily be validated by using an empirical estimation of δt c .computational resources provided by the Aalto University Science-IT project and CSC -IT Center for Science, Finland.
The different types of random networks are generated on the fly in the networks/ directory by the scripts provided in the scripts/ directory.Of the real-world datasets, analyzed in [11], the data for the US air transportation network [12] and Helsinki public transportation network [13] are provided along with the implementation.
The code can be compiled using the command make all.Various C++17 features are used extensively throughout the code and it is only tested to compile with GCC 9.3.0,though it is expected to work with a more recent version of GCC as well.
You can generate a random temporal network using the random_network executable: $ ./random_network--seed 0 --nodes 512 --average-degree 9 --static-model regular \ --temporal-model poisson --mean-dt 1 --max-time 2048 > example/example-small.events or by passing a pre-generated static network to the --static-model, instead of specifying one of the existing a static network model.To estimate quantities for source source limited-time spreading process on a synthetic temporal network with δt = 0.8, you can use a variation of the following command: $ ./random_large_single--seed 0 --dt 0.8 --network example/example-small.events\ --summary /dev/stdout --window-min 0 --window-max 2048 \ --time-bins 16 --time-bins-min 0.25 --time-bins-max 2048 which prints the results, including statistics on cluster mass, volume and lifetime and some basic information on the temporal network, in JSON format on standard output.The quantities are presented as average over all events, as well as binned logarithmically based on the time of the event and the time from the event to the end of the observation window.The lifetime and temporal distribution of events that produce "censored" percolation clusters, clusters that might span longer than the window of observation of the temporal network, are reported separately as well to facilitate estimating the probability of survival using Kaplan-Meier estimators [14].Similarly, if the same set of parameters is used with either of the executables random_large_homogeneous or random_large_homogeneous_alternative, it provides information about the distribution of occupied events in case of starting from homogeneous, fully occupied initial condition.The former accepts an additional parameter --field h to simulate an external field of rate h by spontaneously occupying the events, while the latter uses the much faster algorithm presented in this paper which calculates the average effect of occupying exactly one random event.The values reported for dt-band-homogeneous-occupied and dt-band-external-field-occupied corresponds to the number of normally occupied events in a δt band before the end of the observation time window, and the number of occupied events in the presence of a minor external field that spontaneously occupies exactly one random event.The aggregation script aggregate_data.pycan provide more details on how the output of multiple runs of these executables results are aggregated and transformed into quantities that are reported in this manuscript.Results for the synthetic networks were generated through the following iterative system: We started at an initial estimation for δt c and precision of p = 1, where precision indicates the number of digits of certainty for δt c after the decimal point.At each step we calculated all the characteristic quantities for all possible values of δt with one extra digit after decimal point (δt c − 9 × 10 −(p−1) ≥ δt ≥ δt c + 10 × 10 −(p−1) ) and calculate total error of collapse for P (t) ρ(t) V (t) and M (t) trajectories assuming that value is the critical threshold δt c .If one of the δt c candidates has a significantly higher quality of collapse than the current value, δt c is updated, precision is incremented and the process is repeated.After δt c is determined up to the desired level of precision or after determining that the precision cannot be increased further with the current method, system sizes or the number of realizations, we use all the data produced at every step to plot quantities as a function of δt − δt c and as a function of time at δt = δt c .Please note that the main goal of the current manuscript is not to provide the best estimate of δt c for various synthetic temporal networks, but to validate the hypothesis that limited-time spreading on these networks belong to the same universality class as any directed percolation system.It might be possible to design more efficient numerical methods of estimating δt c , but the current method and the estimates presented would suffice for the purposes of this manuscript.

FULL SUIT OF EMPIRICAL RESULTS FOR GRID LATTICE BASED TEMPORAL NETWORKS
Grid lattices is generated as a line, square-shaped, cubic and hyper-cubic with periodic boundary condition.Each system has same linear size across all dimensions, as presented in Tabs.IV, V, VI and VII.

FIG. 1 .
FIG.1.Two spreading scenarios starting from random events (marked with black circles on (b) and (c)) represented over (a) temporal network (b) δt-limited event graph and (c) reduced event graph of a temporal network built from a one-dimensional grid of 40 nodes (displayed on the left side) with Poisson activation of events with mean inter-event time 1 unit of time, simulated for 20 units of time.The adjacency relations have a maximum waiting time δt = 0.8 unit of time.Spatial volume V can be visualized as the mean size of projection of a spreading cluster on the spatial plain, i.e. the static base network on the vertical axis, whereas survival time T is equivalent to the mean size of projection of the cluster on time (horizontal) axis.While measuring a direct analogue to component mass M , integrating the pair-connectedness function across time and space equivalent to the mean sum of lengths of colored horizontal lines in (a) is not straightforward with the event graph representation.It is possible to show that the mean number of uniquely counted events involved in the spreading process, corresponding to the cardinality of the out-component of the initial event here represented by the total number of colored nodes in the event graph, show the same scaling behavior.(d) Homogeneous, fully-occupied initial condition with the occupied events shown in a darker shade than unoccupied events shows the decline and eventual stabilization of the occupation density as time grows.In this scenario, all nodes are considered occupied for time −∞ < t < 0, which translates to the occupation of all events in period 0 ≤ t < δt and all events in their out-components.

Bursty exponent = 2 5 FIG. 2 .
FIG.2.Sample timestamps from a single realization (activations of a single link) with different temporal dynamics.Each point represents a single activation at a specific time.The points are scattered over the vertical axis to avoid overlaps in the visualization.All timestamps were generated for 256 units of time with parameters or minimum cutoffs that would result in an expected inter-event time of 1. Equation (8) defines the parameters and the intensity function of the Hawkes univariate exponential self-exciting process.

FIG. 5 .
FIG.5.Universal scaling functions for δt limited-time reachability over (a,c,e,g,i,k,m,o) random 9-regular network with bursty (heavy-tail with minimum value cutoff) link activation with mean inter-event time of 1 and exponent γ = 2.8 and (b,d,f,h,j,l,n,p) random Erdős-Rényi network ⟨k⟩ = 8 with Hawkes univariate exponential self-exciting process link activation with parameters µ = 0.2, α = 0.8 and θ = 1.0.The finite size scaling is performed for the following single-source scenarios: (a,b) the mean component mass M as function of δt close to critical point and (c,d) as function of time t at the critical point, (e,f) the mean component volume V as function of δt close to critical point and (g,h) as function of time t at the critical point and (o,p) Survival probability P as function of time t at the critical point.For fully-occupied initial conditions the finite-size scaling is performed for (k,l) the occupation density ρ as function of time t at the criticality, and both (i,j) the static density ρstat and (o,p) susceptibility χ as function of δt close to the critical point.The collapse of the universal scaling functions validates the hypothesis that these systems are governed by the same critical exponents as in directed percolation in the mean-field regime.See Sec.III A for the full definitions of the parameters.

FIG. 7 .
FIG. 7. (a) Theoretically derived value of control parameter τ as a function of δt as given in Eq. (13) for random k-regular networks with Poisson processes link activation with λ = 1.The intersection with the horizontal line at τ = 0 indicates the predicted critical value δtc.(b) The analytical solutions for mean out-component size M = (−τ ) −1 as a function of δt compared to empirical measurement of M (δt) over 256 realizations of large (N = 217 ) finite network for random 9-regular networks in the absorbing phase δt < δtc.Also visualised is the effect of using δt − δtc as an approximation of control parameter τ , which shows similar behavior close to δtc.

1 .
Solution for random k-regular static base networks with Poisson link activation

FULL
SUIT OF EMPIRICAL RESULTS FOR RANDOM ERD ŐS-R ÉNYI BASED TEMPORAL NETWORKS

TABLE IV .
Experimental setup for 1-dimensional grid.

TABLE V .
Experimental setup for 2-dimensional square grid.

TABLE VII .
Experimental setup for 4-dimensional nearest-neighbour square grid.