Here we give a more abstract explanation of how one can construct a generative model within the X-PSI framework, and apply that model to astronomical X-ray photon data.

The overarching concept

X-PSI is designed for the statistical analysis of astronomical X-ray photon data sets acquired through observation of pulsed emission from neutron stars. To achieve this, X-PSI couples code for likelihood function evaluation with existing open-source software for posterior sampling.

X-ray photon data and likelihood functions

The relevant data is acquired via direct astronomical detection of photons in the X-ray regime.

The likelihood function has domain \(\mathbb{R}^{n}\), defined at each point \(\boldsymbol{\theta}\) (optionally with finite local prior probability density) by a parametrised sampling distribution on the space of the data evaluated at the fixed dataset \(\mathcal{D}\). The likelihood is defined by


where \(\mathcal{M}\) is the global (hierarchical) model.

A fundamental assumption upon which X-PSI is based (at least until posterior checking suggests otherwise) is that a model with sufficient complexity to describe data can be constructed from a restricted set of abstract objects that form a modelling language. Crucially, those objects must have software implementations that are sufficiently fast for parameter estimation to be tractable on a high-performance computing system, but are also simple for a modeller to piece together into a bespoke model. These objects often have properties which are discrete representations of continuous mathematical structures, and these continuous structures are imbued with physical meaning in order to attempt to describe observable reality in terms of deterministic and stochastic processes.

In X-PSI, we adopt general relativistic gravity to construct a spacetime manifold associated with a model neutron star. The exterior spacetime of the star is stationary (and thus time-independent) during the acquisition of all data \(\mathcal{D}\).[1] Further, the exterior spacetime is assumed to be non-rotating for the purpose of constructing models that are expected to be more tractable given some resource allowance than models that account for rotation in the exterior metric (this is reasonable for stars rotating at up to a few hundred Hz, see the discussion in Bogdanov et al. 2019). Spherical symmetry is arguably a logical baseline assumption to condition on at the outset of modelling some new data set.

We first need to consider our X-ray photon data set \(\mathcal{D}\): in general let it be constituted by some union of data subsets \(\mathcal{D}=\cup_{i}\mathcal{D}_{i}\), such that \(\mathcal{D}_{i}\cap\mathcal{D}_{j}=\emptyset\) for \(i\neq j\), and the subsets are statistically independent—i.e., the joint sampling distribution conditional on the model is separable. We describe the stochastic generation of this data via observations—with a model X-ray telescope—of a single model star, whose exterior spacetime parameters we are interested in constraining (primarily for the purpose of dense matter study). Further, we assume that the data subsets \(\mathcal{D}_{i}\) are realistically acquired during mutually disjoint time intervals.

Let us define a data subset \(\mathcal{D}_{i}\) as having a parametrised model sampling distribution that is written in terms of a single pulse generated by a rotating non-axisymmetric surface radiation field. Further, let the radiation field be everywhere invariant in surface local comoving frames (which need not necessarily be photospheric rest frames). All photon events constituting \(\mathcal{D}_{i}\) need not be acquired within a single period of rotation of the star—indeed, in reality the stars of interest are of galactic origin but are sufficiently distant that the incident flux of photons is very small.

The photon events which comprise \(\mathcal{D}_{i}\) can span an arbitrarily long time interval, but over that time interval the surface radiation field is assumed to be stable, with any quasi-periodicity in the signal incident on a telescope manifesting due to relative motion the star and the telescope, and/or due to a small evolution in the coordinate rotation frequency of a mode of asymmetry in the surface radiation field. Each subset \(\mathcal{D}_{i}\) is acquired with some model instrument with a response matrix (a photon energy redistribution matrix combined with an energy-dependent effective area vector) which is used to transform a radiation field incident on the telescope into a form which enters directly in the sampling distribution of \(\mathcal{D}_{i}\).


The joint posterior distribution on parameter space is the joint probability density distribution proportional to the integrable product of the likelihood function with a joint prior distribution:


The normalisation is the prior predictive distribution on the space of the data evaluated at \(\mathcal{D}\):

\[\mathcal{P}(\mathcal{D}\;|\;\mathcal{M},\mathcal{I}) =\mathop{\int}\mathcal{P}(\mathcal{D}\;|\;\boldsymbol{\theta},\mathcal{M})\mathcal{P}(\boldsymbol{\theta}\;|\;\mathcal{M},\mathcal{I})d^{n}\boldsymbol{\theta}.\]

Also termed the evidence or fully marginal likelihood, this normalisation can be approximated by a nested sampler such as MultiNest. However, the interpretation and robustness of the evidence in the context of model comparison is problem-dependent a subject of much debate in the literature.

Parallelisation paradigms

X-PSI inherits the MPI parallelisation of an external sampling package. In general it is necessary to run X-PSI posterior sampling/integration on distributed memory architectures (e.g., a cluster or a supercomputer) because likelihood evaluation times are slow, of \(\mathcal{O}(1)\,s\).

In addition, the source code for pulse simulation (required for likelihood evalution) is OpenMP-enabled, meaning that in principle one can explore hybrid parallelisation paradigms—in particular, enabling multithreaded computation on shared memory nodes of distributed architectures.

Direct inference of Equation of State

In principle, support for direct interior (Equation of State) parameter estimation could be added (see Riley et al. 2018). This would involve solving the field equations for stable stationary spacetime solutions and embedding the rotationally deformed stellar surface into an ambient Schwarzschild spacetime with the rotationally deformed gravitational mass monopole moment, whilst neglecting all other rotational \(\ell>0\) metric deformations at \(\mathcal{O}(\Omega)\) (see Hartle, 1967; Hartle & Thorne, 1968).