A Foundation for Gravitational-Wave Population Inference within the LISA Global Fit

Abstract: Population inference in gravitational-wave astronomy allows us to connect individual detections to the astrophysics of compact objects and their environments. Current approaches employed for population inference with LIGO-Virgo-KAGRA data approximate evaluation of the hierarchical population likelihood via post-processing of individual-event posteriors. However, the case of the Laser Interferometer Space Antenna (LISA) will be more complex for two main reasons: the transdimensional "global fit" approach to LISA data analysis which models all signals and noise simultaneously, and the presence of both individually-resolved signals and the unresolved stochastic ``Galactic foreground" arising from the Galactic binary population, which induces a circular dependence between the resolved and unresolved systems and our ability to detect the former. These challenges are not without opportunity; LISA's data will contain every mHz compact binary in the Milky Way -- either individually or within the Galactic foreground -- with great potential for Galactic and stellar astrophysics. We therefore propose an alternative approach: direct evaluation of the full hierarchical population likelihood within the LISA global fit. We develop a statistical formalism for joint inference of individually-resolved gravitational-wave sources, an unresolved stochastic foreground, and a shared, underlying astrophysical population, present PELARGIR, a prototype GPU-accelerated population inference module for the LISA global fit, demonstrate the formalism and PELARGIR via a toy model analysis, and lay out a roadmap towards an astrophysically-motivated LISA global fit with embedded population inference. While we apply the formalism here to the population of LISA Galactic binaries, it is applicable across the gravitational-wave spectrum with use cases in pulsar timing and next-generation terrestrial observatories.

• The paper introduces a hierarchical framework for joint inference of resolved and unresolved gravitational-wave sources in LISA, resolving the circularity between source detectability and foreground modeling.
• It implements an efficient GPU-accelerated module, PELARGIR, to partition millions of Galactic binaries in near real-time, significantly improving computational performance.
• A toy model validates the method by accurately recovering astrophysical hyperparameters, supporting its application to complex future LISA analyses.

Framework for Joint Population Inference in the LISA Global Fit

The paper "A Foundation for Gravitational-Wave Population Inference within the LISA Global Fit" (2604.03390) addresses the structural and statistical challenges inherent to population inference in the context of the Laser Interferometer Space Antenna (LISA). It introduces a general, fully hierarchical formalism for jointly inferring parameters describing both resolved and unresolved gravitational-wave (GW) sources—specifically Galactic binaries (GBs)—together with the properties of their underlying astrophysical population, while concurrently accounting for instrumental and astrophysical noise in a transdimensional analysis context.

---

Motivation and Problem Context

Current hierarchical population inference for ground-based GW detectors proceeds via a two-step method: single-event inference is performed on resolved sources, then population parameters are inferred in a subsequent post-processing phase. However, LISA differs fundamentally owing to two key aspects:

• Global Fit and Transdimensionality: The LISA data stream is highly signal-rich, requiring the simultaneous modeling of all GWs and noise components via a blocked Gibbs framework and reversible-jump Markov Chain Monte Carlo (RJMCMC) for variable event dimensionality.
• Resolved vs. Unresolved Signal Circularity: Tens of millions of compact GBs exist in the Milky Way, but only $\sim$0.1% are individually resolvable. The rest form a stochastic "foreground" that is not separable from instrumental noise, and the ability to resolve individuals is itself contingent on the current foreground level, resulting in unavoidable circular dependence between population description, source detection, and the foreground.

The distinction from the ground-based, event-limited regime is nontrivial: in LISA, the target population is frequency-complete, removing explicit selection effects and altering the mathematical structure of the hierarchical Bayesian inference problem. Prior, in-post inference approaches (including importance sampling over event catalogs) are shown to be inefficient and potentially bias-prone in the high-dimensional, transdimensional LISA setting.

---

Hierarchical Statistical Formalism

The paper develops a formal statistical framework wherein the full hierarchical likelihood is sampled directly within the LISA global fit, using a semi-analytic thresholding model to resolve the circularity between population parameters, source detectability, and the unresolved foreground power spectral density (PSD).

Let $\Lambda$ denote population-level hyperparameters, $N$ the total number of binaries, and $\eta$ describe the instrumental noise model. For the resolved subset, parameters $\{\vec{\theta}_i\}$ specify source properties. The remaining binaries form the unresolved stochastic foreground, characterized by its frequency and sky-dependent PSD, $S(f, \hat{\Omega} | N, \Lambda)$. The formalism constructs the data likelihood via a frequency-domain (Whittle) multivariate Gaussian with covariance

$$C(f; t | N, \Lambda, \eta) = C_n(f; t | \eta) + C_\mathrm{GW}(f; t | N, \Lambda) \,,$$

where $C_\mathrm{GW}$ integrates the astrophysical foreground for the given population.

Population inference proceeds by jointly sampling $(N, \Lambda, \{\vec{\theta}_i\}, S, \eta)$ using hierarchical priors that encode all circular dependencies between source classification (resolved/unresolved), population generation, and the induced foreground. Crucially, no explicit selection term is included, as the partition into resolved and unresolved systems is determined self-consistently for a given population realization and SNR threshold.

---

Efficient GPU-Accelerated Implementation: PELARGIR

A key computational challenge is rapidly partitioning a $\sim 10^7$-sized synthetic GB population at each likelihood evaluation. The authors introduce PELARGIR, a GPU-accelerated inference module implementing an optimized, parallelized thresholding algorithm that sorts binaries by instrumental SNR and iteratively computes which cross the resolution boundary as the composite noise (instrumental plus foreground) is updated.

Figure 1: Result of the thresholding algorithm as applied to the fiducial population synthesis catalogue of Thiele et al. (2023), efficiently separating resolved and unresolved Galactic binaries.

Evaluation at a frequency resolution $\Lambda$0 Hz for a million-source catalog takes 1–6 seconds, enabling practical integration within a high-dimensional global fit.

---

Numerical Demonstration via a Toy Model

To validate the formalism, the authors construct a toy model with a simplified GB population (white dwarf binaries), containing both resolved and unresolved sources. The foreground is modeled as a 1D PSD, neglecting anisotropy for computational tractability. They sample over six astrophysical hyperparameters (mass, separation, disk/bulge structure) and demonstrate that the joint inference formalism can recover true simulated population parameters, the resolved/unresolved source distributions, and the total foreground spectrum, even in the presence of complex population-induced spectral structure.

Figure 2: Posterior distributions of population hyperparameters from the toy model analysis, demonstrating accurate recovery to within the quoted credible intervals.

Figure 3: Comparison of recovered population distributions (credible intervals) with histograms of the true full and resolved GB populations, highlighting the differing mass, separation, and distance characteristics between resolved and unresolved systems.

The method self-consistently captures the number of resolved sources, population-level priors, and jointly infers the spectral shape of the unresolved foreground.

Figure 4: Population-derived prior on the number of resolved binaries, showing consistency with the underlying simulated value.

Figure 5: Posterior on the total PSD (instrumental plus foreground), showing strong informational gain for the astrophysical component without undue bias.

---

Theoretical and Practical Implications

A significant theoretical advance is the elimination of inefficiencies and potential biases intrinsic to two-step or in-post importance sampling approaches—particularly as event numbers and overall parameter-space dimensionality scale up in LISA. The formalism elegantly encapsulates the unavoidable circularity between source detectability and foreground properties, crucial for both astrophysical and cosmological population inference. Practically, the population-inferred priors on resolved GBs can be propagated to source-level parameter estimation, improving both search sensitivity and confidence in source classification, while the same framework provides an objective, data-driven model of the complex Galactic foreground for subsequent LISA analyses.

Additionally, with minor adaptation, this framework is extensible to multiple, coexisting source classes (e.g., EMRIs, GW backgrounds from dwarf galaxies), and analogous structure is directly applicable to hierarchically informed population inference in the PTA SMBH binary context and in next-generation ground-based GW detectors.

---

Prospects for Future Developments

Key directions for further development include:

• Astrophysical Realism & Complexity: Incorporation of anisotropy, realistic waveform models, and heterogeneous populations (white dwarf, NS, BH binaries) will be required for the full LISA data analysis application.
• Accelerated Inference via Flow Emulation: The authors argue for the use of machine learning surrogates, particularly normalizing flows trained on thresholded populations, to directly emulate the high-dimensional mapping between hyperparameters and observables, further reducing computational demands and allowing for joint inference in even more complicated parameter spaces.
• Cross-Domain Applications: As LISA, PTAs, and third-generation terrestrial detectors approach the confusion-limited regime, this direct in-situ hierarchical approach will be the standard for unbiased, physically meaningful population inference.

---

Conclusion

This work provides the first general, practical statistical and computational foundation for fully joint inference of resolved and unresolved GW source populations within the LISA global fit. By directly embedding the hierarchical population model, the approach naturally resolves the intricate circularity between source detection and astrophysical foregrounds, respects the transdimensional nature of the LISA data, and is efficiently realizable with modern hardware acceleration. This framework will underpin both astrophysical characterizations of compact binary populations and robust subtraction and modeling of the LISA Galactic foreground, with immediate relevance for Galactic structure, stellar evolution, and the search for cosmological GW backgrounds.