Efficient Gillespie algorithms for spreading phenomena in large and heterogeneous higher-order networks

Higher-order interactions, where groups of nodes interact collectively rather than pairwisely, are central to many complex systems, from neural and ecological networks to social contagion. However, simulating dynamical processes on such higher-order structures remains computationally challenging due to the combinatorial growth of possible interactions. Here, we develop efficient and statistically exact Gillespie algorithms for Markovian spreading dynamics on large and heterogeneous hypergraphs. By incorporating phantom processes  - events that advance time without altering the system's state - , we drastically reduce the computational complexity of standard algorithms ( $O\left(N^2\right)$ ), achieving up to linear scaling with system size. Relying on the susceptible-infected-susceptible model with critical mass thresholds as a benchmark, we show that the optimized algorithms outperform standard approaches by several orders of magnitude, enabling simulations of networks with millions of nodes and broad heterogeneity in both degree and interaction order. Efficient sampling methods, needed to overcome the bottlenecks imposed by either a high maximum order or number of interactions, and other dynamical processes on higher-order networks are tackled. These results establish a general framework for scalable, continuous-time simulations of higher-order contagion and related dynamical processes.