Interacting Hawkes Processes Overview

• Interacting Hawkes processes are counting processes whose intensities depend on historical events via excitation kernels, network weights, and auxiliary state variables.
• They integrate random graph formulations, age-dependent dynamics, and state-dependent feedback to capture complex interactions in neural and social systems.
• The framework offers rigorous well-posedness, propagation of chaos, and long-time asymptotic results that support practical high-dimensional inference and simulation.

Searching arXiv for recent and foundational papers on interacting Hawkes processes. Interacting Hawkes processes are systems of counting processes in which the predictable intensity of each component depends on the past of other components through excitation kernels, network weights, and, in many variants, auxiliary state variables such as age, membrane potential, or latent ancestor labels. In the graph-based formulation of Delattre, Fournier and Hoffmann, each node $i$ has intensity

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$

while later work replaces the complete graph by random possibly diluted and inhomogeneous graphs, introduces age-dependent and variable-length memory mechanisms, and allows state-dependent or ancestor-dependent excitation (Delattre et al., 2014, Agathe-Nerine, 2021, Chevallier, 2015, Ross et al., 4 May 2026). The subject therefore lies at the intersection of interacting particle systems, point-process theory, graph limits, renewal equations, stochastic PDE limits, and high-dimensional statistical inference.

1. Network formulations and model classes

A standard construction uses independent Poisson random measures and thinning. On a directed graph $G=(\mathcal S,\mathcal E)$, a Hawkes process with parameters $(G,\mu,h)$ is defined by

$$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$

This formulation already permits countably many interacting components and nearest-neighbour interaction structures on $\mathbb Z^d$ (Delattre et al., 2014).

Random-graph formulations introduce spatial positions and inhomogeneous connectivity. In the model of neurons located at $x_1,\dots,x_N\in I\subset\mathbb R^d$, one draws independent Bernoulli adjacency variables

$$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$

and defines

$$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$

Here $f$ may be linear, $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$0, or a general Lipschitz nonnegative function; $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$1 is the synaptic memory kernel; and $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$2 is a bounded continuous baseline drive (Agathe-Nerine, 2021).

Mean-field age-dependent models add the predictable age

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$3

to the intensity map. In the Age-Dependent Random Hawkes Process of Chevallier and collaborators, the rate is

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$4

so the process depends simultaneously on time since the last event and on the interaction field generated by the population (Chevallier, 2015).

The state-dependent generalization of Morariu-Patrichi and Pakkanen places the process on a product mark space $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$5, with marks $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$6, and uses the product-form intensity

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$7

A common special case is

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$8

which realizes a coupled Hawkes–Markov chain and recovers the classical multivariate Hawkes process by choosing $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$9 (Morariu-Patrichi et al., 2017).

2. Well-posedness, non-explosion, and exact construction

For finite but heterogeneous networks, existence and uniqueness are typically established by thinning, Picard iteration, and convolution Grönwall estimates. Under Hypothesis 2.3 in the random-graph model, if $G=(\mathcal S,\mathcal E)$0 is Lipschitz, $G=(\mathcal S,\mathcal E)$1, $G=(\mathcal S,\mathcal E)$2 is continuous in $G=(\mathcal S,\mathcal E)$3, Lipschitz in $G=(\mathcal S,\mathcal E)$4 uniformly in $G=(\mathcal S,\mathcal E)$5, and bounded, and the graph weights $G=(\mathcal S,\mathcal E)$6 satisfy uniform bounds, then for each fixed $G=(\mathcal S,\mathcal E)$7 and graph realization there is a unique adapted solution $G=(\mathcal S,\mathcal E)$8 satisfying

$G=(\mathcal S,\mathcal E)$9

Uniqueness is obtained by comparing two solutions through the total-variation difference $(G,\mu,h)$0, while existence follows from Picard iteration (Agathe-Nerine, 2021).

On a possibly infinite directed graph, Delattre–Fournier–Hoffmann assume weighted Lipschitz and integrability conditions involving constants $(G,\mu,h)$1, $(G,\mu,h)$2, and a locally integrable $(G,\mu,h)$3. Their Theorem 6 yields a unique pathwise solution $(G,\mu,h)$4 such that

$(G,\mu,h)$5

again via Picard iteration and convolution-type bounds (Delattre et al., 2014).

Hybrid marked point processes admit a general non-explosion theory. Morariu-Patrichi and Pakkanen prove strong existence and uniqueness for product-form intensities under either a sublinearity condition indexed by a nondecreasing sequence $(G,\mu,h)$6 with $(G,\mu,h)$7, or a Hawkes-type domination condition

$(G,\mu,h)$8

together with

$(G,\mu,h)$9

The construction proceeds by thinning a Poisson random measure $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$0 through

$$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$1

and uniqueness follows by a coupling/enumeration argument (Morariu-Patrichi et al., 2017).

In variable-length memory models, bounded spiking rates supply an especially transparent non-explosion criterion. If $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$2 is nonincreasing and

$$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$3

then $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$4, so each coordinate has finitely many jumps on finite intervals. The same framework underlies the graphical construction and perfect simulation algorithm for the stationary process (Goncalves et al., 2022).

3. Mean-field, graphon, and fluctuation limits

A central theme is propagation of chaos: as the number of interacting components grows, finite subsets become asymptotically independent after coupling with a nonlinear limit process. In the age-dependent mean-field setting, one couples the $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$5-particle system $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$6 to i.i.d. copies $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$7 of a McKean–Vlasov limit driven by the same Poisson random measures and the same past. Under the Lipschitz condition $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$8 and square-integrability assumptions on $$Z^i_t=\int_0^t\int_0^\infty \mathbf{1}_{\{z\le \lambda^i_{s-}\}}\,\pi_i(ds,dz), \qquad \lambda^i_t=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ^j_s.$$9 and $\mathbb Z^d$0, Theorem 4.1 gives

$\mathbb Z^d$1

The limit intensity is deterministic and continuous, and the law of the age process solves the Pakdaman–Perthame–Salort age-structured PDE

$\mathbb Z^d$2

with

$\mathbb Z^d$3

This provides a rigorous micro–macro link between interacting Hawkes particles and an age-PDE description (Chevallier, 2015).

For random spatial graphs, the macroscopic limit is encoded by the nonlinear convolution equation

$\mathbb Z^d$4

where $\mathbb Z^d$5 is the limit graphon and $\mathbb Z^d$6 is the limiting spatial distribution. A quenched coupling uses the same Poisson random measures to build both the finite system and an inhomogeneous Poisson comparison process $\mathbb Z^d$7 with intensity $\mathbb Z^d$8. Under graph-convergence in cut norm and moment bounds on $\mathbb Z^d$9, Theorem 3.5 yields

$x_1,\dots,x_N\in I\subset\mathbb R^d$0

while a stronger operator-norm hypothesis gives uniform convergence in $x_1,\dots,x_N\in I\subset\mathbb R^d$1. The empirical measure

$x_1,\dots,x_N\in I\subset\mathbb R^d$2

also converges in bounded-Lipschitz distance to the law of a Poisson process with intensity $x_1,\dots,x_N\in I\subset\mathbb R^d$3 under $x_1,\dots,x_N\in I\subset\mathbb R^d$4 (Agathe-Nerine, 2021).

Second-order asymptotics are now well developed in several regimes. For mean-field interacting age-dependent Hawkes processes, the fluctuation field

$x_1,\dots,x_N\in I\subset\mathbb R^d$5

converges to a Gaussian limit characterized by a linear stochastic system driven by Gaussian noise rather than Poisson noise. The proof combines Poisson thinning coupling, a Hilbert-space approach with weighted Sobolev spaces, and Rebolledo’s Martingale CLT (Chevallier, 2016).

A distinct limit appears in the diffusive $x_1,\dots,x_N\in I\subset\mathbb R^d$6-regime. Erny, Löcherbach and Loukianova consider a fully symmetric system in which each spike delivers a random kick of size $x_1,\dots,x_N\in I\subset\mathbb R^d$7 to a common membrane potential. The extended generator converges to that of the diffusion

$x_1,\dots,x_N\in I\subset\mathbb R^d$8

so the common intensity converges in distribution in Skorohod space to a CIR-type diffusion. For any fixed $x_1,\dots,x_N\in I\subset\mathbb R^d$9, the point processes $$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$0 converge to conditionally independent counting processes with hazard $$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$1; the authors describe this as conditional propagation of chaos (Erny et al., 2019).

4. Spatial structure, long-time asymptotics, and stability regimes

Long-time behavior depends sharply on the interaction spectrum. In the linear random-graphon model,

$$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$2

the integral operator

$$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$3

has spectral radius

$$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$4

If $$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$5, then there is a unique bounded continuous limit $$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$6 solving

$$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$7

and $$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$8 for each $$\xi_{ij}^{(N)}\sim \mathrm{Bernoulli}(W_N(x_i,x_j)), \qquad A_{ij}^{(N)}:=\kappa_i^{(N)}\xi_{ij}^{(N)},$$9. The limit admits the Neumann-series expansion

$$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$0

If $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$1 and $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$2, then $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$3, recovering the mean-field result; any nonuniform indegree $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$4 imprints spatial variation in the long-time firing rate. By contrast, if $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$5 and $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$6 satisfies a mild irreducibility/primitivity condition, then $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$7 exponentially fast, with growth rate $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$8 determined by $$\lambda_i^{(N)}(t)=f\!\Bigl(u_0(t,x_i)+\frac1N\sum_{j=1}^N A_{ij}^{(N)}\int_0^{t-} h(t-s)\,dZ_j^{(N)}(s)\Bigr).$$9. Spatial inhomogeneity thus enters through the spectrum and eigenfunctions of $f$0 (Agathe-Nerine, 2021).

For random spatial graphs with exponential memory

$f$1

the synaptic current

$f$2

obeys a finite-dimensional Markovian dynamics. If $f$3 is the spectral radius of the operator

$f$4

and $f$5, then the deterministic neural-field ODE has a unique stationary solution $f$6 and the finite system remains close to it over polynomial times. More precisely, for every sufficiently small $f$7 there exist constants $f$8 and a burn-in time $f$9 such that, with probability tending to $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$00 as $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$01,

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$02

The proof balances exponential contraction of the linearized semigroup against concentration inequalities for Bernoulli edges and martingale bounds for the Poisson noise (Agathe-Nerine, 2022).

On the circle, translation symmetry produces a qualitatively different large-time picture. For neurons at equally spaced points on $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$03, cosine interaction, exponential synaptic memory, and sigmoid rate

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$04

the neural-field limit admits a one-parameter manifold

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$05

of stationary solutions. The finite-$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$06 synaptic voltage enters an $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$07-tube around $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$08 by time $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$09 and stays there up to arbitrarily large polynomial times $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$10. On the slower timescale $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$11, the phase $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$12 converges in law to Brownian motion on the circle with diffusion coefficient

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$13

This yields a rigorous phase-diffusion description of wandering bumps (Agathe-Nerine, 2023).

Power-law spatial interactions lead to a different asymptotic regime. In the long-range model on $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$14,

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$15

the subcritical condition $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$16 yields

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$17

If $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$18, there is a unique $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$19 such that $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$20, and

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$21

The same paper states a heuristic dichotomy for fluctuations based on the tail exponent $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$22: for $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$23, the spatial kernel is tied to symmetric $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$24-stable scaling rather than Gaussian scaling (Belmabrouk, 6 Mar 2026).

5. Memory structure, state dependence, and heterogeneous excitation

Interacting Hawkes processes need not depend on the entire past through fixed kernels. In the variable-length memory model, each neuron $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$25 carries a membrane-potential process

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$26

so the relevant history resets at the last spike time. The spiking intensity is

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$27

with $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$28 nonincreasing and bounded between $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$29 and $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$30. A graphical construction splits the dominating Poisson measure into “sure” jumps and “possible” jumps, and perfect simulation is performed through the backward “clan of ancestors.” There exists a critical threshold $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$31 such that if

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$32

then the clan dies out almost surely in finite time, uniformly over $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$33, yielding a unique stationary and ergodic Hawkes process with variable-length memory (Goncalves et al., 2022).

Age dependence and state dependence furnish two distinct mechanisms for memory modulation. In age-dependent Hawkes models, the elapsed time since the last spike enters the intensity map directly and leads, in the mean-field limit, to a nonlinear age-PDE of von Foerster–McKendrick type (Chevallier, 2015). In hybrid marked point processes, by contrast, the marks carry the post-jump state, and the state feeds back into both the event rate and the transition mechanism. Morariu-Patrichi and Pakkanen emphasize feedback loops and regime-switching as characteristic phenomena of this framework, and recover a pure Markov jump chain by setting $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$34 while retaining state-dependent jump probabilities (Morariu-Patrichi et al., 2017).

Two further extensions modify how one event influences later ones. In the two-population model with multiplicative inhibition, the excitatory intensity is

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$35

for the inhibitory population. In the fully coupled case with $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$36, Theorem 3.1 shows that population $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$37 is never supercritical; if $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$38, the long-time behavior is governed by the self-map $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$39, and failure of $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$40 to reduce to a singleton is associated numerically with persistent oscillations (Duval et al., 2021).

The Ancestor Hawkes process of Ross and Deutsch introduces latent labels $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$41 that distinguish immigrant events from triggered events. Its intensity in dimension $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$42 is

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$43

so first-generation and later-generation events can have different excitation matrices. In a standard Hawkes process one would set $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$44; the ancestor formulation separates the two mechanisms. The stability condition is $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$45, and the group-chat application shows that immigrant-message influence $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$46 is systematically larger than triggered-message influence $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$47, whereas a standard Hawkes fit “blurs” these mechanisms (Ross et al., 4 May 2026).

6. Statistical methodology, approximation, and applications

Current statistical work treats interacting Hawkes processes as flexible models for high-dimensional event data rather than only as analytically tractable linear systems. One approximation route starts from a univariate marked Hawkes process and partitions the mark space into $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$48 cells. The resulting $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$49-variate unmarked representation has intensities

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$50

and the marked intensity is reconstructed by

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$51

Under continuity and compactness assumptions, the multivariate representation approximates the true marked intensity in $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$52 as the partition is refined, inherits stationarity when the original process is stationary, and has identifiable conditional-intensity parameters under assumptions (A1)–(A4) (Davis et al., 2024).

For high-dimensional generalized nonlinear Hawkes processes, variational Bayes methods now supply scalable nonparametric inference. The model

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$53

is equipped with spike-and-slab priors on the kernels, basis expansions $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$54, and a mean-field variational approximation. In the sigmoid case, Pólya–Gamma augmentation and a thinned Poisson process yield Gaussian variational factors with closed-form updates, and a two-step sparsity-inducing procedure thresholds the posterior means of the kernel $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$55-norms

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$56

The algorithm is parallelisable and computationally efficient in high-dimensional setting, and the variational posterior concentrates in $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$57-norm at the same rate as the full posterior under the stated prior-mass, entropy, and approximation conditions (Sulem et al., 2022).

Sparse-network estimation with covariates and common drivers is addressed by the model

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$58

Kreiss, Mammen and Polonik estimate $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$59 by penalized least squares with $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$60-penalties on the rows of $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$61, prove first-stage consistency and a de-biased asymptotic normality result

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$62

and obtain per-actor convergence rates for the local parameters under compatibility conditions and sparsity assumptions (Kreiss et al., 4 Apr 2025).

A different methodological issue is incomplete observation. Bayesian inference from aggregated data treats the exact point pattern $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$63 and latent branching structure $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$64 as missing data, and samples from

$$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$65

subject to the observed bin counts. The identifiability results cover temporal, spatio-temporal, and mutually exciting Hawkes processes under general specifications, and the MCMC updates combine conjugate Gibbs steps for $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$66, $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$67, and $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$68 with Metropolis–Hastings updates for nonconjugate parameters and latent event times and locations (Zhou et al., 2022).

Discrete-time interacting Hawkes models can also be fitted nonparametrically. Browning and collaborators model the multivariate discrete-time kernel $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$69 as a random histogram with an unknown number of bins $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$70, unknown knot locations, and unknown heights. Reversible-jump MCMC performs within-model knot moves and birth–death updates for $$\lambda_t^i=\mu_i+\sum_{j\to i}\int_0^{t-} h_{ij}(t-s)\,dZ_s^j,$$71, allowing the triggering kernel to take the form of any step function rather than a fixed parametric family (Browning et al., 2022).

Taken together, these developments show that interacting Hawkes processes now encompass exact and perfect simulation, graphon and mean-field asymptotics, age- and state-structured dynamics, long-time stability and oscillation theory, and statistically identifiable high-dimensional inference under sparsity, aggregation, or nonparametric kernel uncertainty. A plausible implication is that the classical linear complete-graph model is now best viewed as one analytically convenient corner of a much broader interacting point-process theory.