January 01, 1970
We establish a connection between the stochastic self-propagating star-formation model and spatio-temporal point processes by showing that, under suitable discretisation, the SSPSF update law can be represented by a separable spatio-temporal Hawkes process. Building on this connection, we propose a spatio-temporal log-Gaussian Cox-Hawkes process as a continuous point-process model for stochastic star formation. The model represents star-formation events as point patterns driven jointly by deterministic galactic structure, latent spatio-temporal background variation, and dependence on past events. Its key feature is that the deterministic mean field, latent Gaussian random field, and history-dependent interaction field enter through a single log-intensity. This log-scale construction differs from additive Cox-Hawkes formulations and allows the history effect to be signed: past events may either increase or decrease future local intensity while the conditional intensity remains positive. The resulting framework provides an interpretable point-process model for representing latent clustering, self-excitation, local inhibition, and event-driven propagation in stochastic star formation. Beyond linking SSPSF to spatio-temporal point-process theory, it offers a continuous stochastic formulation for analysing the propagation of star formation in galaxies and for interpreting observational surveys of star-forming regions within a unified statistical model.
galaxy formation; galaxy evolution; Cox process; Hawkes process; self-excitation
The stochastic self-propagating star formation (SSPSF) model, introduced by [1] and further developed by [2], [3], provides a probabilistic framework for modelling the formation and propagation of star-forming regions across galactic disks. Unlike traditional density-wave explanations of spiral structure, such as those discussed by [4]–[6], the SSPSF framework explains a range of galactic morphologies through stochastic local triggering, feedback, and differential rotation. In this model, newly formed massive stars or star-forming regions may induce subsequent star formation in neighbouring regions, leading to clustered and self-propagating stellar patterns. Several extensions have refined the original SSPSF mechanism. [7] allowed the triggered star-formation probability to vary with galactic radius, representing the decline of star-formation efficiency in the outer disk. [8] extended the model to three dimensions by stacking multiple two-dimensional layers and adding vertical neighbour connections. Percolation-based extensions proposed by [9] further emphasised the importance of local connectivity and threshold-induced activation. Later, [10] introduced anisotropic propagation rules, in which the triggering probability is distributed elliptically around a star-forming cell, with the orientation and eccentricity linked to the galactic rotation curve. Together, these developments show that SSPSF models provide a flexible rule-based framework for studying how local interactions, stochastic triggering, and feedback processes can generate large-scale star-formation patterns. The SSPSF framework is also conceptually related to sandpile-type models of self-organised criticality. The Abelian sandpile model, introduced by [11] and later developed analytically by [12], shows how simple local rules, threshold activation, and redistribution can generate complex cascading behaviour without global control. In this sense, sandpile models provide a useful analogue for SSPSF, since both frameworks describe how local activation may propagate through neighbouring regions and produce large-scale spatial structure.
Spatio-temporal point processes provide a natural framework for analysing events observed in both space and time [13]–[17]. Building on general point process theory [18], early work developed explicit space-time models for applications such as forest dynamics [19] and introduced second-order methods for detecting space-time interaction and clustering [20]. These developments established spatio-temporal point processes as flexible tools for studying event patterns shaped by spatial structure, temporal variation, and event history [16], [21].
Self-exciting models are widely used for clustered event data. The Hawkes process was introduced as a temporal self-exciting point process by [22], with its cluster representation developed by [23]. It was later extended to spatio-temporal settings, particularly through earthquake models and related branching formulations [13], [24]–[26]. These models decompose event risk into a background component and a history-dependent triggering component, thereby separating exogenous variation from events induced by past events. This interpretable structure has led to broad use in earthquake modelling, crime forecasting, disease modelling, and other clustered event applications [27], [28]. However, classical spatio-temporal Hawkes models are mainly excitatory, since past events are assumed to increase future event risk. They are therefore less suitable when previous events may suppress occurrence nearby in space and time. Several recent works have extended Hawkes processes beyond the purely excitatory setting. [29] considered non-linear Hawkes processes with Gaussian process self-effects, allowing past events to have either excitatory or inhibitory effects. [30] studied log-linear Hawkes processes, where the rate function is exponential and the conditional intensity is a log-linear function of past events. [31] developed inference methods for multivariate exponential Hawkes processes with inhibition. These works study signed or non-linear Hawkes-type dynamics, but they do not consider a model in which a latent spatio-temporal log-Gaussian Cox field and a signed Hawkes-type history field are combined within a single log-intensity.
On the other hand, Log-Gaussian Cox processes [32], [33] provide a complementary framework for modelling clustering through latent background variation. Introduced by [34] as Cox processes with log-Gaussian random intensities, they were later extended to spatio-temporal prediction with stochastic variation over both space and time [35]. Conditional on the latent intensity surface, events follow an inhomogeneous Poisson process, with the log intensity usually represented by covariates and a latent Gaussian field [17], [36]. These models are well suited to clustering driven by unobserved background field [37], but they do not explicitly model event interaction and therefore cannot directly separate latent clustering from triggering by previous events. Neural spatio-temporal point processes further increase modelling flexibility by replacing restrictive parametric assumptions with neural representations. Examples include neural ODE formulations for continuous time event distributions [38], diffusion based models for joint time and location dependence [39], and latent neural intensity models such as DeepSTPP for non-parametric spatio-temporal intensities [40]. These methods can improve prediction and capture complex event patterns, but their learned dependence structures are often less directly interpretable in terms of background variation, excitation, or inhibition.
A closely related development is the Cox-Hawkes model, which combines a stochastic Cox process background with Hawkes-type triggering in order to separate latent spatio-temporal variation from event-induced excitation [41]. However, because the triggering component is usually additive and nonnegative, such models primarily capture excitation and do not naturally represent inhibitory history effects.
While the classical stochastic self-propagating star formation model and spatio-temporal point processes have both been widely studied, their relationship has not been explicitly formulated. In this work, we first build a bridge between the classical SSPSF model and spatio-temporal point processes by showing that, under suitable discretisation, the SSPSF update law can be represented by a separable spatio-temporal Hawkes process. The detailed description of the classical SSPSF model, the relevant spatio-temporal point-process background, the equivalence between the SSPSF update law and a discretised spatio-temporal Hawkes process, and the corresponding discrete likelihood derivation are provided in Appendix 9.
Building on this connection and motivated by the limitations of additive Cox-Hawkes formulations [41], we then propose a new continuous stochastic model for star formation: a spatio-temporal log-Gaussian Cox-Hawkes process with signed history dependence. The proposed model combines deterministic galactic structure, latent spatio-temporal background variation, and history-dependent interaction within a single log-intensity. Unlike additive Cox-Hawkes models, the history effect enters on the log-intensity scale and can therefore be either excitatory or inhibitory while the conditional intensity remains positive. Thus, the model is not simply a continuous version of the classical SSPSF model, but a new point-process framework for modelling stochastic star formation with latent background variation, self-excitation, local inhibition, and event-driven propagation. While our existence and explosion arguments are inspired by the Poisson embedding approach for log-linear Hawkes processes studied by [30], the present model has a different purpose and structure.
The remainder of this work is structured as follows. Section 2 introduces the proposed spatio-temporal log-Gaussian Cox-Hawkes process with signed history dependence. Section 3 discusses its basic properties and special cases. Sections 4 and 5 study existence, non-explosion, and possible explosion of the proposed process. Section 6 introduces the recovery-time version of the model and proves its non-explosion. Section 7 gives a physical interpretation of the recovery-time model for stochastic star formation. Section 8 concludes with a discussion and outlook. Appendix 9 provides the detailed background on the classical SSPSF model, spatio-temporal point processes, the connection between the SSPSF update law and a discretised spatio-temporal Hawkes process, and the corresponding discrete likelihood derivation.
Let \(\mathcal{S}\subset\mathbb{R}^2\) be a bounded spatial domain and let \([0,T]\) be the observation window. We observe a spatio-temporal point pattern \(\mathcal{W}=\{(s_i,t_i):i=1,\ldots,n\}\), where \(s_i\in\mathcal{S}\) and \(0<t_1<...<t_n\leq T\). Let \(\mathcal{H}_t=\{(s_i,t_i):t_i<t\}\) denote the history before time \(t\). We model the process through its conditional intensity \(\lambda(s,t\mid\mathcal{H}_t)\) and propose a log-Gaussian Cox-Hawkes process of the form \[\lambda(s,t\mid\mathcal{H}_t) = \exp\{\mu(s,t)+Y(s,t)+D(s,t\mid\mathcal{H}_t)\},\] where \(\mu(s,t)\) is a deterministic mean field, \(Y(s,t)\) is a latent spatio-temporal Gaussian random field, and \(D(s,t\mid\mathcal{H}_t)\) is a history induced dynamic field. The latent field captures residual spatio-temporal variation, while the dynamic field captures the effect of previous events on future event intensity. We assume that \(Y(s,t)\) is a mean zero Gaussian process, \(Y(s,t)\sim \mathrm{GP}\{0,C_Y\}\), with separable exponential covariance function \[C_Y((s_1,t_1),(s_2,t_2)) = \sigma_Y^2 \exp\left\{ -\frac{\|s_2-s_1\|}{\phi_Y} -\theta_Y |t_2-t_1| \right\}.\] Here \(\sigma_Y^2>0\) is the marginal variance, \(\phi_Y>0\) controls the spatial correlation range, and \(\theta_Y>0\) controls the temporal rate of decorrelation. The dynamic field is defined as \(D(s,t\mid\mathcal{H}_t)=\sum_{i:t_i<t}g(t-t_i,s-s_i)\), where \(g(u,h)\) is a causal spatio-temporal interaction function satisfying \(g(u,h)=0\) for \(u\leq 0\). We use the separable specification \[g(u,h) = \alpha \exp(-\beta u)\varphi_{\Sigma}(h)\mathbf{1}\{u>0\},\] where \(\alpha\in\mathbb{R}\) controls the direction and magnitude of the history effect, \(\beta>0\) controls temporal decay, and \(\varphi_{\Sigma}(h)\) is a Gaussian spatial kernel, \[\varphi_{\Sigma}(h) = \frac{1}{2\pi|\Sigma|^{1/2}} \exp\left\{ -\frac{1}{2}h^\top\Sigma^{-1}h \right\},\] with positive definite covariance matrix \(\Sigma\). In this work, we define the spatial covariance matrix in the interaction kernel as \[\Sigma = \sigma_\Sigma^2 \begin{pmatrix} 1 & \eta_\Sigma\\ \eta_\Sigma & 1 \end{pmatrix},\] where \(\sigma_\Sigma>0\) controls the marginal spatial scale of the history-dependent interaction and \(\eta_\Sigma\in(-1,1)\) controls the correlation between the two spatial coordinate directions. When \(\eta_\Sigma=0\), this reduces to the isotropic specification \(\Sigma=\sigma_\Sigma^2 I_2\). Hence, \[D(s,t\mid\mathcal{H}_t) = \alpha \sum_{i:t_i<t} \exp\{-\beta(t-t_i)\} \varphi_{\Sigma}(s-s_i),\] and the full conditional intensity becomes \[\lambda(s,t\mid\mathcal{H}_t) = \exp\left[ \mu(s,t) + Y(s,t) + \alpha \sum_{i:t_i<t} \exp\{-\beta(t-t_i)\} \varphi_{\Sigma}(s-s_i) \right].\]
For a spatio-temporal point process observed on \(\mathcal{S}\times[0,T]\), the conditional log likelihood associated with the intensity \(\lambda(s,t\mid\mathcal{H}_t)\) is \[\ell = \sum_{i=1}^{n} \log \lambda(s_i,t_i\mid\mathcal{H}_{t_i}) - \int_0^T\int_{\mathcal{S}} \lambda(s,t\mid\mathcal{H}_t)\,ds\,dt.\] Under the proposed log-Gaussian Cox-Hawkes model, conditional on the latent field \(Y(s,t)\), this becomes \[\begin{align} \ell(\Theta) &= \sum_{i=1}^{n} \left[ \mu(s_i,t_i) + Y(s_i,t_i) + \alpha \sum_{j:t_j<t_i} \exp\{-\beta(t_i-t_j)\} \varphi_{\Sigma}(s_i-s_j) \right] \\ &\quad - \int_0^T\int_{\mathcal{S}} \exp\left[ \mu(s,t) + Y(s,t) + \alpha \sum_{j:t_j<t} \exp\{-\beta(t-t_j)\} \varphi_{\Sigma}(s-s_j) \right] \,ds\,dt . \end{align}\] Let \(\gamma_\mu\) denote the finite dimensional parameters in the deterministic mean field \(\mu(s,t)\). We collect the unknown parameters as \(\Theta = \left( \gamma_\mu, \alpha, \beta, \sigma_\Sigma, \eta_\Sigma, \sigma_Y, \phi_Y, \theta_Y \right),\) where \(\alpha\in\mathbb{R}\) controls the direction and strength of the history effect, and \(\beta>0\) controls its temporal decay. The parameters \(\sigma_\Sigma\) and \(\eta_\Sigma\) determine the spatial covariance matrix \(\Sigma\) in the interaction kernel. The parameters \(\sigma_Y>0\), \(\phi_Y>0\), and \(\theta_Y>0\) determine the covariance structure of the latent Gaussian field, controlling its marginal variability, spatial correlation range, and temporal decorrelation rate, respectively.
The conditional intensity is defined as
\[\lambda(s,t \mid \mathcal{H}_t) = \exp \{\mu(s,t)+Y(s,t)+D(s,t \mid \mathcal{H}_t)\}.\]
Since the exponential function is strictly positive, the conditional intensity remains positive for all parameter values. This property holds regardless of whether the history-dependent effect is positive or negative.
When \(\alpha = 0,\) the history-dependent field vanishes, \(D(s,t \mid \mathcal{H}_t)=0,\) and the model reduces to \[\lambda(s,t) = \exp \{\mu(s,t)+Y(s,t)\}.\]
This is the standard spatio-temporal log-Gaussian Cox process. In this case, clustering is entirely driven by the latent Gaussian field and the deterministic mean structure.
When the latent Gaussian field is absent, \(Y(s,t)=0,\) the conditional intensity becomes
\[\lambda(s,t \mid \mathcal{H}_t) = \exp \{\mu(s,t)+D(s,t \mid \mathcal{H}_t) \}.\]
In this case, the model reduces to a spatio-temporal log-linear Hawkes-type process. The dependence on past events enters through the log-intensity rather than additively on the intensity scale. The process no longer contains latent background variation, and local dependence is generated by the history-dependent field.
The sign of the parameter \(\alpha\) determines the direction of the history effect. When \(\alpha > 0,\) previous events increase the future local log-intensity, producing self-excitation. When \(\alpha < 0,\) previous events decrease the future local log-intensity, producing local inhibition. When \(\alpha = 0,\) the process contains no history dependence and reduces to a log-Gaussian Cox process.
Cox-Hawkes models typically take the additive form \[\lambda(s,t\mid\mathcal{H}_t) = \exp\{f(s,t)\} + \sum_{i:t_i<t}g(t-t_i,s-s_i).\]
In these models, the Cox background and the Hawkes triggering term are added on the intensity scale, and the triggering function is usually nonnegative. In contrast, the proposed model places the deterministic mean field, latent Gaussian field, and history-induced dynamic field inside a single log-intensity. As a result, excitation and inhibition can be represented within the same framework, while positivity of the intensity is automatically preserved.
Following [30], this section gives a Poisson embedding construction of the proposed process, conditional on a realisation of the latent Gaussian field. This construction shows that the model is well defined up to its explosion time, \(T_{\infty}=\lim_{n\to\infty}T_n.\)
The idea is to introduce a background homogeneous Poisson random measure on the three-dimensional space \([0,\infty)\times S\times \mathbb{R}_{\geq 0},\) with coordinates \((t,s,z)\), where \(t\) is the time coordinate, \(s\) is the spatial coordinate, and \(z\) is an auxiliary height coordinate. The Poisson points are spread uniformly in this enlarged space. The proposed point process is then obtained by accepting only those points whose height satisfies \(z\leq \lambda(s,t).\) Thus, the Poisson embedding construction turns the intensity function \(\lambda(s,t)\) into an acceptance threshold for points from the background Poisson random measure.
Let \(S\subset\mathbb{R}^2\) be bounded, and let \(Y(s,t)=y(s,t)\) be a fixed realisation of the latent Gaussian field. A spatio-temporal point process \(N\) on \(S\times[0,\infty)\) with explosion time \(T_{\infty}\in(0,\infty]\) is called a log-Gaussian Cox-Hawkes process with signed history dependence if its conditional intensity is \[\lambda^y(s,t\mid\mathcal{H}_t) = \exp\{\mu(s,t)+y(s,t)+D(s,t\mid\mathcal{H}_t)\} \mathbf{1}\{t<T_{\infty}\},\] where \[D(s,t\mid\mathcal{H}_t) = \alpha \sum_{i:t_i<t} \exp\{-\beta(t-t_i)\} \varphi_{\Sigma}(s-s_i).\] The process is called non-explosive if \(T_{\infty}=\infty\) almost surely, and explosive if \(T_{\infty}<\infty\) with positive probability.
Let \(\mathcal{N}\) be a homogeneous Poisson random measure on \([0,\infty)\times S\times \mathbb{R}_{\geq 0}\) with intensity measure \(dt\,ds\,dz\). Thus, the expected number of Poisson points in any measurable region is equal to the Lebesgue measure of that region.
Let \((\mathcal{F}_t^{\mathcal{N}})_{t\geq 0}\) be the filtration generated by the Poisson random measure up to time \(t\), where \[\mathcal{F}_t^{\mathcal{N}} = \sigma\left( \mathcal{N}((a,b]\times B\times C): 0<a<b\leq t,\;B\in\mathcal{B}(S),\;C\in\mathcal{B}(\mathbb{R}_{\geq 0}) \right).\] This \(\sigma\)-field contains all information from the background Poisson process up to time \(t\), namely all candidate points whose time coordinate is at most \(t\). Let \(\mathcal{G}_t = \sigma(\mathcal{F}_t^{\mathcal{N}},\mathcal{F}_0^N), \: t\geq0,\) where \(\mathcal{F}_0^N\) contains the initial history.
Let \(N\) be a spatio-temporal point process on \(S\times\mathbb{R}\) with initial condition \(N^0\) on \(S\times\mathbb{R}_{\leq 0}\), independent of \(\mathcal{N}\). Let \[\mathcal{F}_t^N = \sigma\{N(A):A\in\mathcal{B}(S\times(-\infty,t])\}\]be the natural filtration of \(N\). The process \(N\) has two parts. The first part is the initial history \(N^0\) before time zero. The second part consists of future events on \(S\times(0,\infty)\), which will be constructed from the background Poisson random measure \(\mathcal{N}\). The independence assumption means that the past history before time zero is independent of the future Poisson randomness used to generate new events. Suppose that the intensity \(\lambda(s,t)\) is a nonnegative predictable process with respect to the filtration \((\mathcal{G}_t)_{t\geq0}\). Predictability means that, at time \(t\), the value of \(\lambda(s,t)\) is determined by information available just before \(t\). It cannot depend on future Poisson points. Otherwise, the acceptance rule would be allowed to look into the future, and it would not define a valid point process intensity.
For the future part of the process, assume that \[N(A) = \int_0^\infty \int_S \int_{\mathbb{R}_{+}} \mathbf{1}\{(s,t)\in A\} \mathbf{1}\{z\leq \lambda(s,t)\} \mathcal{N}(dt,ds,dz),\] for all \(A\in\mathcal{B}(S\times[0,\infty))\).
Here, \(N(A)\) counts the background Poisson points whose space-time coordinates lie in \(A\) and whose auxiliary height \(z\) is below the threshold \(\lambda(s,t)\). Equivalently, the process \(N\) is obtained by accepting candidate points \((t,s,z)\) whenever \(z\leq \lambda(s,t).\) Thus, \(\lambda(s,t)\) acts as an acceptance surface over space and time.
Then \(\lambda(s,t)\) is the conditional intensity of \(N\) with respect to \((\mathcal{G}_t)_{t\geq0}\). If \(N\) is constructed by accepting Poisson points below the predictable threshold \(\lambda(s,t)\), then \(\lambda(s,t)\) is indeed the conditional intensity of the accepted process. In short, Poisson embedding together with a predictable acceptance rule gives an accepted point process with intensity \(\lambda(s,t)\). Furthermore, if \(\lambda(s,t)\) is predictable with respect to the natural filtration \((\mathcal{F}_t^N)_{t\geq 0}\) of the accepted process \(N\), then \(\lambda(s,t)\) is the conditional intensity of \(N\) with respect to \((\mathcal{F}_t^N)_{t\geq0}\).
There are two information structures. The first one is the larger filtration generated by the background Poisson random measure and the initial condition \(N^0\). The second one is the smaller natural filtration generated only by the accepted process \(N\). The first part of the lemma shows that \(\lambda(s,t)\) is an intensity with respect to the larger filtration. If \(\lambda(s,t)\) can be computed from the history of \(N\) itself, then it is also an intensity with respect to the natural history of \(N\). This is important because the proposed intensity should depend only on the observed history of the process, not on rejected candidate Poisson points.
Proof. Recall that \(\mathcal{G}_t=\sigma(\mathcal{F}_t^{\mathcal{N}},\mathcal{F}_0^N)\) is the filtration generated by the background Poisson random measure up to time \(t\), and \(\mathcal{F}_0^N\) contains the initial history.
The first statement follows from the standard compensator characterisation of point processes constructed from Poisson random measures with predictable acceptance functions [30], [42]. The process \(N\) is obtained by accepting background Poisson candidate points \((t,s,z)\) whenever \(z\leq \lambda(s,t)\). Since the acceptance function is predictable with respect to \((\mathcal{G}_t)_{t\geq0}\), the compensator of the accepted random measure is \(\lambda(s,t)\,dt\,ds\). Therefore, for all \(B\in\mathcal{B}(S)\) and \(0\leq a<b\), \[\mathbb{E}\{N(B\times(a,b])\mid \mathcal{G}_a\} = \mathbb{E} \left[ \int_a^b \int_B \lambda(s,t)\,ds\,dt \mid \mathcal{G}_a \right].\]
It remains to show the second statement. Since \(N\) is constructed from \(\mathcal{N}\), we have \(\mathcal{F}_a^N\subseteq \mathcal{G}_a\). Using the tower property of conditional expectation, for \(0\leq a<b\) and \(B\in\mathcal{B}(S)\), \[\mathbb{E}\{N(B\times(a,b])\mid \mathcal{F}_a^N\} = \mathbb{E} \left[ \mathbb{E}\{N(B\times(a,b])\mid \mathcal{G}_a\} \mid \mathcal{F}_a^N \right].\] By the first part of the lemma, \[\mathbb{E}\{N(B\times(a,b])\mid \mathcal{G}_a\} = \mathbb{E} \left[ \int_a^b \int_B \lambda(s,t)\,ds\,dt \mid \mathcal{G}_a \right].\] Hence \[\mathbb{E}\{N(B\times(a,b])\mid \mathcal{F}_a^N\} = \mathbb{E} \left[ \mathbb{E} \left[ \int_a^b \int_B \lambda(s,t)\,ds\,dt \mid \mathcal{G}_a \right] \mid \mathcal{F}_a^N \right].\] Applying the tower property again gives \[\mathbb{E}\{N(B\times(a,b])\mid \mathcal{F}_a^N\} = \mathbb{E} \left[ \int_a^b \int_B \lambda(s,t)\,ds\,dt \mid \mathcal{F}_a^N \right].\] Therefore \(\lambda(s,t)\) is the conditional intensity of \(N\) with respect to the natural filtration \((\mathcal{F}_t^N)_{t\geq0}\). This step moves from the larger information structure \(\mathcal{G}_a\), which contains the background Poisson information and the initial condition, to the smaller natural filtration \(\mathcal{F}_a^N\), which contains only the history of accepted events up to time \(a\). It shows that, when \(\lambda(s,t)\) is predictable from the accepted process itself, the same function is also the intensity with respect to the natural history of \(N\). ◻
Let \(S\subset \mathbb{R}^2\) be bounded, and let \(y(s,t)\) be a fixed realisation of the latent Gaussian field \(Y(s,t)\). Suppose that \[b(s,t)=\mu(s,t)+y(s,t)\] is Borel measurable and locally bounded on \(S\times[0,\infty)\). Let \(N^0\) be an initial spatio-temporal point process supported on \(S\times(-\infty,0]\), independent of the background Poisson random measure \(\mathcal{N}\). Represent the initial process \(N^0\) by the points \(\{(S_i^0,T_i^0):i\leq 0\}, \: T_i^0\leq 0.\) For \(t>0\), define the initial history by \[\mathcal{H}_t^{(0)} = \{(S_i^0,T_i^0):i\leq 0,\;T_i^0<t\}.\] Assume the initial history-dependent part \[D_0(s,t) = D(s,t\mid\mathcal{H}_t^{(0)})\] is finite and locally bounded on \(S\times[0,\infty)\). For a history \(\mathcal{H}_t=\{(s_i,t_i):t_i<t\}\) for which the following sum is well defined, define the history-dependent component \[D(s,t\mid\mathcal{H}_t) = \alpha \sum_{i:t_i<t} \exp\{-\beta(t-t_i)\} \varphi_{\Sigma}(s-s_i),\] and the conditional intensity \[\lambda^y(s,t\mid\mathcal{H}_t) = \exp\{b(s,t)+D(s,t\mid\mathcal{H}_t)\}.\] Let \(\mathcal{N}\) be a homogeneous Poisson random measure on \([0,\infty)\times S\times\mathbb{R}_{\geq0}\) with intensity measure \(dt\:ds\:dz\), independent of \(N^0\). We construct the future events after time zero recursively. Set \(T_0=0\). For \(u>0\), define the initial conditional intensity \[\lambda_0(s,u) = \lambda^y(s,u\mid\mathcal{H}_u^{(0)}) = \exp\{b(s,u)+D_0(s,u)\}.\] The first future event time is \[T_1 = \inf \left\{ t>0: \int_{(0,t]\times S\times\mathbb{R}_{\geq0}} \mathbf{1}\{z\leq\lambda_0(s,u)\} \mathcal{N}(du,ds,dz)>0 \right\}.\] If \(T_1<\infty\), let \(S_1\) be the spatial coordinate of the corresponding accepted point. Suppose that \((S_1,T_1),\ldots,(S_n,T_n)\) have been constructed, with \(0<T_1<\cdots<T_n<\infty.\) For \(u>T_n\), define the history \[\mathcal{H}_u^{(n)} = \mathcal{H}_u^{(0)} \cup \{(S_i,T_i):i=1,\ldots,n,\;T_i<u\}.\] Since \(u>T_n\), all first \(n\) future events are included in the history. Hence \[D(s,u\mid\mathcal{H}_u^{(n)}) = D_0(s,u) + \alpha \sum_{i=1}^{n} \exp\{-\beta(u-T_i)\} \varphi_{\Sigma}(s-S_i).\] Set \[\lambda_n(s,u) = \lambda^y(s,u\mid\mathcal{H}_u^{(n)}) = \exp\{b(s,u)+D(s,u\mid\mathcal{H}_u^{(n)})\}.\] The next future event time is \[T_{n+1} = \inf \left\{ t>T_n: \int_{(T_n,t]\times S\times\mathbb{R}_{\geq0}} \mathbf{1}\{z\leq\lambda_n(s,u)\} \mathcal{N}(du,ds,dz)>0 \right\}.\] If \(T_{n+1}<\infty\), let \(S_{n+1}\) be the spatial coordinate of the corresponding accepted point. If \(T_{n+1}=\infty\), set \(T_k=\infty\) for all \(k>n+1\), and the construction stops. Recall that the explosion time is \(T_{\infty} = \lim_{n\to\infty}T_n.\) The constructed counting measure is \[N(A) = N^0(A\cap(S\times(-\infty,0])) + \sum_{i:T_i<T_{\infty}} \mathbf{1}\{(S_i,T_i)\in A\}, \: A\in\mathcal{B}(S\times\mathbb{R}).\] Then, conditional on \(Y(s,t)=y(s,t)\), the constructed process is a spatio-temporal point process up to the explosion time \(T_{\infty}\). For \(t<T_{\infty}\), its conditional intensity is \[\lambda^y(s,t\mid\mathcal{H}_t) = \exp\{\mu(s,t)+y(s,t)+D(s,t\mid\mathcal{H}_t)\}.\]
Proof. Before the explosion time, only finitely many future events have been constructed. On each interval \((T_n,T_{n+1})\), the past events consist of the initial history together with the first \(n\) future events. Hence the history used by the constructed process is \(\mathcal{H}_u^{(n)}\). Therefore, \[\lambda^y(s,u\mid\mathcal{H}_u) = \lambda_n(s,u), \: u\in(T_n,T_{n+1}).\]
We first check that the recursive construction is locally well defined before explosion. The Gaussian spatial kernel is \[\varphi_{\Sigma}(h) = \frac{1}{2\pi|\Sigma|^{1/2}} \exp\left\{ -\frac{1}{2}h^\top\Sigma^{-1}h \right\}.\] Its maximum is attained at \(h=0\), and hence \[C_{\Sigma} = \sup_{h\in\mathbb{R}^2}\varphi_{\Sigma}(h) = \varphi_{\Sigma}(0) = \frac{1}{2\pi|\Sigma|^{1/2}} <\infty.\] Since \(\beta>0\) and \(u>T_i\) for each previously observed future event, we have \[\exp\{-\beta(u-T_i)\}\leq 1.\] Hence, for any finite \(t>T_n\) and any \(u\in[T_n,t]\), \[\left| \alpha \sum_{i=1}^{n} \exp\{-\beta(u-T_i)\} \varphi_{\Sigma}(s-S_i) \right| \leq |\alpha|\,n\,C_{\Sigma}.\] Since \(b(s,u)\) and \(D_0(s,u)\) are locally bounded on \(S\times[0,\infty)\), there exists a finite constant \(Q_{n,t}\) such that \[b(s,u)+D_0(s,u) \leq Q_{n,t}, \: (s,u)\in S\times[T_n,t].\] Therefore, we have \[\lambda_n(s,u) \leq \exp\{Q_{n,t}+|\alpha|\,n\,C_{\Sigma}\}.\] As \(S\) is bounded, \[\int_{T_n}^{t}\int_S \lambda_n(s,u)\:ds\:du \leq |S|(t-T_n) \exp\{Q_{n,t}+|\alpha|\,n\,C_{\Sigma}\} <\infty.\] Thus the accepted region under \(\lambda_n\) has finite Poisson measure on every finite interval before the next event time.
Because time is continuous, the probability that two background Poisson candidate points have exactly the same time coordinate is zero. Therefore, two accepted candidate points also do not occur at exactly the same time with probability one. Hence, whenever \(T_{n+1}<\infty\), there is almost surely a unique accepted point at time \(T_{n+1}\), and its spatial coordinate \(S_{n+1}\) is well defined.
For every \(t<T_{\infty}\), the piecewise recursive construction gives \[N(A\cap(S\times(0,t])) = \int_0^t \int_S \int_{\mathbb{R}_{\geq0}} \mathbf{1}\{(s,u)\in A\} \mathbf{1}\{z\leq\lambda^y(s,u\mid\mathcal{H}_u)\} \mathcal{N}(du,ds,dz).\] The integrand depends only on events before time \(u\), and is therefore predictable with respect to the natural history of the constructed process. By Lemma 4.2, the accepted point process has conditional intensity \[\lambda^y(s,t\mid\mathcal{H}_t) = \exp\{\mu(s,t)+y(s,t)+D(s,t\mid\mathcal{H}_t)\}\] for \(t<T_{\infty}\). ◻
We now give a simple sufficient condition for non-explosion (\(T_{\infty} = \infty\)) in the non-excitatory case, where \(\alpha\leq0\). The result is again stated conditional on a realisation of the latent Gaussian field.
Let \(S\subset\mathbb{R}^2\) be bounded, and let \(y(s,t)\) be a fixed realisation of the latent Gaussian field. Suppose that \(b(s,t)=\mu(s,t)+y(s,t)\) is Borel measurable and locally bounded on \(S\times[0,\infty)\). Let the process be constructed by the Poisson embedding construction of Lemma 4.3, with an initial history for which the history-dependent field is well defined. If \(\alpha\leq0,\) then the constructed process is non-explosive. That is, \(T_{\infty}=\infty\) almost surely, given the fixed realisation \(Y(s,t)=y(s,t)\).
Proof. Since the Gaussian spatial kernel is nonnegative and the temporal decay term is positive, we have \[\exp\{-\beta(t-t_i)\}>0, \quad \varphi_{\Sigma}(s-s_i)\geq0.\] Therefore, when \(\alpha\leq0\), every event in the history gives a non-positive contribution to the history-dependent field. Hence, the history-dependent component \[D(s,t\mid\mathcal{H}_t) = \alpha \sum_{i:t_i<t} \exp\{-\beta(t-t_i)\} \varphi_{\Sigma}(s-s_i) \leq0.\] It follows that the conditional intensity satisfies \[\lambda^y(s,t\mid\mathcal{H}_t) = \exp\{b(s,t)+D(s,t\mid\mathcal{H}_t)\} \leq \exp\{b(s,t)\}.\]
Fix any finite time \(T>0\). Since \(b(s,t)\) is locally bounded on \(S\times[0,T]\) and \(S\) is bounded, we have \[\int_0^T\int_S \exp\{b(s,t)\}\,ds\,dt <\infty.\] Using the same Poisson random measure \(\mathcal{N}\) as in the embedding construction, define \[\widetilde{N}_T = \int_0^T \int_S \int_{\mathbb{R}_{\geq0}} \mathbf{1}\{z\leq \exp\{b(s,t)\}\} \mathcal{N}(dt,ds,dz).\] Then \(\widetilde{N}_T\) is a Poisson random variable with finite mean \(\int_0^T\int_S \exp\{b(s,t)\}\,ds\,dt.\) Thus \(\widetilde{N}_T<\infty\), then the constructed process is non-explosive conditional on \(Y(s,t)=y(s,t)\). Since \[\lambda^y(s,t\mid\mathcal{H}_t) \leq \exp\{b(s,t)\},\] every accepted future point of the proposed process on \(S\times(0,T]\) is also contained in the dominating Poisson process. Therefore, \[N(S\times(0,T]) \leq \widetilde{N}_T <\infty\] almost surely given \(Y(s,t)=y(s,t)\). Then, we need to connect this finite window bound to the explosion time. If \(T_{\infty}<\infty\), then infinitely many future events must occur before some finite integer time \(k\). Hence \[\{T_{\infty}<\infty\} \subseteq \bigcup_{k=1}^{\infty} \{N(S\times(0,k])=\infty\}.\]
For each integer \(k\), the argument above gives \[N(S\times(0,k]) \leq \widetilde{N}_k <\infty\] almost surely conditional on \(Y(s,t)=y(s,t)\). Hence \[\mathbb{P}\left(N(S\times(0,k])=\infty\mid Y(s,t)=y(s,t)\right)=0\] for every \(k\). Therefore, \[\mathbb{P}\left(T_{\infty}<\infty\mid Y(s,t)=y(s,t)\right)=0,\] and hence \[T_{\infty}=\infty\] almost surely given \(Y(s,t)=y(s,t)\). ◻
We now show that, in contrast to the non-excitatory case, the excitatory case may lead to explosion. The argument follows the same intuition as the positive near-zero memory case in [30] that repeated events occurring in a sufficiently small time window can increase the log-intensity fast enough to produce infinitely many events in finite time.
Suppose that \(\alpha>0\). Let \(B\subset S\) be a measurable spatial set with \(|B|>0\). Assume that there exists \(\delta>0\) such that \(m= \inf_{(s,t)\in B\times[0,\delta]} b(s,t) > -\infty.\) Let \(c_B = \inf_{s,s'\in B} \varphi_{\Sigma}(s-s').\) Since \(\varphi_{\Sigma}\) is a positive Gaussian kernel and \(B\) is bounded, we have \(c_B>0\). Then, conditional on \(Y(s,t)=y(s,t)\), the process has positive probability of explosion before time \(\delta\). That is, \[\mathbb{P}\left(T_{\infty}\leq \delta \mid Y(s,t)=y(s,t)\right)>0.\]
Proof. Let \(N_B(t)=N(B\times(0,t])\) be the number of events occurring in the spatial region \(B\) up to time \(t\). Let \(\tau_k\) be the time of the \(k\)-th event occurring in \(B\), and set \(\tau_0=0\). For \(s,s_i\in B\) and \(0<t-t_i\leq\delta\), we have \[\exp\{-\beta(t-t_i)\} \geq \exp\{-\beta\delta\}\] and the Gaussian kernel \(\varphi_{\Sigma}(s-s_i) \geq c_B.\) Then, each previous event in \(B\) occurring within the time interval \((0,\delta]\) contributes at least \[a = \alpha \exp\{-\beta\delta\}c_B > 0\] to the history-dependent component. On the event that \(N_B(t)=n\) for some \(t\leq\delta\), all \(n\) previous events in \(B\) contribute non-negatively, and the events in \(B\) alone give the lower bound \[D(s,t\mid\mathcal{H}_t) \geq an, \: s\in B.\] Since \(b(s,t)\geq m\) on \(B\times[0,\delta]\), it follows that the total conditional rate of events in \(B\) satisfies \[\int_B \lambda^y(s,t\mid\mathcal{H}_t)\,ds \geq |B|\exp\{m+an\}.\]
We need to show there is a positive probability that the waiting times between events become so small that infinitely many events occur before time \(\delta\).
Choose \(\varepsilon>0\) such that \(\varepsilon \sum_{k=1}^{\infty} \frac{1}{k^2} < \delta.\) Define \[E_n = \bigcap_{k=1}^{n} \left\{ \tau_k-\tau_{k-1} < \frac{\varepsilon}{k^2} \right\}.\] On \(E_n\), the first \(n\) waiting times between successive events in \(B\) are all sufficiently small. In particular, on \(E_{k-1}\), we have \(\tau_{k-1}<\delta\), and before the \(k\)-th event in \(B\), there are already \(k-1\) events in \(B\). Hence the conditional rate of events in \(B\) is bounded from below by \[r_k = |B|\exp\{m+a(k-1)\}.\] Thus, before the \(k\)-th event, the process in \(B\) is at least as active as a Poisson process with rate \(r_k\). Therefore, using the Poisson embedding construction and the probability of at least one Poisson point in a time interval, the probability that the next event in \(B\) occurs within time \(\varepsilon/k^2\) is bounded below by \[\mathbb{P}\left(\tau_k-\tau_{k-1}<\frac{\varepsilon}{k^2}\mid E_{k-1}\right)\geq 1-\exp\left\{-r_k\frac{\varepsilon}{k^2}\right\}.\] Substituting \(r_k = |B|\exp\{m+a(k-1)\}\) gives \[\mathbb{P}\left(\tau_k-\tau_{k-1}<\frac{\varepsilon}{k^2}\mid E_{k-1}\right)\geq 1-\exp\left\{-\frac{\varepsilon}{k^2}|B|\exp\{m+a(k-1)\}\right\}.\] Consequently, \[\mathbb{P}\left(\bigcap_{k=1}^{\infty}\left\{\tau_k-\tau_{k-1}<\frac{\varepsilon}{k^2}\right\}\mid Y(s,t)=y(s,t)\right) \geq \prod_{k=1}^{\infty}\left[1-\exp\left\{-\frac{\varepsilon}{k^2}|B|\exp\{m+a(k-1)\}\right\}\right].\] The infinite product is strictly positive because \[\sum_{k=1}^{\infty}\exp\left\{-\frac{\varepsilon}{k^2}|B|\exp\{m+a(k-1)\}\right\}<\infty.\] Hence the event \(\bigcap_{k=1}^{\infty}\left\{\tau_k-\tau_{k-1}<\frac{\varepsilon}{k^2}\right\}\) has positive conditional probability, that is \[\mathbb{P}\left(\bigcap_{k=1}^{\infty}\left\{\tau_k-\tau_{k-1}<\frac{\varepsilon}{k^2}\right\}\mid Y(s,t)=y(s,t)\right) > 0.\]
Let \(A=\bigcap_{k=1}^{\infty}\left\{\tau_k-\tau_{k-1}<\frac{\varepsilon}{k^2}\right\}.\) The argument above shows that \(\mathbb{P}(A\mid Y(s,t)=y(s,t))>0\). On \(A\), the time of the \(n\)-th event in \(B\) satisfies \(\tau_n = \sum_{k=1}^{n}(\tau_k-\tau_{k-1}) < \sum_{k=1}^{n}\frac{\varepsilon}{k^2}.\) Letting \(n\to\infty\) gives \[\lim_{n\to\infty}\tau_n \leq \sum_{k=1}^{\infty}\frac{\varepsilon}{k^2} < \delta.\] So on \(A\), infinitely many events occur in \(B\times(0,\delta]\). Therefore \(A\subseteq\{T_{\infty}\leq\delta\}\), and so \[\mathbb{P}\left(T_{\infty}\leq\delta\mid Y(s,t)=y(s,t)\right) \geq \mathbb{P}(A\mid Y(s,t)=y(s,t)) > 0.\] ◻
Let \(S\subset\mathbb{R}^2\) be bounded, and let \(Y(s,t)=y(s,t)\) be a fixed realisation of the latent Gaussian field. Let \(\rho>0\) be a fixed recovery time. A spatio-temporal point process \(N\) on \(S\times[0,\infty)\) with explosion time \(T_{\infty}\in(0,\infty]\) is called a log-Gaussian Cox-Hawkes process with recovery time \(\rho\) if its conditional intensity is \[\lambda_{\rho}^y(s,t\mid\mathcal{H}_t) = \exp\{\mu(s,t)+y(s,t)+D_{\rho}(s,t\mid\mathcal{H}_t)\} \mathbf{1}\{t<T_{\infty}\},\] where \[D_{\rho}(s,t\mid\mathcal{H}_t) = \alpha \sum_{i:t_i<t-\rho} \exp\{-\beta(t-t_i)\} \varphi_{\Sigma}(s-s_i).\] Here \(\rho\) represents a recovery time: an event at time \(t_i\) has no effect on the intensity during the interval \((t_i,t_i+\rho]\), and it can affect the future intensity only after time \(t_i+\rho\).
Let \(S\subset\mathbb{R}^2\) be bounded, and let \(y(s,t)\) be a fixed realisation of the latent Gaussian field. Suppose that \(b(s,t)=\mu(s,t)+y(s,t)\) is Borel measurable and locally bounded on \(S\times[0,\infty)\). Let \(\rho>0\) be fixed, and define \[D_{\rho}(s,t\mid\mathcal{H}_t)=\alpha\sum_{i:t_i<t-\rho}\exp\{-\beta(t-t_i)\}\varphi_{\Sigma}(s-s_i).\] Assume that the initial history contribution \(D_{\rho,0}(s,t)\) is finite and locally bounded on \(S\times[0,\infty)\). Then, conditional on \(Y(s,t)=y(s,t)\), the Poisson embedding construction defines a log-Gaussian Cox-Hawkes process with recovery time up to its explosion time \(T_{\infty}\). For \(t<T_{\infty}\), its conditional intensity is \[\lambda_{\rho}^y(s,t\mid\mathcal{H}_t)=\exp\{\mu(s,t)+y(s,t)+D_{\rho}(s,t\mid\mathcal{H}_t)\}.\]
Proof. This follows from the same Poisson embedding construction as Lemma 4.3, with \(D\) replaced by \(D_{\rho}\). Indeed, before explosion, the history contains only finitely many future events. Since \(\varphi_{\Sigma}\) is bounded and the initial contribution \(D_{\rho,0}(s,t)\) is locally bounded, the intensity is locally bounded on every finite time interval before the next accepted event. Hence the recursive construction is locally well defined. The resulting process satisfies the Poisson embedding representation with predictable acceptance function \(\lambda_{\rho}^y(s,t\mid\mathcal{H}_t)\). By Lemma 4.2, this accepted process has conditional intensity \(\lambda_{\rho}^y(s,t\mid\mathcal{H}_t)\) for \(t<T_{\infty}\). ◻
Under the assumptions above, the log-Gaussian Cox-Hawkes process with recovery time \(\rho>0\) is non-explosive for any \(\alpha\in\mathbb{R}\). That is, \(T_{\infty}=\infty\) a.s., conditional on \(Y(s,t)=y(s,t)\).
Proof. Events cannot affect the intensity during the first \(\rho\) units of time after they occur. Fix an interval \(I_j=(j\rho,(j+1)\rho],\: j=0,1,2,\ldots.\) For \(t\in I_j\), any event occurring after time \(j\rho\) satisfies \(t_i>j\rho\geq t-\rho\), and therefore it is not included in the sum defining \(D_{\rho}(s,t\mid\mathcal{H}_t)\). Hence events born inside \(I_j\) cannot increase the intensity inside the same interval.
We prove by induction that only finitely many future events occur in each interval \(I_j\). For \(j=0\), the intensity on \(I_0=(0,\rho]\) depends only on the initial history through \(D_{\rho,0}(s,t)\). Since \(b(s,t)\) and \(D_{\rho,0}(s,t)\) are locally bounded on \(S\times[0,\rho]\), there exists a finite constant \(M_0\) such that \[\lambda_{\rho}^y(s,t\mid\mathcal{H}_t)\leq \exp\{M_0\},\: (s,t)\in S\times I_0.\] Since \(S\) is bounded, \[\int_{I_0}\int_S \lambda_{\rho}^y(s,t\mid\mathcal{H}_t)\,ds\,dt\leq |S|\rho\exp\{M_0\}<\infty.\] Thus only finitely many events occur in \(S\times I_0\) a.s..
Now suppose that only finitely many events have occurred in \(S\times(0,j\rho]\). For \(t\in I_j\), the recovery time condition implies that \(D_{\rho}(s,t\mid\mathcal{H}_t)\) can depend only on the initial history and on the finitely many future events before time \(j\rho\). Since \(\varphi_{\Sigma}\) is bounded and \(D_{\rho,0}(s,t)\) is locally bounded, there exists a finite constant \(M_j\) such that \[\lambda_{\rho}^y(s,t\mid\mathcal{H}_t)\leq \exp\{M_j\},\: (s,t)\in S\times I_j.\] Therefore, \[\int_{I_j}\int_S \lambda_{\rho}^y(s,t\mid\mathcal{H}_t)\,ds\,dt\leq |S|\rho\exp\{M_j\}<\infty.\] Hence only finitely many events occur in \(S\times I_j\) a.s.. By induction, each interval \(I_j\) contains only finitely many future events. Now fix any finite time \(T>0\). Choose \(J\) such that \(T\leq J\rho\). Then \[N(S\times(0,T])\leq \sum_{j=0}^{J-1}N(S\times I_j)<\infty\] a.s., the process has only finitely many events on every finite time interval. Therefore, \[\mathbb{P}\left(T_{\infty}<\infty\mid Y(s,t)=y(s,t)\right)=0,\] and so \(T_{\infty}=\infty\) a.s., conditional on \(Y(s,t)=y(s,t)\). ◻
Motivated by stochastic self-propagating star formation models and recent connections between SSPSF and spatio-temporal point processes, we propose a new spatio-temporal log-Gaussian Cox-Hawkes process for stochastic star formation. Recall that the log-Gaussian Cox-Hawkes process with recovery time provides a continuous stochastic formulation of self-propagating star formation. The model is specified through the conditional intensity \[\lambda_\rho^y(s,t\mid \mathcal{H}_t) = \exp \{\mu(s,t)+y(s,t)+D_\rho(s,t\mid \mathcal{H}_t)\},\] where \[D_\rho(s,t\mid \mathcal{H}_t) = \alpha \sum_{i:t_i<t-\rho} \exp\{-\beta(t-t_i)\}\varphi_\Sigma(s-s_i).\] Here, \(u=t-t_i\) is the age of a previous star formation event, \(h=s-s_i\) is the spatial displacement from that star formation event and \(\rho\) is the recovery time.
The mean trend \(\mu(s,t)\) represents the deterministic background tendency for star formation. Physically, this component may describe large-scale galactic structure across the galactic disk. It acts as the role of spontaneous star formation in the classical SSPSF framework, but allows the spontaneous component to vary continuously over space and time. This is important because the background tendency for star formation in a galaxy is expected to vary continuously across space and over time, rather than remain fixed at a constant probability.
The latent field \(y(s,t)\) represents unobserved galactic background variation. Star formation is affected by physical conditions that may not be fully observed, such as local gas density, metallicity, pressure, magnetic fields, or dynamical structure. The log-Gaussian Cox component allows these hidden environmental effects to generate additional spatio-temporal clustering beyond that explained by the deterministic background mean field. The history-dependent field \(D_\rho(s,t\mid \mathcal{H}_t)\) represents the self-propagation of star-formation activity. It describes the effect of previous star-formation events on the future log-intensity of star formation nearby in space and time. The parameter \(\alpha\) controls the sign and strength of this history effect. When \(\alpha>0\), previous events increase the future local intensity, corresponding to positive feedback or triggered star formation. When \(\alpha<0\), previous events decrease the future local intensity, corresponding to local suppression caused by feedback that temporarily reduces the availability or suitability of nearby gas for further star formation. The parameter \(\beta>0\) controls the temporal decay of the history effect in star formation. Larger values of \(\beta\) imply that the influence of a previous event decreases more quickly, while smaller values imply a longer-lasting influence. The spatial kernel \(\varphi_\Sigma(s-s_i)\) controls the spatial range and shape of the propagation effect. The covariance matrix \(\Sigma\) determines whether this influence is isotropic or anisotropic, and therefore allows the propagation of star formation to spread differently in different spatial directions.
The recovery time \(\rho>0\) should be interpreted as a delay before a previous event can affect the future intensity. An event occurring at time \(t_i\) does not contribute to the history-dependent field during the interval \((t_i,t_i+\rho]\). It can only affect the future intensity after its age exceeds \(\rho\). It represents a finite feedback delay or source-side recovery time, rather than a cell refractory period. Physically, this delay represents the time required for massive stars to evolve, for supernova feedback or for feedback from a star-forming region to become capable of influencing neighbouring gas. Mathematically, the delayed interaction prevents events from immediately increasing the intensity after they occur, thereby avoiding instantaneous self-triggering cascades and supporting the non-explosion property of the model.
This work developed a unified spatio-temporal point-process framework for stochastic star formation. First, we established a connection between the classical stochastic self-propagating star-formation model and spatio-temporal point processes by showing that, under suitable discretisation, the SSPSF update law can be represented through a separable spatio-temporal Hawkes process. This connection provides a statistical interpretation of the SSPSF mechanism and shows how spontaneous star formation and triggered star formation can be understood in terms of background intensity and history-dependent triggering.
Building on this connection and motivated by the Cox-Hawkes framework of [41], we then introduced a new spatio-temporal log-Gaussian Cox-Hawkes process with signed history dependence. The model places the deterministic mean field, latent Gaussian field, and history-dependent interaction field inside a single log-intensity. This log-scale construction allows past events to have either excitatory or inhibitory effects while keeping the conditional intensity positive. In this way, the proposed model extends beyond additive Cox-Hawkes formulations and provides a new interpretable point-process model for stochastic star formation. In particular, it can represent latent clustering, self-excitation, local inhibition, and event-driven propagation within a single continuous spatio-temporal framework.
We also introduced a recovery time version of the log-Gaussian Cox-Hawkes process as a stochastic star-formation model with a direct physical interpretation. The recovery time prevents immediate reactivation of recently active regions and reflects the delay required before a region can again contribute to subsequent star formation. This formulation is therefore not simply a continuous version of the classical SSPSF model, but a new stochastic point-process framework for understanding, simulating, modelling, and performing inference for star-formation activity in galaxies.
Several methodological and theoretical questions remain open. On the theoretical side, following [30], this work establishes existence up to explosion, non-explosion in the non-excitatory case, possible explosion in the excitatory case, and non-explosion for the recovery-time model. However, a full stability theory is not developed here. The interaction between the latent Gaussian field and the history-dependent component may also raise further questions about parameter identifiability and asymptotic behaviour. On the methodological side, the conditional log-likelihood contains a spatio-temporal integral involving both the latent Gaussian field and the history-dependent component, which is generally not available in closed form. As a result, inference is likely to require numerical approximation together with treatment of the latent field. Possible approaches include discretisation-based likelihood approximation, Monte Carlo likelihood methods [43]–[46], and Bayesian inference [47]–[51]. Future work will focus on stability theory, parameter identifiability, estimation methods, simulation studies, and applications to simulated or observed galaxy data.
The author declares no conflict of interest.
In this section, we build a bridge between the STPP and SSPSF frameworks, organised into two parts. The notation and formulation are self-contained for each part. First, we show the equivalence between the spatio-temporal Hawkes process and the SSPSF in their update rule. Second, we derive a log-likelihood for a simplified STPP that can be used to estimate SSPSF parameters.
In this work, we focus mainly on the classical Stochastic Self-Propagating Star Formation (SSPSF) model introduced by [2]. This model represents a galactic disk as a two-dimensional array of concentric rings, each subdivided azimuthally into equal-area cells, which represent regions of the galaxy where stars can form. The disk rotates differentially according to a prescribed rotation curve, such that the set of nearest neighbours of a given cell changes over time as the rings shear past one another. Each cell of the disk can contain at most one star at a given time step, where a star refers to a star-forming region such as a young stellar cluster that is capable of producing massive stars that can trigger further star formation in surrounding regions. A supernova denotes a cell whose massive stars produce feedback strong enough to initiate star formation in its surrounding cells.
This discrete, probabilistic stochastic star formation framework couples local star formation propagation to the global galaxy disk, generating persistent spiral structure without invoking any density wave. The model evolves according to the following rules:
Initialisation of the model: randomly assign \(K\) cells with stars before simulation. Every star becomes a supernova in the very next time step after it is formed.
Triggered star: If a cell has at least one surrounding supernova cell, it can form a star with a triggered probability (\(0 \leq P_{t} \leq 1\)) , unless it is within a recovery period \(t_{r}\) following recent star formation. \[\begin{align} P_t = \begin{cases} P_{t} \dfrac{ t_{a}}{t_r}, & 0 \le t_{a} < t_r, \\ P_{t}, & t_{a} \ge t_r, \end{cases} \label{eq:refractory} \end{align}\tag{1}\] where \(t_{a}\) is the age in time steps since the last star formed in the cell.
Spontaneous star: Any cell can spontaneously form a star with a spontaneous probability (\(0 \leq P_{s} \leq 1\)), independent of its neighbours.
Aging and Triggering: the region ages by one time step at each iteration after a star was born. Massive stars trigger surrounding cells in the time step immediately after their creation.
Rotation: After each time step, rings are rotated according to the pre-designed galaxy’s rotation curve.
Define \(C_{i,j}(t)\) as the state of the cell in the ring \(i\) and the azimuth \(j\) at time \(t\) (\(0\) = empty, \(1\) = star). If \(C_{i,j}(t) = 0\), the probability of a new star forming in this region is \[\Pr(C_{i,j}(t+1) = 1) = 1 - (1 - P_{s}) \prod_{n \in N_{i,j}(t)} ( 1 - P_{t}(A_{i,j}(t))),\] where \(N_{i,j}(t)\) is the set of nearest neighbours of \(C_{i,j}(t)\) and \(A_{i,j}(t)\) = \(t_a\) time after the star born in \(C_{i,j}(t)\) at time \(t\) (0 if empty). If \(C_{i,j}(t) = 1\), the region simply ages: \(A_{i,j}(t+1) = A_{i,j}(t) + 1\) for each time step. For \(A_{i,j}(t) \ge t_r\), reset it to 0 where \(t_r\) is the refractory period in time steps. The probability of triggered star formation denoted \(P_{t}\) and the probability of spontaneous star formation \(P_{s}\) are two main parameters used to control the star formation process.
A spatio-temporal point process (STPP) on \(S\times(a,b]\) is a stochastic framework for modelling events indexed by spatial location and time. Given the event history \(\mathcal{H}_t\), the conditional intensity function of an STPP is defined informally by \[\lambda^{*}(t,s) = \lim_{\Delta t\rightarrow 0, |ds|\rightarrow 0} \frac{ \mathbb{E}\{N(ds\times[t,t+\Delta t))\mid\mathcal{H}_t\} }{ |ds|\Delta t }.\] The conditional intensity \(\lambda^{*}(t,s)\) plays a central role in spatio-temporal point processes, as it describes the instantaneous rate at which events occur near location \(s\) and time \(t\), conditional on the past history.
The Hawkes process (see e.g. [27], [28]) is a self-exciting point process in which each past event increases the conditional intensity of future events around it. In the spatio-temporal setting, events occur at times \(t \in \mathbb{R}^{+}\) and locations \(\mathbf{s} \in \mathbb{R}^{d}\) (e.g. \(d=2\) for a two-dimensional plane), and the process is specified by its conditional intensity function (e.g. [13], [27]): \[\begin{align} \lambda^*(t,\mathbf{s}) = \lambda(t,\mathbf{s} \mid \mathcal{H}_t) = \mu(\mathbf{s}) + \sum_{t_i < t} g( t - t_i, \, \mathbf{s} - \mathbf{s}_i), \end{align}\] where \(\mathcal{H}_t\) denotes the history of all events before \(t\), \(\mu(\mathbf{s})\) is a background rate that can vary across locations and be modeled by various kernels, such as a Gaussian kernel, and \(g(\Delta t, \Delta\mathbf{x})\) is the triggering kernel that describes how past events in \((t_i,\mathbf{s}_i)\) affect the rate of occurrence at \((t,\mathbf{s})\). To simplify the model, the kernel \(g\) is typically separable into temporal and spatial components: \[\begin{align} g( t - t_i, \, \mathbf{s} - \mathbf{s}_i) = g(\Delta t, \Delta\mathbf{x}) = g_1(\Delta t) \, g_2(\Delta\mathbf{x}), \\ \lambda^*(t,\mathbf{s}) = \mu(\mathbf{s}) + \sum_{t_i < t} g_1( t - t_i) g_2(\mathbf{s} - \mathbf{s}_i), \end{align}\] with \(g_1(\Delta t)\) controlling the decay of influence with time (e.g. an exponential decay function \(g_1(\Delta t) = \alpha e^{-\beta \Delta t}\) for \(\Delta t > 0\) and normally \(\alpha >0, \beta > 0\)) and \(g_2(\Delta\mathbf{x})\) describing the spatial distribution of triggered events (e.g. a Gaussian kernel).
Let events occur in continuous time \(t \in \mathbb{R}^+\) and in continuous space \(\mathbf{s}\in\mathbb{R}^2\). The spatio-temporal Hawkes process with separable kernel has conditional intensity: \[\lambda^*(t,\mathbf{s})=\mu(\mathbf{s}) +\sum_{t_i < t} g_1(t - t_i)\, g_2(\mathbf{s} - \mathbf{s}_i), \label{eq:hawkes-intensity}\tag{2}\] where \(\mu(\mathbf{s}) \ge 0\) is the background rate and \(g_1(t-t_i)\ge 0\), \(g_2(\mathbf{s} - \mathbf{s_i})\ge 0\) are the temporal and spatial trigger kernels, respectively. We also assume \(g_1(t - t_i)=0\) for \(t-t_i \le 0\).
Discretize space into equal-area cells \(\{S_{i,j}\}\) (ring \(i\), azimuth \(j\)) and time into steps of length \(\Delta \tau\), and denote discrete times as \(t_k = k\,\Delta\tau\). Let historical events \(\mathcal{H}_{t} = \{(t_i,\mathbf{s}_i): t_i < t\}\) and the event count \(N(S)\) in a spatial region \(\mathcal{S}\subset\mathbb{R}^+\times\mathbb{R}^2\). For the spatio-temporal Hawkes process, by the defining property of conditional intensities for point processes conditional on history \(\mathcal{H}_{t_k}\), the integrated intensity can be defined as: \[\begin{align} \Lambda_{i,j}(t_k) = \int_{t_k}^{t_{k+1}} \int_{S_{i,j}} \lambda^*(t,\mathbf{s}) d\mathbf{s}dt=\Lambda^{(1)}_{i,j}(t_k)+ \Lambda^{(2)}_{i,j}(t_k), \end{align}\] where \[\begin{align} \Lambda^{(1)}_{i,j}(t_k) &= \int_{t_k}^{t_{k+1}} \int_{S_{i,j}} \mu(\mathbf{s}) d\mathbf{s} dt \\ \Lambda^{(2)}_{i,j}(t_k) &= \int_{t_k}^{t_{k+1}} \int_{S_{i,j}} \sum_{t_i < t} g_1(t-t_i) g_2(\mathbf{s}-\mathbf{s}_i) d\mathbf{s} dt. \end{align}\] Considering the number of events in a cell \(S_{i,j} \times [t_k,t_{k+1})\) conditioned on \(\mathcal{H}_{t_k}\), we have \(\Pr(N(S_{i,j} \times [t_k,t_{k+1})) = 0 | \mathcal{H}_{t_k}) = \exp(-\Lambda_{i,j}(t_k)),\) since the count of events in the cell is conditionally Poisson with mean equal to the integrated intensity defined above. \[\begin{align} \Pr(N(S_{i,j} \times [t_k,t_{k+1}))\ge 1 | \mathcal{H}_{t_k}) = 1 -\exp(-\Lambda_{i,j}(t_k)). \label{eq:atleastone} \end{align}\tag{3}\] Then decompose \(\Lambda_{i,j}(t_k)\) as a sum of background and trigger contributions. For each neighbour cell \(S_{m,n}\in\mathcal{N}_{i,j}(t_k)\) that contains a supernova at some time \(t_l \le t_k\), with \(m,n\) denoting the neighbour indices: \[\Lambda_{i,j}(t_k) = \Lambda^{(1)}_{i,j}(t_k) + \sum_{S_{m,n}\in\mathcal{N}_{i,j}(t_k)} \Lambda^{(2)}_{i,j \leftarrow m,n}(t_k). \label{eq:Lambdasum}\tag{4}\] Using 3 and 4 , and assuming each cell can have at most one star, we have the following: \[\begin{align} \Pr(C_{i,j}(t_{k+1})=1 | C_{i,j}(t_k)=0, \mathcal{H}_{t_k}) &= 1 - \exp(-\Lambda^{(1)}_{i,j}(t_k)) \exp(- \sum_{S_{m,n}\in \mathcal{N}_{i,j}(t_k)} \Lambda^{(\mathrm{2})}_{i,j \leftarrow m,n}(t_k))\\ &= 1 - \exp(-\Lambda^{(1)}_{i,j}(t_k)) \prod_{S_{m,n}\in\mathcal{N}_{i,j}(t_k)} \exp(-\Lambda^{(2)}_{i,j\leftarrow m,n}(t_k)). \label{eq:sthp} \end{align}\tag{5}\] Define spontaneous and triggered probabilities by \[\begin{align} P_s(i,j,t_k) = 1 - \exp(-\Lambda^{(1)}_{i,j}(t_k)), \: \exp(-\Lambda^{(1)}_{i,j}(t_k)) = 1 - P_s(i,j,t_k) \end{align}\] \[\begin{align} P_t(i,j\leftarrow m,n;t_k) = 1 - \exp(-\Lambda^{(\mathrm{2})}_{i,j\leftarrow m,n}(t_k)), \: \exp(-\Lambda^{(\mathrm{2})}_{i,j\leftarrow m,n}(t_k)) = 1 - P_t(i,j\leftarrow m,n;t_k). \end{align}\] Then, the probability of the spatio-temporal Hawkes process becomes \[\begin{align} \Pr(C_{i,j}(t_{k+1})=1 | C_{i,j}(t_k)=0, \mathcal{H}_{t_k}) = 1 - (1 - P_s(i,j,t_k)) \prod_{S_{m,n}\in\mathcal{N}_{i,j}(t_k)} (1 - P_t(i,j\leftarrow m,n;t_k)), \end{align}\] which is exactly the SSPSF probability for updating stars. If \(\Lambda^{(1)}_{i,j}\) and \(\Lambda^{(2)}_{i,j}\) are treated as constant within cells, and the temporal kernel is multiplied by a function, e.g. \(r(t_a,t_r) = \min(1, \frac{t_a}{t_r})\) to introduce a refractory ramp, and define the refractory kernel as \(\tilde{g}_1(\Delta t) = g_1(\Delta t)\, r(t_a, t_r).\) Replacing \(g_1\) with \(\tilde{g}_1\) makes triggering impossible immediately after an event, then allows it to increase linearly up to full strength by time \(t_r\), matching the refractory behavior of the SSPSF. Then, by adding a differential rotation map to the spatio-temporal Hawkes process when generating offspring events in the Stochastic Self-Propagating Star Formation (SSPSF) model, we obtain an equivalent formulation. We show that the one-step update probability of the SSPSF model is the probability induced by a separable spatio-temporal Hawkes process on the corresponding space-time cell, when the Hawkes intensity is integrated over that cell and one discrete time step, with both refractory modulation and differential rotation. Moreover, this formulation offers a more flexible framework for modelling the SSPSF, since the kernel structures can be specified with greater generality.
Let events \(\mathcal{H}_{t} = \{(t_n,\mathbf{s}_n)_{n=1}^N\}\) be observed on \([0,T]\times S\), with separable spatio-temporal intensity \(\lambda^*(t,\mathbf{s})\). The continuous log-likelihood of the spatio-temporal point processes (see e.g., [14], [27]) is \[L =\sum_{n=1}^N \log \lambda^*(t_n,\mathbf{s}_n) -\int_0^T\!\!\int_S \lambda^*(t,\mathbf{s})\,d\mathbf{s}\,dt.\] To pass from continuous space and time to discrete space and time, let \(\{B_{i,j,k}\}\) be a measurable partition of \([0,T]\times S\) into disjoint time space cells where \((i,j)\) are spatial indices and \(k\) is time index, then \(B_{i,j,k}=S_{i,j}\times[t_k,t_{k+1}), t_k=k\,\Delta\tau\) and let \(\Delta\tau = 1\) in our case. For each cell define the integrated intensity, \[\Lambda_{i,j}(t_k)=\int_{B_{i,j,k}}\lambda^*(t,\mathbf{s})d\mathbf{s}dt,\] and let \(N_{i,j,k}=N(B_{i,j,k})\) be the observed number of events in a space-time cell \(B_{i,j,k}\). Regrouping the event term by cells and decomposing the integral over the partition for the continuous log-likelihood, we have \[\begin{align} L &=\sum_{n=1}^N \log \lambda^*(t_n,\mathbf{s}_n) -\int_0^T\int_S \lambda^*(t,\mathbf{s})d\mathbf{s}dt\\ & = \sum_{i,j,k}\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log \lambda^*(t_n,\mathbf{s}_n) - \sum_{i,j,k} \int_{B_{i,j,k}} \lambda^*(t,\mathbf{s})d\mathbf{s}dt \\ & = \sum_{i,j,k}\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log \lambda^*(t_n,\mathbf{s}_n) - \sum_{i,j,k} \Lambda_{i,j}(t_k) \end{align}\] For any event, \((t_n,\mathbf{s}_n)\in B_{i,j,k}\), \[\log \lambda^*(t_n,\mathbf{s}_n) =\log(\Lambda_{i,j}(t_k)\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)})=\log \Lambda_{i,j}(t_k) +\log\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)}.\] Then, summing all events in the cell \(B_{i,j,k}\), \[\begin{align} \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log \lambda^*(t_n,\mathbf{s}_n) &=\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}}\log \Lambda_{i,j}(t_k)+\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)} \\ &=N_{i,j,k}\log \Lambda_{i,j}(t_k) +\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)} \end{align}\] Assume \(\log \Lambda_{i,j}(t_k)\) is constant within the cell and the number of events in the cell is \(N_{i,j,k}\). The log-likelihood becomes \[\begin{align} L &= \sum_{i,j,k}\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log \lambda^*(t_n,\mathbf{s}_n) - \sum_{i,j,k} \Lambda_{i,j}(t_k)\\ &= \sum_{i,j,k} N_{i,j,k}\log \Lambda_{i,j}(t_k) + \sum_{i,j,k}\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)} - \sum_{i,j,k} \Lambda_{i,j}(t_k)\\ & =\sum_{i,j,k}( N_{i,j,k}\log\Lambda_{i,j}(t_k) - \Lambda_{i,j}(t_k)) +\sum_{i,j,k}\sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)} \end{align}\] Define the parameter \(p = 1 - \exp(-\Lambda_{i,j}(t_k))\), then \(\Lambda_{i,j}(t_k) = -\log(1-p)\), \[\begin{align} N_{i,j,k}\log\Lambda_{i,j}(t_k) - \Lambda_{i,j}(t_k) &= N_{i,j,k}\log(-\log(1-p)) +\log(1-p). \end{align}\] Now, assume \(\Lambda_{i,j}(t_k)\) is small and approximately constant within a cell, and hence \(\lambda^*(t_n,s_n) \approx \Lambda_{i,j}(t_k)\) in a cell. Then, \(p = 1 - \exp(-\Lambda_{i,j}(t_k))\) can be treated as approximately constant and since \(\exp(x)\approx 1+x\) for small \(x\), then \(\Lambda_{i,j}(t_k) \approx \exp(\Lambda_{i,j}(t_k)) - 1\). We have the following approximation of \(\lambda^*(t_n,s_n)\), \[\begin{align} \lambda^*(t_n,s_n) &\approx \exp(\Lambda_{i,j}(t_k)) - 1 \\ &=\frac{1}{\exp(-\Lambda_{i,j}(t_k)) } - \frac{\exp(-\Lambda_{i,j}(t_k)) }{\exp(-\Lambda_{i,j}(t_k)) } \\ &= \frac{1-\exp(-\Lambda_{i,j}(t_k)) }{\exp(-\Lambda_{i,j}(t_k)) } =\frac{p}{1-p}. \end{align}\] Substituting the approximation of \(\lambda^*(t_n,s_n)\) and based on our assumptions, we have \[\begin{align} \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log\frac{\lambda^*(t_n,\mathbf{s}_n)}{\Lambda_{i,j}(t_k)} & = \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log(\lambda^*(t_n,\mathbf{s}_n)) - \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log(\Lambda_{i,j}(t_k))\\ & = \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log(\lambda^*(t_n,\mathbf{s}_n)) - \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log(-\log(1-p))\\ & = \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log(\frac{p}{1-p}) - \sum_{(t_n,\mathbf{s}_n)\in B_{i,j,k}} \log(-\log(1-p))\\ & = N_{i,j,k} \log(\frac{p}{1-p}) - N_{i,j,k} \log(-\log(1-p)). \end{align}\] Putting all together, we obtain the following log-likelihood function for the simplified spatio-temporal point process of the reduced stochastic self-propagating star formation: \[\begin{align} L &= \sum_{i,j,k}\{N_{i,j,k}\log(-\log(1-p)) +\log(1-p) + N_{i,j,k} \log(\frac{p}{1-p}) - N_{i,j,k} \log(-\log(1-p))\}\\ &= \sum_{i,j,k}\{\log(1-p) + N_{i,j,k} log(\frac{p}{1-p})\} \\ &= \sum_{i,j,k}\{\log(1-p) + N_{i,j,k} log(p) - N_{i,j,k} log(1-p)\}\\ &= \sum_{i,j,k}\{N_{i,j,k} \log(p) + (1- N_{i,j,k}) \log(1-p)\} \end{align}\] which is the log-likelihood function of a Bernoulli distribution ([52], [53], etc.) with the parameter \(p=1-\exp(-\Lambda_{i,j}(t_k))\) and \(\Lambda_{i,j}(t_k)\) should be specified according to the STPP model being used, e.g. the spatio-temporal Hawkes process (STHP). This simple log-likelihood can be used to estimate parameters from star formation data. More generally, the spatiotemporal likelihood provides a flexible way to adapt other SSPSF variants by modifying or adding components, and it provides a path toward a fully continuous SSPSF formulation within a spatio-temporal point-process framework.