Flexible Spectral Separation of Multiple Isotropic and Anisotropic Stochastic Gravitational Wave Backgrounds in LISA

TL;DR

This paper introduces BLIP 2.0, a GPU-accelerated framework for separating multiple stochastic gravitational wave backgrounds in LISA data, enabling spectral analysis of diverse astrophysical and cosmological signals.

Contribution

The paper presents a novel flexible framework for spectral separation of multiple isotropic and anisotropic SGWBs, including a demonstration with LISA data.

Findings

First spectral separation of Galactic foreground, LMC SGWB, and extragalactic binaries.

Proof-of-concept for upper limits on cosmological SGWB detection.

Framework accelerates Bayesian analysis of complex SGWB combinations.

Abstract

The Laser Interferometer Space Antenna (LISA) will observe mHz gravitational waves from a wide variety of astrophysical sources. Of these, some will be characterizable as individual deterministic signals; the remainder will overlap to create astrophysical confusion noise. These sources of confusion noise are known as stochastic gravitational wave backgrounds (SGWBs). LISA data is expected to include several such astrophysical SGWBs, including the notable Galactic binary foreground, SGWBs from white dwarf binary populations in satellite galaxies of the Milky Way, and the SGWB from extragalactic stellar-origin binary black holes far from merger. To characterize these astrophysical signals and attempt to seek out possible underlying backgrounds of cosmological origin, it will be necessary to separate the contribution of each SGWB from that of the others. Crucially, several of these SGWBs are expected to be highly anisotropic on the sky, providing a powerful tool for spectral separation. To this end, we present BLIP 2.0: a flexible, GPU-accelerated framework for simulation and Bayesian analysis of arbitrary combinations of isotropic and anisotropic SGWBs. We leverage these capabilities to demonstrate for the first time spectral separation of the Galactic foreground, the Large Magellanic Cloud SGWB, and the SGWB from extragalactic stellar-origin binaries, and show a proof-of-concept for placing upper limits on the detection of an underlying isotropic cosmological SGWB in the presence of multiple astrophysical foregrounds.

Tables

Table 1. Table 1: Summary of current models available in BLIP , both spectral and spatial. For further details, see the listed reference or section; for parameter definitions refer to Table  2 . ( Note: the other parameters in the analytic foreground spectral model are fixed for a given observing duration. See III.2.2 for details.)

Name	Type	Parameters	Formula	Reference
Power Law	Spectral	${\Omega }_{\mathrm{ref}},\alpha$	${\Omega }_{\mathrm{GW}}\left(f\right)={\Omega }_{\mathrm{ref}}{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{\alpha }$	Banagiri et al. [23]
Broken Power Law	Spectral	${\Omega }_{\mathrm{ref}}$, ${\alpha }_1$, ${\alpha }_2$, $f_{\mathrm{break}}$, $\delta$	${\Omega }_{\mathrm{GW}}\left(f\right)={\Omega }_{\mathrm{ref}}{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{{\alpha }_1}{\left(1+{\left(\frac{f}{f_{\mathrm{break}}}\right)}^{\frac{1}{\delta }}\right)}^{\delta \left({\alpha }_1-{\alpha }_2\right)}$	§III.2.1
Truncated Power Law	Spectral	${\Omega }_{\mathrm{ref}}$, $\alpha$, $f_{\mathrm{cut}}$, $f_{\mathrm{scale}}$	${\Omega }_{\mathrm{GW}}\left(f\right)=\frac{1}{2}{\Omega }_{\mathrm{ref}}{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{\alpha }\left(1+\tanh \left(\frac{f_{\mathrm{cut}}-f}{f_{\mathrm{scale}}}\right)\right)$	Criswell et al. [46]
Analytic Foreground	Spectral	${\Omega }_{\mathrm{ref}}$ (see note)	${\Omega }_{\mathrm{GW}}\left(f\right)={\Omega }_{\mathrm{ref}}{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{2/3}{e}^{-f^{{\alpha }_{\mathrm{C17}}}-\beta f\sin \left(\kappa f\right)}\times \left[1+\tanh \left(\gamma \left(f_k-f\right)\right)\right]$	§III.2.2
Spherical Harmonic	Spatial	${\overrightarrow{b}}_{\ell m}$	Spherical harmonic expansion of $\sqrt{𝒫\left(𝐧\right)}$; see ref.	Banagiri et al. [23]
Analytic Galaxy	Spatial	$r_h$, $z_h$	Pixel-basis template created from a simple bulge + exponential disk density distribution; see ref.	Criswell et al. [46]
Analytic Satellite	Spatial	$d$, $r$, RA, DEC	Pixel-basis template created from a simple spherical density distribution.	§III.2.3
Point Source (Radiometer)	Spatial	${\{{\theta }_i,{\varphi }_i\}}_{i\in 1\dots N}$	Pixel-basis template created from a variable number of point sources.	§III.2.4
popmap	Spatial	—	Pixel-basis template created from any DWD population synthesis catalogue; see ref.	Criswell et al. [46]

Table 2. Table 2: Parameter definitions for the models discussed in this work.

Parameter	Definition	Reference
${\Omega }_{\mathrm{ref}}$	Amplitude at $f_{\mathrm{ref}}$, in dimensionless GW energy density.	Banagiri et al. [23]
$f_{\mathrm{ref}}$	Reference frequency. 1 mHz unless specified otherwise.	Banagiri et al. [23]
$\alpha$	Power law slope. Subscripts (e.g. ${\alpha }_1$, ${\alpha }_2$) denote multiple power laws.	Banagiri et al. [23]
$f_{\mathrm{break}}$	Break frequency for the broken power law model.	§III.2.1
$\delta$	Smoothing parameter the broken power law model.	§III.2.1
$f_{\mathrm{cut}}$	Truncation frequency for the truncated power law model.	Criswell et al. [46]
$f_{\mathrm{scale}}$	Truncation scaling parameter for the truncated power law model.	Criswell et al. [46]
${\overrightarrow{b}}_{\ell m}$	Coefficients of the spherical harmonic expansion of $\sqrt{𝒫\left(𝐧\right)}$	Banagiri et al. [23]
$r_h$	Radial scale height for the analytic Milky Way disk model of Breivik et al. [33].	Criswell et al. [46]
$z_h$	Vertical scale height for the analytic Milky Way disk model of Breivik et al. [33].	Criswell et al. [46]
$d$	Distance to the center of the analytic MW satellite.	§III.2.3
$r$	Radius of the (assumed spherical) analytic MW satellite.	§III.2.3
RA, DEC	Right ascension and declination of the analytic MW satellite.	§III.2.3
${\theta }_i$, ${\varphi }_i$	Sky location of the $i^{\mathrm{th}}$ point source, for $i\in 1\dots N$. Can also be given in RA/DEC.	§III.2.4

Table 3. Table 3: Summary of simulation components.

Component	Spectral	Parameters	Value	Spatial	Parameters	Values	Notes
Noise	Instrumental Noise	$N_a$	$3.6\times {10}^{-49}$ Hz-4	—	—	—	Per Amaro-Seoane et al. [15].
		$N_p$	$9\times {10}^{-42}$
MW		${\Omega }_{\mathrm{ref}}$	$2.197\times {10}^{-9}$				Thin disk MW model of Breivik et al. [33].
	Tanh-truncated	$\alpha$	0.667	Analytic	$r_h$	2.9 kpc
	Power Law	$f_{\mathrm{cut}}$	2.08 mHz	Galaxy	$z_h$	0.3 kpc
		$f_{\mathrm{scale}}$	0.5 mHz
LMC		${\Omega }_{\mathrm{ref}}$	$7.6\times {10}^{-12}$		$r$	2.15 kpc	Based on LMC spectrum of Rieck et al. [91].
	Broken Power Law	${\alpha }_1$	0.667	Analytic Satellite	$d$	50 kpc
		${\alpha }_2$	2.65		RA	5h23m34s
		$f_{\mathrm{break}}$	3.83 mHz		DEC	-69°45’22”
		$\delta$	0.1
SOBBH	Power Law	${\Omega }_{\mathrm{ref}}$	$2.11\times {10}^{-12}$	Isotropic	—	—	Upper range estimate of Babak et al. [21].
		$\alpha$	0.667
CGWB	Power Law	${\Omega }_{\mathrm{ref}}$	$1.0\times {10}^{-13}$	Isotropic	—	—
		$\alpha$	1.0

Table 4. Table 4: Summary of priors and model choices.

Submodel	Spectral	Parameters	Prior/Fixed Value	Spatial	Parameters	Template Values
Noise	Instrumental Noise	${\log }_{10}{N}_a$	$𝒰(-51,-46)$	—	—	—
		${\log }_{10}{N}_p$	$𝒰(-44,-39)$
MW		${\log }_{10}{\Omega }_{\mathrm{ref}}$	$𝒰(-10,-7)$
	Tanh-truncated	$\alpha$	$𝒰(0,2)$	Analytic	$r_h$	2.9 kpc
	Power Law	$f_{\mathrm{cut}}$	$𝒰(0.8,4)$ mHz	Galaxy	$z_h$	0.3 kpc
		$f_{\mathrm{scale}}$	$𝒰(0.1,10)$ mHz
LMC		${\log }_{10}{\Omega }_{\mathrm{ref}}$	$𝒰(-9,-12)$		$r$	2.15 kpc
	Broken Power Law	${\alpha }_1$	$𝒰(0,2)$	Analytic Satellite	$d$	50 kpc
		${\alpha }_2$	$𝒰(0,4)$		RA	5h23m34s
		$f_{\mathrm{break}}$	$𝒰(1,10)$ mHz		DEC	-69°45’22”
		$\delta$	0.1
SOBBH	Power Law	${\log }_{10}{\Omega }_{\mathrm{ref}}$	$𝒰(-11,-13)$	Isotropic	—	—
		$\alpha$	0.667
CGWB	Power Law	${\log }_{10}{\Omega }_{\mathrm{ref}}$	$𝒰(-21,9)$	Isotropic	—	—
		$\alpha$	$1$

Equations

$$\mathcal{L}(d|\vec{\theta })=\prod _{f,t}\frac{1}{2\pi T_{\mathrm{seg}}\det \left(C_{IJ}(f,t|\vec{\theta })\right)}\times \exp \left[-\frac{2}{T_{\mathrm{seg}}}\sum _{I,J}{\tilde{d}}_I^{\ast }(f,t){\left(C(f,t|\vec{\theta })\right)}_{IJ}^{-1}{\tilde{d}}_J(f,t)\right],$$

$$C_{IJ}(f,t)=S_{IJ}^n(f)+S_{IJ}^{\mathrm{GW}}(f,t),$$

$$S_{\mathrm{GW}}(f,\mathbf{n})=S_{\mathrm{GW}}(f)\mathcal{P}(\mathbf{n}),$$

$$\int \mathcal{P}(\mathbf{n})d^2\mathbf{n}=1.$$

$$S_{IJ}^{\mathrm{GW}}(f,t)=\int \mathcal{P}(\mathbf{n}){\mathcal{R}}_{IJ}(f,t,\mathbf{n})d^2\mathbf{n}\times S_{\mathrm{GW}}(f).$$

$$\mathcal{L}_{\mathrm{SI}}(d|\vec{\theta})=\prod_{t,f}\left[\frac{1}{2\pi T_{\mathrm{seg}}\det\left(C_{IJ}(t,f|M_{1...N};\vec{\theta})\right)}\times\exp\Biggl{[}-\frac{2}{T_{\mathrm{seg}}}\sum_{I,J}\tilde{d}_{I}^{*}(t,f)\left(C_{IJ}(t,f|M_{1...N};\vec{\theta})\right)^{-1}\tilde{d}_{J}(t,f)\Biggr{]}\right],$$

$$C_{IJ}(t,f|M_{\mathrm{1...}N};\vec{\theta })=\sum _{k=1}^N{C}_{IJ}^{M_k}(t,f|{\vec{\theta }}_k),$$

$$\begin{split}\vec{\theta}&=\bigl{\{}\vec{\theta}_{n},\vec{\theta}_{\mathrm{PL}},\vec{\theta}_{\mathrm{TPL}}\bigr{\}}\\ &=\bigl{\{}\vec{\theta}_{n},\vec{\theta}_{\mathrm{PL}},\{\vec{\theta}_{\mathcal{S},\mathrm{TPL}},\vec{\theta}_{f,\mathrm{TPL}}\}\bigr{\}}\\ &=\bigl{\{}\{N_{a},N_{p}\},\{\alpha_{\mathrm{PL}},\Omega_{\mathrm{ref,PL}}\},\{\vec{b}_{\ell m},\alpha_{\mathrm{TPL}},\Omega_{\mathrm{ref,TPL}},f_{\mathrm{break}},f_{\mathrm{scale}}\}\bigr{\}},\end{split}$$

$$\begin{array}{rl}{S}_{XX}(t,f|M_{\mathrm{1...}N};\vec{\theta })=& S_{XX}^{n,M_n}(t,f|{\vec{\theta }}_n)\\ & +{\mathcal{R}}_{XX}(t,f)S_{\mathrm{GW}}^{M_{\mathrm{PL}}}(f|{\vec{\theta }}_{\mathrm{PL}})\\ & +\int {\mathcal{R}}_{XX}(t,f,\mathbf{n})\mathcal{P}(\mathbf{n}|{\vec{\theta }}_{\mathcal{S},\mathrm{TPL}})d^2\mathbf{n}\times S_{\mathrm{GW}}^{M_{\mathrm{TPL}}}(f|{\vec{\theta }}_{f,\mathrm{TPL}}).\end{array}$$

$$\text{model=noise+powerlaw\_isgwb+truncatedpowerlaw\_sph}.$$

$$\text{model=noise+powerlaw\_isgwb-1+powerlaw\_isgwb-2}.$$

$${\Omega }_{\mathrm{GW}}(f)=\frac{2{\pi }^2}{3H_0^2}{f}^3{S}_{\mathrm{GW}}(f).$$

$$\Omega_{\mathrm{GW}}(f)=\Omega_{\mathrm{ref}}\left(\frac{f}{f_{\mathrm{ref}}}\right)^{\alpha_{1}}\bigg{(}1+\left(\frac{f}{f_{\mathrm{break}}}\right)^{\frac{1}{\delta}}\bigg{)}^{\delta(\alpha_{1}-\alpha_{2})},$$

$$S_h(f)=A{f}^{-7/3}{e}^{-f^{{\alpha }_{C17}}-\beta f\sin (\kappa f)}\left[1+\tanh (\gamma (f_k-f))\right]{\mathrm{Hz}}^{-1},$$

$${\Omega }_{\mathrm{GW}}(f)={\Omega }_{\mathrm{ref}}{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{2/3}{e}^{-f^{{\alpha }_{C17}}-\beta f\sin (\kappa f)}\left[1+\tanh (\gamma (f_k-f))\right]$$

$$\begin{split}S_{p}(f)&=N_{p}\left[1+\left(\frac{2\,\mathrm{mHz}}{f}\right)^{4}\right]\,\mathrm{Hz}^{-1},\\ S_{a}(f)&=\left[1+\bigg{(}\frac{0.4\,\mathrm{mHz}}{f}\bigg{)}^{2}\right]\,\left[1+\bigg{(}\frac{f}{8\,\mathrm{mHz}}\bigg{)}^{4}\right]\frac{N_{a}}{(2\pi f)^{4}}\,\mathrm{Hz}^{-1},\end{split}$$

Full text

I Introduction

The Laser Interferometer Space Antenna [LISA; 15, 37] is a joint ESA-NASA mission with a current launch date of 2035. It will be a spaceborne interferometric gravitational wave (GW) observatory in the mHz band, formed by a triangular constellation of three spacecraft 2.5 million km apart in a trailing (or leading) Earth orbit. As it occupies an entirely different frequency band to current-day GW observatories, LISA will probe a vast array of astrophysical phenomena as-of-yet unobserved in GWs. These include massive black hole binary (MBHB) mergers [98], extreme mass-ratio inspirals (EMRIs;[57, 101, 14]), stellar-origin compact binaries far from merger [49, 53, 97, 100, 21], tens of millions of Galactic and extragalactic double white dwarf binaries (DWDs) [79, 51, 95], and more; see Amaro-Seoane et al. [16] for a review.

In principle, any of the astrophysical GW sources discussed above could give rise to a stochastic gravitational wave background (SGWB) in LISA. In contrast to individually-resolvable, bright signals, SGWBs are formed from the superposition of many faint GW sources which collectively produce a GW confusion “noise”. SGWBs have been searched for in terrestrial gravitational-wave observatories [2, 3, 4, 5], and evidence for a nHz SGWB has been observed by pulsar timing arrays (PTAs) [10, 11, 12, 17, 18, 88, 89, 115, 111, 75]. Depending on the nature of the astrophysical or cosmological phenomena underpinning a given SGWB, the resulting signal can have different spectral characteristics, spatial morphology (i.e. anisotropy or lack thereof), and polarization content.111This work focuses on the spectral and spatial distributions of LISA SGWBs, and leaves consideration of polarization as a compelling direction of future work. As we will see, leveraging these facets of SGWB signals can prove useful for distinguishing the respective contributions of several such signals — a feat known as spectral separation. The ability to perform this kind of spectral separation will be crucial for SGWB science with LISA, as it is expected that the LISA datastream will contain several such signals.

Foremost among these is the stochastic GW contribution from every unresolved mHz DWD system in the Milky Way (MW), known as the Galactic foreground [79, 51, 95].222So named due to its prominence above even the LISA instrumental noise. For clarity and conciseness, we take the term “SGWB” to be inclusive of the Galactic foreground throughout this work, despite its “foreground” status. The Galactic foreground dominates the LISA band below a few mHz, where it rapidly truncates due to the subtraction of resolvable sources, source discreteness, and stellar evolutionary effects [79, 51, 95, 39, 65]. Due to its large amplitude, accurate characterization of the foreground signal is of paramount importance for detection and characterization of not just other SGWBs, but resolvable signals as well. In addition to its characteristic truncated spectrum, the Galactic foreground is highly anisotropic, tracing the spatial distribution of the Milky Way’s DWD population on the sky. This not only proves to be a valuable asset for spectral separation between the foreground and underlying SGWBs, but can also provide us with valuable insights into our Galaxy [e.g., 33, 113, 9].

The Milky Way is not, however, the only significant nearby population of DWD systems. Rieck et al. [91] show that the unresolved DWD systems in the Large Magellanic Cloud (LMC) produce a significant anisotropic SGWB (ASGWB) in LISA. Other MW satellites may contribute their own, lower-amplitude ASGWBs [87].

LISA is also expected to observe the astrophysical SGWB from extragalactic stellar-origin compact binaries. This background consists of binaries that will eventually merge in terrestrial GW detectors like that of the LIGO-Virgo-KAGRA (LVK) Collaborations [1, 7, 13]. However, these systems are far from merger while still in the mHz; as such, this SGWB is expected to possess an $\alpha=2/3$ power law spectrum [83].333$\alpha=2/3$ in $\Omega_{\rm GW}$; alternatively, $-7/3$ in PSD or $-2/3$ in characteristic strain. As the dominant contribution is expected to be from binary black hole mergers [e.g., 5], this signal is often referred to as the stellar-origin binary black hole (SOBBH) background; we will adopt this terminology moving forward. Current LVK rate measurements of the binary black hole (BBH) population [e.g., 6] can be used to predict the amplitude of the SOBBH SGWB. For a $\alpha=2/3$ power law with $f_{\mathrm{ref}}=25$ Hz, the SOBBH SGWB amplitude is expected to be $\Omega_{\mathrm{ref}}=7.3^{+4.5}_{-4.6}\times 10^{-10}$, as estimated for the LISA band by Babak et al. [21]; this estimate is consistent with those of other recent studies [e.g., 36, 48, 81, 5, 67]. The SOBBH SGWB will likely manifest as an isotropic SGWB (ISGWB) as seen in LISA.

Beyond these likely signals, there remain many other, more uncertain prospects. As discussed earlier, any resolved LISA source can have a corresponding background arising from its unresolved counterparts. However, whether the rate of the source in question is sufficient to generate a SGWB (let alone a detectable SGWB) in the LISA band is a different question. It is possible that LISA may observe a SGWB from EMRIs [30, 78, 86, 84, 112], MBHB mergers, or other populations of unresolved DWD systems beyond the MW and its satellites [60, 102]. Finally, a broad array of proposed cosmological SGWBs (i.e., inherently stochastic fields from the early universe) may be detectable (or at least able-to-be-constrained) by LISA; see Auclair et al. [20] for a review.

I.1 LISA Data Analysis

This plethora of sources carries simultaneously both great scientific potential and immense technical difficulty. In contrast to current terrestrial GW observatories which observe individual GW signals interspersed among long stretches of detector noise, the LISA datastream will instead be signal-dominated due to the sheer number of GW sources emitting constantly in the LISA band. In short, we will not know what LISA data looks like without GW signals in it, which lends considerable difficulty to the task of disambiguating GW signals from instrumental noise. This aspect of the LISA data analysis problem has motivated a “fit everything, everywhere, all at once” approach, dubbed the Global Fit [42, 47, 71, 70, 72].

The majority of current Global Fit prototypes [72, 64, 50] — with the exception of the maximum-likelihood approach of Strub et al. [103] — take the form of $\mathcal{O}(100,000)$-parameter transdimensional Bayesian analyses with a blocked structure. Their transdimensional nature allows for inference from the data of not just individual source parameters, but also the presence/absence and overall number of GW signals in the data. This is accomplished through Reverse Jump Markov Chain Monte Carlo (RJMCMC) sampling; such Global Fits are therefore necessarily Bayesian in nature. The blocked nature of the Global Fit structure comprises a high-level Gibbs sampler which allows for each source category (DWDs, MBHBs, EMRIs, SGWBs, etc.) to be treated by a purpose-made module, with information passing from module to module at every cycle of the sampler. Analyses which seek to constrain a single LISA source category must themselves be Bayesian as well, if they are to be integrated into a Global Fit setting of this kind.

Finally, due to otherwise-excessive amounts of laser phase noise, LISA data analysis will be performed using Time-Delay Interferometry (TDI) channels constructed virtually from the individual Doppler-tracking data of each LISA link [107]. The analyses considered in this work use the AET channels [107, 106].

I.2 LISA Multi-SGWB Analyses

Building on previous single-signal works [38, 23, 35, 108, 66, 104, 90, 40, 43, 8], a number of prototype analyses have been developed to analyze multiple SGWB signals in the LISA datastream. These include: SGWB searches in Global Fit residuals, following characterization of the foreground by the Global Fit itself [94]; two-component spectral separation of the Galactic foreground and a single background signal [29, 58, 69, 74, 91]; two-component spectral separation of an astrophysical and cosmological background [29]; three-component separation of a simplified Galactic foreground (i.e. sans treatment of anisotropy or its induced time-modulation), the SOBBH astrophysical background, and an underlying cosmological signal [85, 52, 32, 68]. These analyses, while promising in many respects, largely neglect to model anisotropies in either the Galactic foreground or potential underlying anisotropic signals.444The exceptions being Rieck et al. [91], which used BLIP to analyze the LMC and employed an early prototype of the framework presented here; Hindmarsh et al. [58], which considers a cosmological signal and a time-modulated Galactic foreground; and Pozzoli et al. [87], which leverages the cyclostationary statistics of the induced time-modulation of anisotropic SGWBs to separate the contributions from simple models of the Galactic foreground and the LMC. In doing so, they discard crucial information that can not only be used to characterize the MW itself [e.g., 33], but can possess significant utility when attempting to separate the anisotropic Galactic foreground from underlying isotropic and anisotropic SGWBs. There is a pressing need in the field for a SGWB analysis for LISA that is 1) Bayesian, so as to be compatible with the Global Fit; 2) capable of characterizing the anisotropic nature of the Galactic foreground and other anisotropic SGWBs like that of the Large Magellanic Cloud; 3) built for flexible spectral separation of multiple SGWBs in the LISA datastream; and 4) designed with eventual transdimensional RJMCMC SGWB analyses in mind.

In this work, we build on the Bayesian LISA Inference Package [BLIP; 23] to develop a flexible simultaneous inference formalism for both isotropic and anisotropic SGWBs in LISA. We leverage the anisotropic analysis capabilities of BLIP detailed in Banagiri et al. [23], Criswell et al. [46] and fully refactor the BLIP codebase to implement a fully modular structure — BLIP 2.0 — that not only enables flexible construction of spectral separation models for any number of isotropic and/or anisotropic SGWBs alongside simulation of the same, but has also been designed with an eventual transdimensional RJMCMC analysis in mind. The statistical formalism of this spectral separation infrastructure is developed in §II. Its implementation as BLIP 2.0 is described in §III. We provide proofs-of-concept for the BLIP 2.0 infrastructure via analyses of simulated LISA multi-SGWB data, including the Galactic foreground, the ASGWB from white dwarf binaries in the LMC, the isotropic SGWB from extragalactic stellar-origin binary black holes, and an underlying isotropic cosmological SGWB. The simulation procedure is discussed in §IV, and the results of the corresponding BLIP 2.0 spectral separation analyses are shown in §V. Discussion of the results and directions of future work are explored in §VI.

II Simultaneous Inference Formalism

II.1 BLIP Single-SGWB Analyses

We present in brief the BLIP likelihood for a single signal, as discussed in Banagiri et al. [23], Criswell et al. [46]; these derivations in turn built on those of previous works [40, 43, 8, 93]. For further details, refer to Banagiri et al. [23], Criswell et al. [46].

BLIP analyses divide the LISA time-domain datastream into segments of duration $T_{\mathrm{seg}}=10^{5}$ s, and perform a short-time Fourier transform with Hann windowing on each segment. Each segment is assumed independent and combined in a variant multidimensional complex Gaussian likelihood [8]:

[TABLE]

where $C_{IJ}(f,t|\vec{\theta})$ is the time- and frequency-dependent channel covariance matrix as a function of some model parameters $\vec{\theta}$. $f$ and $t$ here index over the $f^{th}$ frequency bin/$t^{th}$ time segment, respectively. $I$ and $J$ index over LISA’s three TDI channels (XYZ or AET). Unless otherwise noted, we use the AET channels. $\tilde{d}_{I}$ denotes the Fourier-domain data in channel $I$.

For a single SGWB in the presence of noise, we can write

[TABLE]

where $S^{n}_{IJ}$ and $S^{\mathrm{GW}}_{IJ}$ denote the noise and SGWB cross-/autopower spectral densities (PSDs), respectively. The former depends solely on the LISA instrument; we detail the simple LISA noise model used in this work in §III.2.5. The latter, however, accounts for both the intrinsic nature of the SGWB in question, as well as its interaction with the LISA detector response. In general, the SGWB power can be described in terms of its spectral and spatial distribution as $S_{\mathrm{GW}}(f,\mathbf{n})$, where $\mathbf{n}$ denotes direction on the sky, i.e., from LISA to the GW origin. We assume that the SGWB spectrum $S_{\mathrm{GW}}(f)$ and spatial distribution $\mathcal{P}(\mathbf{n})$ are separable; that is,

[TABLE]

where $\mathcal{P}(\mathbf{n})$ is normalized such that

[TABLE]

For an isotropic SGWB, $\mathcal{P}(\mathbf{n})=1/4\pi$. BLIP is capable of modelling anisotropic SGWBs in terms of a spherical harmonic expansion [23] or a pixelated skymap [46]. Depending on the nature of the an/isotropic model, the SGWB signal is convolved with the appropriate time/frequency/directionally-dependent LISA response $\mathcal{R}_{IJ}(f,t,\mathbf{n})$ such that

[TABLE]

For full derivations, we refer the reader to Banagiri et al. [23] for the spherical harmonic case, Criswell et al. [46] for the pixel-basis case, and Cornish [40], Cornish & Larson [43], Adams & Cornish [8], Romano & Cornish [93] for the general case.

II.2 Extension to Multiple Signals with Diverse Morphologies

We now consider the case of simultaneous inference of multiple SGWBs in the presence of instrumental noise, allowing for diverse SGWB spectral and spatial morphologies. The likelihood for this case is a straightforward extension of the single-signal likelihood of Eq. (1). As the covariance contributions of multiple (uncorrelated) signals in each channel are additive, we can write the likelihood for simultaneous inference of $N$ signal/noise (sub)models $M_{1...N}$ as

[TABLE]

with $\vec{\theta}=\{\vec{\theta}_{k}\}$ for $k\in[1...N]$ such that $\vec{\theta}_{k}$ are the parameters of signal/noise (sub)model $M_{k}$. “$\mathrm{SI}$” denotes the simultaneous inference case. The cumulative model covariance is defined as

[TABLE]

Note that only the spatial and spectral distributions of individual ASGWBs are assumed to be separable; the overall (i.e., cumulative) spatial and spectral distributions of SGWB power are no longer separable as defined by Eq. (3).

II.3 Example: Covariance for Three Components

For example, consider a simultaneous inference model consisting of the LISA instrumental noise and two SGWB signals: an isotropic power law — e.g., a cosmological or SOBBH ISGWB — and an anisotropic tanh-truncated power law with a spherical harmonic parameterized spatial model, motivated by simplified models of the Galactic foreground; see §III.2 for detailed discussion of these models. The overall parameter vector is then

[TABLE]

where PL and TPL refer to the power law and truncated power law (sub)models, respectively. A given channel’s model covariance contribution is then (using TDI X channel auto-power as an example):555Note that $\mathcal{R}_{XX}(t,f)$ in Eq. (9) is the integral of $\mathcal{R}_{XX}(t,f,\mathbf{n})$ for the isotropic case; see Eq. (5) and surrounding discussion.

[TABLE]

This result generalizes to all cross- and autopower channels to construct the full covariance matrix.

III BLIP 2.0 : a Flexible Simultaneous Inference Framework

As elegant in its simplicity as the formulation described in the previous section may be, its practical implementation in the general case is deceptively complex. While it is possible to implement individual cases of interest in the original BLIP infrastructure [as seen in 91], doing so repeatedly requires considerable duplicated effort and is ultimately far too inefficient for general use. As discussed in §I, there exist a wide variety of potential LISA SGWBs; these span varied spectral and spatial morphologies which can of course be described by any number of potential parametric (or even nonparametric) models. Moreover, flexibility in choice of both spectral and spatial models is highly valuable in a Global Fit context, where the RJMCMC nature of the Global Fit could in principle be used to perform SGWB model selection during the analysis. The ideal solution is then a framework in which the simultaneous inference formalism of §II is flexibly implemented for an arbitrary number and combination of signal models — as well as one in which implementing new signal models is simple and straightforward.

To this end, the BLIP codebase has been refactored from the ground up to enable full flexibility in both simulation and analysis of isotropic and anisotropic SGWBs. This is accomplished through a fundamentally modular design, dubbed BLIP 2.0 . All simulation and analysis routines have been rewritten to interact with a unified submodel class, along with the Injection and Model classes that marshal one or more submodels for use in simulation or analysis, respectively. Each instantiated submodel (with the exception of the LISA instrumental noise) consists of a spectral and spatial model, and carries with it all relevant information and functions related to one SGWB signal: its spectrum, its spatial distribution, the LISA response to it, its model parameters and prior distributions (in the case of an inference submodel), aesthetic details and names for plotting purposes, and so on. The submodel for the LISA instrumental noise is structured in a similar fashion, but does not contain a response function or spatial distribution. The Model class constructs the overall likelihood and prior distribution from a collection of submodels, and is the primary object (along with the data) that is passed to the sampler. For a single signal + instrumental noise, this likelihood is the one discussed in §II.1. In the case of multiple signals in the presence of instrumental noise, the Model class implements the simultaneous inference likelihood discussed in §II.2. The Injection class performs a similar function for simulations, combining the contributions from multiple submodels. Refer to Fig. 1 for a visual outline of the BLIP 2.0 infrastructure.

Within this new framework, to add a new spectral or spatial model to BLIP, all one must do is provide a name, a function governing its spectral or spatial distribution, the parameters for that function, its prior distribution, and — in the case of a spatial model — designate which of the LISA response routines it should connect to (i.e., isotropic responses for an ISGWB signal, anisotropic responses for an ASGWB). Anisotropic spatial models can be represented in either the pixel or spherical harmonic basis. Any combination of spectral and spatial models can be combined to form a submodel; this is accomplished easily by the end-user via a string within the BLIP input (“params”) file. For example, to simulate and analyze a power law ISGWB in the presence of LISA instrumental noise, one simply enters

[TABLE]

for both the Model and Injection items in the params file. These specifications need not be matched, and — as will be demonstrated in §IV-V — can be formed of any arbitrary combination of submodels. Relevant information for each submodel, such as the true parameter values for a simulated signal, can be passed by the user to BLIP via dictionaries in the params file where each submodel string operates as a key. This framework has enabled the implementation of a large number of new spectral and spatial models within BLIP; these models are detailed in §III.2.

Additionally, the Model/Injection/submodel structure not only allows relevant information to be seamlessly carried around BLIP itself for simulation, analysis, and plotting of results, but also enables data availability by saving the Model and Injection objects (and all submodels therein) alongside the posterior samples, plots, and so on. As a result, all information pertaining to a BLIP simulation/analysis run can be easily accessed after it is complete and uploaded to open data repositories such as Zenodo.

III.1 Dynamic Assembly of Simultaneous Inference Models in BLIP 2.0

The modular  Model / submodel infrastructure of BLIP 2.0 is uniquely suited to support simultaneous inference analyses. The BLIP Model object’s likelihood and covariance routines are designed to flexibly implement the simultaneous inference formalism described above with arbitrary combinations of submodels, across any number of signals.666With the following disclaimers: subject to available computing power; poorly- or over-specified models may not converge. Just because one can create a model with an arbitrary number of signals, doesn’t necessarily mean one should. As always, discretion is key. To the author’s knowledge at time of writing, BLIP is the only data analysis software for LISA with the capability to flexibly characterize arbitrary combinations of isotropic and anisotropic SGWBs.

To elucidate how this is accomplished in practice, consider again the example of §II.3. To perform an analysis using that simultaneous inference model — i.e., a power-law ISGWB paired with a truncated-power law ASGWB described by the spherical harmonic spatial model — we simply set

[TABLE]

BLIP will then construct a Model from these submodels by computing the appropriate LISA response functions for each,777With the exception of the LISA detector noise, which draws on its own set of unique routines and does not require the calculation of a response function. specifying the form of each submodel’s covariance contribution (and its dependence on that submodel’s spatial and spectral parameters), and combining those covariance contributions within the formalism of §II.2. Throughout this process, BLIP automatically handles all the requisite miscellanea to ensure that the sampler treats each submodel’s prior draws appropriately, parameter ordering is maintained throughout the code, duplicate model specifications are clearly differentiated, and so on. The model string given above can comprise of any number of ISGWB and ASGWB submodels. Duplicate submodels can be described by, e.g.,

[TABLE]

Every included submodel is ascribed a unique identifier, labeled as such on all plots produced at the end of an analysis run, and color-coded accordingly. This structure allows BLIP to seamlessly support almost any simultaneous inference analysis of interest, as will be seen in §V.

III.2 Currently Available Submodels

Table 1 summarizes the currently available spectral and spatial submodels within BLIP 2.0. The implementation within BLIP of many of these have been described elsewhere; see the listed references where relevant. The SGWB spectral models are parameterized in terms of the dimensionless GW energy density $\Omega_{\mathrm{GW}}$, but they are ultimately transformed into the (equivalent) strain PSD $S_{\rm GW}$. These quantities are related via

[TABLE]

Models not already described elsewhere are elaborated on below. Each model can be used flexibly for simulation and/or inference, and the spatial models can be implemented in either the pixel- or spherical harmonic basis. Only the spherical harmonic spatial model currently allows for inference of spatial parameters when used as an inference model; pixel-basis spatial models use the specified distribution as a fixed spatial template as described at length in Criswell et al. [46]. Enabling parameter inference for some of these latter spatial models is an ongoing effort beyond the scope of this work.

III.2.1 Broken Power Law (Spectral)

The broken power law spectral model,

[TABLE]

consists of a power law spectrum with amplitude $\Omega_{\mathrm{ref}}$ and slope $\alpha_{1}$, which turns over at the break frequency $f_{\mathrm{break}}$ to a second power law with slope $\alpha_{2}$. The sharpness of the turnover is mediated by a smoothing parameter $\delta$, which is usually fixed to a fiducial parameter rather than inferred. The amplitude $\Omega_{\mathrm{ref}}$ is normalized such that it is the projected amplitude of the first power law at $f_{\mathrm{ref}}$ Hz. This model is a good approximation for the LMC DWD ASGWB established in Rieck et al. [91].

III.2.2 Analytic Foreground (Spectral)

This model is a direct adaptation of the form given in Cornish & Robson [39], originally proposed by Babak as noted in [39]:

[TABLE]

which gives the Galactic foreground PSD as a function of an amplitude $A$ and a set of shape parameters $\{\alpha_{\mathrm{C17}},\beta,\kappa,\gamma,f_{k}\}$888$\alpha_{\mathrm{C17}}$ is given as $\alpha$ in Cornish & Robson [39]; the “C17” subscript included here is to disambiguate between this variable and the more general power law slope $\alpha$. which evolve over time as more individual DWD systems are resolved out of the Galactic foreground; see Cornish & Robson [39] for further details. As implemented in BLIP, this spectral model is recast in terms of $\Omega_{\mathrm{GW}}$:

[TABLE]

and adjusts the shape parameters based on the user’s choice of observing duration.

III.2.3 Analytic MW Satellite (Spatial)

The Analytic MW Satellite submodel was initially implemented as a simplified model of the LMC, but its parameters can be modified to serve as an analog of any approximately spherical dwarf galaxy satellite of the MW.999i.e., any one of the nearby dwarf galaxies which orbit the MW: the Large and Small Magellanic Clouds, Sagittarius Dwarf, etc.. It is an appropriate model for any DWD ASGWB that would be better represented by a spatially extended distribution rather than a point source. It functions similarly to the Analytic Galaxy model in that it creates a parameterized 3D density distribution which is then projected onto the sky as seen by LISA. It creates a filled sphere of evenly spaced points centered at the location specified by its RA and DEC parameters, given in degrees and corresponding to standard ICRS right ascension and declination; at the specified distance $d$ from the Solar System Barycenter given in kpc. The radius of the sphere is set by the parameter $r$, given in kpc. GW power is evenly distributed among these points and projected onto a Healpix skymap.

III.2.4 Point Source/Radiometer (Spatial)

The point source spatial submodel consists of a variable number of individual bright pixels on a Healpix [54, 116] skymap; these can be contiguous or non-contiguous as desired. All pixels within a single point source submodel share a spectrum. The user can create multiple spectrally-independent point sources by instantiating several point source submodels, each containing a single source with its own individualized spectral submodel. When used for inference, this submodel effectively becomes a radiometer analysis [22, 76] and searches for stochastic GW power emanating from a particular sky direction or set of directions.

III.2.5 Noise Spectral Model

The current LISA instrumental noise model supported by BLIP101010This noise model has been detailed in the context of BLIP elsewhere [23, 91, 46, 27], but is included here for completeness. is the two-parameter model given in Amaro-Seoane et al. [15]:

[TABLE]

where $S_{p}(f)$ and $S_{a}(f)$ describe the contributions due to position and acceleration noise, parameterized respectively by the noise amplitudes $N_{p}$ and $N_{a}$. This noise model is used for both noise simulation and inference. The actual LISA instrumental noise is expected to be significantly more complex than what is considered here. Detailed, realistic models of the LISA instrumental noise are under development [72, 24, 56, 80] and should be incorporated in future; however, such efforts are beyond the scope of this work.

III.2.6 Spectral Models with Fixed Parameters

For the extant models discussed above, it is often beneficial to fix (i.e., hold constant and not sample over) one or more parameters. This can be desirable for several reasons. From an astrophysical perspective, we expect (for example) the SOBBH SGWB to have a slope of $\alpha=2/3$; it may behoove us to search for a power law ISGWB in LISA with that slope. Similarly, DWD SGWB spectra are expected to exhibit a $2/3$ slope at low frequencies before the effects of discreteness and binary interactions become relevant. Provided the Universe conforms to our expectations, fixing one or more parameters to their physically-motivated values can make our search more efficacious. From a pure modelling perspective, one may wish to prescribe that certain parameters like the truncated power law’s $f_{\mathrm{scale}}$ — which have little astrophysical meaning and tend to not vary much within a given use case — be fixed to some fiducial value. BLIP ’s submodel structure flexibly allows for the creation of a submodel variant to any of these spectral models which fixes any number of model parameters. Some pertinent cases, like those with a fixed $\alpha=2/3$ slope, are hardcoded and available as out-of-the-box spectral submodels.

III.3 Optimization

Many BLIP analyses previously required a prohibitive amount of computational resources and time to perform. For instance, some anisotropic analyses — which must convolve the LISA responses with a model skymap while sampling — could take several months to converge within the original BLIP framework (assuming they did in fact converge). Even before the MCMC sampling step, computing the LISA responses for every time segment, across the full LISA frequency band, for every sky direction, and (in the spherical harmonic case) every spherical harmonic mode up to some truncation $\ell^{a}_{\mathrm{max}}$ to simulate a SGWB — let alone multiple — for the nominal LISA mission duration of 4 years was beyond the computational resources and time allocation available to the authors.

Many avenues have been pursued to push the limits of what BLIP can accomplish with the computational power available. At a base level, more efficient trade-offs between memory usage and computation time have been implemented in the response function calculations, and the new modular, streamlined structure of BLIP 2.0 alleviates many instances of repeated calculations. Beyond these basic coding improvements, more involved acceleration techniques have been implemented within BLIP, each falling into one of three categories: parallelization, just-in-time compilation with JAX, and GPU-acceleration.

III.3.1 Parallelization

Parallelization has been implemented across several modules of the BLIP codebase. The LISA response function calculations — far and away the most computationally-taxing component of BLIP — are now fully parallelized along their frequency axis. This has reduced the required computation time for the time-dependent LISA response by over an order of magnitude for the most complex analyses considered in this work.111111With an even greater reduction for less intensive analyses; total available RAM remains a fundamental limitation on the number of active parallel processes. Additionally, each submodel can be computed in parallel to the others,121212(although this is mutually exclusive with the usually-superior response function parallelization — the latter being an improvement on the former — due to Python’s multiprocessing limitations) and parallel sampling routines have been implemented to work with the dynesty and Numpyro samplers.

III.3.2 Just-in-Time Compilation with JAX

All parts of the BLIP code that interact with its sampling routines (i.e., the Model and submodel objects and everything attached to them) are now fully compatible with JAX. JAX [31] is a Python library that — modulo several specific design requirements for user code131313See the JAX documentation. — serves as an approximately 1:1 swap in for Numpy [55] with several powerful additions. First among these is just-in-time compilation (JIT) of Python code. Python is not a compiled coding language; as such, it gains significant flexibility in exchange for reduced efficiency. JAX’s JIT compiles code just as it is about to be executed; the resulting speed-up is usually negligible for one-time operations, but for functions that will be called repeatedly,141414A pertinent example being the likelihood and prior given to a sampler. the improvement can be dramatic. Second, JAX enables automatic differentiation, which allows the gradient to be taken of any JAX-able Python function; crucially, this includes the BLIP likelihood, allowing its gradient with respect to its (potentially many) parameters to be evaluated and enabling the use of Hamiltonian Monte Carlo (HMC) sampling with BLIP . Due to BLIP’s new modular structure, any subsequent spatial or spectral models added to BLIP in future will automatically be JAX-compliant. Finally, JAX is also designed to function seamlessly across both CPU and GPU settings.

III.3.3 GPU-Accelerated Sampling with Numpyro

These latter points lead to the final major speed improvement for BLIP analyses: the implementation of GPU-accelerated HMC sampling with Numpyro [82]. Numpyro is built with JAX in mind, using the latter’s automatic differentiation to resolve the central challenge of HMC (i.e., needing to know the gradient of the likelihood) and leveraging JAX’s innate GPU compatibility to do so at breathtaking speeds. On Nvidia A100 GPUs, convergence times for ISGWB and low-resolution spherical harmonic/pixel-basis ASGWB analyses for one year of LISA data were reduced from order days to weeks to order minutes to hours. As will be seen throughout the rest of this work, most analyses over the nominal LISA mission data duration of 4 years have been brought in reach, taking $\sim$a day on an Nvidia A100 and $\mathcal{O}$(tens of minutes — a few hours) on an Nvidia GH200 superchip. The sole exception to this new computational horizon is high-resolution, parameterized anisotropic searches such as that of the spherical harmonic analysis at high $\ell^{a}_{\mathrm{max}}$. The final response functions for such analyses carry RAM requirements ranging from hundreds of GBs to a few TBs; a far cry from being able to fit within the available 80GB of RAM on a A100 GPU, or even the 141GB of an Nvidia H200 GPU. The GH200 superchip’s shared memory presents promising solution to this limitation and can in principle enable GPU-acceleration of parameterized, high-resolution anisotropic analyses. At present, however, JAX does not yet support JIT with arrays distributed across shared memory; making use of this powerful resource will require further development that we leave for future work. Finally, computation of the time-dependent directional LISA response function and time-domain data has also been GPU-accelerated, with simulation of a multi-component, 4-yr dataset reduced to $\mathcal{O}$(hours); shorter and/or less-complex datasets can take as little as a few minutes for the entire data simulation process.

The advances highlighted here represent a factor of $\mathcal{O}(10^{4})-\mathcal{O}(10^{5})$ improvement in efficiency, and have achieved a major milestone in BLIP’s development: the runtime of most BLIP analyses is now computationally feasible in a Global Fit setting, fitting many times over within, e.g., the $\sim$week runtime of the similarly GPU-accelerated Erebor [64]. The stage is set for inference of multiple isotropic and anisotropic SGWBs within the Global Fit itself.

III.4 Utility Features

In addition to the full code restructure, a number of utility and quality-of-life features have been added to BLIP. These include:

Checkpointing of all analyses at a user-defined interval to preserve progress in the event a long run is interrupted. 2. 2.

The ability to resume sampling of chains from the final state of a completed run. This is useful for MCMC samplers like Numpyro that do not have the well-defined convergence criteria of nested sampling. 3. 3.

The ability to produce and save a simulation for later analysis, as opposed to the previous requirement of a joint simulation + analysis run; this saves a not-insignificant amount of computational expense, especially for repeated analyses of the same dataset. Doing so stores not only the produced time-domain dataset, but also the entire Injection object, allowing BLIP to carry forward all pertinent simulation information for evaluation of analysis results. 4. 4.

All computed simulated and estimated spectra and skymaps used to create the final set of BLIP plots are saved to a separate compressed file, the contents of which can be easily accessed outside of the BLIP structure and its dependencies. This feature greatly eases the process of creating and reformatting plots of BLIP results. 5. 5.

Many other small aesthetic and ease-of-use improvements, such as: dynamic sizing of corner and spectral plots; the ability to specify a skymap color scheme for anisotropic recoveries; the ability to alias simulation submodels to those used in the analysis for purpose of plotting, carrying around “true” (simulated) parameter values, etc.; and more.

IV Simulations

The modular framework described in §III also enables simulation of LISA datasets containing — in principle — any number and combination of isotropic and anisotropic SGWBs. Within BLIP, these are referred to as simulation (or injection) “components”. Realistically, the number of components is of course limited somewhat by scaling computational expense, as well as our expectations of what will actually be present in LISA data. From a computational perspective, the primary limitation is the need to compute the full time-dependent directional LISA response functions $\mathcal{R}_{IJ}(f,t,\mathbf{n})$ for every simulated SGWB. At appreciable timescales and high fidelity of the simulated Healpix skymaps, calculating the LISA response carries with it the vast majority of the computational expense of a given simulation. To mitigate this, BLIP dynamically shares information between individual simulation components to avoid repeated calculations and ensure that only the necessary calculations are performed. For example, if one wishes to simulate two point source SGWBs, each comprising a single pixel on the sky, BLIP will only compute the underlying satellite orbits, transfer functions, etc. once, and will innately limit the directional calculations to the two pixels with power (as the remaining pixels have zero power and will sum to zero when convolving the directional LISA response functions with the SGWB skymaps; see Eqn. (5)). At present, BLIP is able to simulate multi-SGWB datasets spanning a reasonable expectation of what will be present within LISA data151515With the exception of a potential non-stationary, non-Gaussian SGWB from extreme mass-ratio inspirals [e.g., 34, 112, 78], as BLIP does assume both Gaussianity and stationarity; see further discussion in §VI. for a duration of 4 years and with a Healpix nside (pixel resolution) of 16. This capability is demonstrated as follows.

IV.1 Simulation Components

We use the BLIP 2.0 framework to simulate several multi-SGWB datasets for use as test cases. These consist of the following components; for each, see §III.2, Table 1, and references therein for detailed descriptions of the underlying spectral and spatial models. Model parameter values for the simulated SGWBs discussed here are given in Table 3.

IV.1.1 Simple Galactic Foreground (“MW”)

We simulate a simple Galactic foreground, combining the smooth Tanh-truncated Power Law spectrum with a parametric model for the Galactic spatial density distribution (Analytic Galaxy) consisting of a simple bulge and exponential disk. The simulated Galaxy has radial and vertical scale heights of $r_{h}=2.9$ kpc and $z_{h}=0.3$ kpc, respectively (the “thin disk” model of [33], based on previous models of Binney et al. [25], Bissantz & Gerhard [26], McMillan [73]). BLIP is capable of simulating and modeling more-realistic Galactic foregrounds produced from binary population synthesis catalogues (e.g., popmap in Table 1 as discussed in [46, 91]). However, performing detailed spectral separation between a realistic Galactic foreground and underlying SGWBs will comprise a significant effort beyond the scope of this work. A simple Galactic foreground better suits the current purpose of demonstrating BLIP’s simulation and analysis capabilities — as the simulated and inferred parameter values can be directly compared — but we will explore the challenge of spectral separation with a realistic foreground in future work (see §VI).

IV.1.2 Isotropic SOBBH SGWB (“SOBBH”)

We simulate the SGWB from extragalactic stellar-origin binaries as an isotropic signal with a power-law slope of $\alpha=2/3$ and the upper range amplitude given in Babak et al. [21]. This latter work provides projections of the LISA SOBBH SGWB based on the population of LVK-observed compact binary mergers.

IV.1.3 Simple LMC (“LMC”)

We simulate a simple version of the anisotropic SGWB from white dwarf binaries in the LMC [91]. This is an analytic simulation of the LMC SGWB for the same reasons discussed above with respect to the Galactic foreground. We simulate the LMC SGWB spectrum as a broken power law, due to its notable turnover but lack of rapid truncation; see, e.g., Figure 2 of Rieck et al. [91]. The LMC spatial distribution is simulated as a uniform-density sphere with the observed radius of the LMC.

IV.1.4 Cosmological SGWB (“CGWB”)

We additionally simulate a generic cosmological SGWB. It is isotropic, lower in amplitude than the SOBBH SGWB, and possesses a flat power spectral density — representative of a generic scale-invariant inflationary background. Its primary purpose is to demonstrate BLIP’s ability to place upper limits on such a signal in the presence of a variety of higher-amplitude isotropic and anisotropic SGWBs.

IV.2 Simulated Datasets

We consider two example datasets, each consisting of combinations of the SGWB components described above. Each dataset spans 4 years, corresponding to $1.262\times 10^{8}$ s of simulated data. Anisotropy is simulated using Healpix skymaps at an nside (pixel resolution) of 16. The reference frequency used throughout is $f_{\rm ref}=1$ mHz. We assume Gaussian noise of the spectral form given in §III.2.5 with parameters $N_{a}$ and $N_{p}$ as given in Amaro-Seoane et al. [15] (see Table 3). We simulate time-domain LISA data by splicing together inverse Fourier transforms of the frequency-domain Gaussian contribution from each simulation component, convolved with the time-dependent directional LISA response across each considered TDI channel, on overlapping, Hann-windowed splice segments of duration $10^{-4}$ s. Further details can be found in Banagiri et al. [23].

The first dataset, “MW + LMC + SOBBH” combines a set of likely SGWB signals: the anisotropic Galactic foreground, the highly localized ASGWB from white dwarf binaries in the LMC, and the ISGWB from extragalactic stellar-origin binaries (SOBBH). The second dataset “MW + LMC + SOBBH + CGWB” additionally includes an low-amplitude isotropic cosmological background to investigate our ability to place upper limits on an underlying SGWB signal of cosmological origin in the presence of multiple foreground signals. The data simulation process is fully GPU-accelerated, and takes $\sim 5$ hours to calculate the full $3\times 3\times$ frequency $\times$ time $\times$ sky direction LISA response tensor and simulate time-domain data for a multi-component, 4-year dataset at the chosen skymap resolution on a single Nvidia GH200.

V Results

V.1 MW + LMC + SOBBH

We analyze the 3-component dataset as follows. We use the Tanh-truncated Power Law spectral model for the MW spectrum, and the Analytic Galaxy pixel-basis template for the MW spatial model [see Table 1 and 46]. The MW model Galaxy morphology is matched to that of the simulation, i.e. with scale height parameters set to $r_{h}=2.9$ kpc and $z_{h}=0.3$ kpc; we can make this assumption as the spatial distribution of the DWD population will be well-measured by the resolved DWDs [e.g., 110]. The LMC spectrum is modelled with a broken power law; the LMC spatial distribution is modelled with the Analytic Satellite pixel-basis template (see §III.2.3). We use a power law ISGWB to model the SOBBH ISGWB, with a fixed slope of $\alpha=2/3$ to target the search for the expected shape of the SOBBH spectrum. We use the noise model described in §III.2.5. The priors for each search are astrophysically-motivated so as to target their respective signal; these priors are given in Table 4. The inferred posterior distributions are shown in Fig. 2, and the recovered spectra are shown in Fig. 3.

The Galactic foreground spectrum is precisely and accurately recovered; this is unsurprising given its prominence, the matching functional form of the simulated spectrum, and the strong prior enforced by the fixed-template spatial analysis. However, the agreement of the inferred parameters with the simulated values as seen in Fig. 2 provides an important consistency test. That being said, the recoveries of the lower-amplitude signals are more interesting. First, the LMC is overall recovered well. The templated anisotropic search allows for recovery of this more-complex spectral model compared to the power law + spherical harmonic analysis of Rieck et al. [91], and recovery of the LMC spectrum’s low-frequency amplitude and slope is much improved compared to that initial study. Its (log) amplitude is measured to within 10%, indicating an intriguing potential for comparison of the LMC and Galactic SGWB spectra — and therefore the masses, number, and typical distances of their underlying DWD populations. Adding to this potential is the precise and accurate inference of the turn-over frequencies for both the MW and LMC spectra. For realistic signals, this phenomenological quantity is driven by both characteristics of these astrophysical populations and detector-dependent effects (vis-à-vis resolving out individual DWDs). If these effects can be separated, the different turnovers of the LMC and MW spectra may provide additional angles through which to investigate the underlying astrophysics of their separate DWD populations.

Conversely, we see significant spectral mixing and a corresponding bias in the recovery of the second LMC power law below the turnover frequency $f_{\rm break}$ (i.e., for $f\gtrsim 5$ mHz). As can be seen in Fig. 3, the spectrum is accurately inferred around the knee, until the noise level raises substantially and the amplitudes of the LMC and SOBBH signals become comparable. Further investigation is needed to determine if this bias presents a serious issue at higher frequencies, especially when considering realistic spectra from actual (simulated) DWD populations which are generally much more complex at higher frequencies. We note this as a direction for future work.

Finally, the simulated SOBBH SGWB log amplitude is inferred to be $\log_{10}\Omega_{\rm GW}=-11.8^{+0.1}_{-0.2}$ (95% C.I.), in agreement with the simulated value of $-11.68$. Even when considering multiple astrophysical foregrounds, this constraint improves upon the expectations derived from observations of compact binary mergers as of GWTC-3; c.f. the range calculated in Babak et al. [21] of $-12.0^{+0.3}_{-0.5}$ (converted here to units of $\log_{10}\Omega_{\rm GW}(1\ \mathrm{mHz})$ for the sake of comparison). Continued LVK observations will of course continue to narrow these predictions; agreement (or lack thereof) between LISA and LVK observations of these systems will provide an exciting window into the evolution of stellar-origin compact binaries. Conversely, projections from LVK-observed mergers [e.g., 5] are necessarily model-dependent; combining information across both LISA and LVK data is a promising route for future multi-band science. It is worth noting that both these and previous constraints on the SOBBH, LMC, and MW spectra depend on accurate characterization of the LISA instrumental noise. Here, we assume knowledge of the functional form of the noise spectrum; this will not be the case when LISA flies, and as such these estimates can be expected to degrade somewhat under realistic conditions [56, 77].

V.2 MW + LMC + SOBBH + CGWB

We now analyze the four-signal MW + LMC + SOBBH + CGWB dataset. We make identical modelling choices as in the previous analysis, with the inclusion of an additional isotropic power law model. We perform a targeted search for a power law of fixed spectral index $\alpha=1$ with a broad prior on the amplitude; all priors used are given in Table 4. We place an upper limit on the amplitude of the simulated CGWB in the presence of astrophysical foregrounds of $\log_{10}\Omega_{\rm GW}(1{\rm\ mHz})\leq-12.0$ ($97.5\%$ credible upper bound). This is roughly consistent with the result of e.g., Boileau et al. [28], who are able to recover a CGWB of amplitude $\log_{10}\Omega_{\rm GW}(1{\rm\ mHz})=-12.097$ (although the comparisons are not precisely 1:1 due to different data durations, slopes, amplitudes of astrophysical foregrounds, inclusion of the LMC, and so on). However, detailed consideration of CGWB detection thresholds in the presence of multiple astrophysical foregrounds warrants a dedicated investigation beyond the scope of this work.

It is worth noting that due to the isotropic nature of both signals, the inclusion of the CGWB increases the uncertainty of the inferred SOBBH amplitude. While the uncertainty of the MW and LMC estimates are unchanged by inclusion of the underlying CGWB, the SOBBH estimate degrades slightly: from $\Omega_{\rm GW}(1{\rm\ mHz})=-11.8^{+0.1}_{-0.2}$ (95% C.I.) in the three-signal case to $\Omega_{\rm GW}(1{\rm\ mHz})=-11.8^{+0.1}_{-0.6}$ (95% C.I.) in the four-signal case. Power is able to be freely traded between the two isotropic signal models, and the SOBBH inference suffers as a result. However, it is important to note that while the inclusion of the CGWB increases the uncertainty of the SOBBH amplitude estimate, it does not induce a bias.

V.3 Effect of Neglecting Anisotropy in Spectral Separation

We additionally, consider the impact of neglecting the anisotropy of the MW and LMC signals. We repeat the MW+LMC+SOBBH analysis of §V.1 keeping all aspects of the analysis identical save that the MW, LMC, and SOBBH SGWBs are all modelled with isotropic spatial distributions. The resulting inferred spectra are shown in Fig. 5. Failing to model the intrinsic anisotropy greatly increases spectral mixing and induces significant systematic biases in the inferred spectra. The MW, LMC, and SOBBH amplitudes are overestimated by 2%, 830%, and 376%, respectively, with each inferred amplitude completely excluding the corresponding simulated value. These biases arise from a combination of unmodeled low-frequency power due to the anisotropy-induced time-modulation and absorption of instrumental noise into the SGWB power, as the acceleration noise is underestimated by 86% whereas the position noise estimate is unbiased (the latter being dominant at higher frequencies). The severity of these effects dramatically underlines the crucial importance of anisotropy for accurate SGWB spectral separation with LISA.

VI Discussion and Conclusions

We present BLIP 2.0, a flexible framework for Bayesian inference of multiple isotropic and anisotropic SGWBs in LISA. Its modular structure allows for assembly of inference models composed of arbitrary numbers and morphologies of SGWB signals, up to computational limitations. The same structure likewise enables simulation of multi-component LISA SGWB datasets. Accounting for signal morphology enhances spectral separation of multiple signals, and allows for spectro-spatial separation of multiple ASGWBs. This is of particular importance for LISA due to the prominent anisotropic Galactic foreground and other local-Universe ASGWBs such as the LMC. We present two multi-component simulated datasets and Bayesian SGWB analyses thereof to demonstrate the capabilities of the framework: the first consisting of the Galactic foreground, the LMC SGWB, and the SOBBH ISGWB, and the second additionally including an underlying cosmological SGWB. We briefly demonstrate the adverse effects of neglecting anisotropy when attempting spectral separation of LISA SGWBs. The BLIP 2.0 framework is the first of its kind, and greatly expands our analysis capability for ASGWB analyses and spectro-spatial separation of multiple isotropic and anisotropic SGWBs with LISA. Additionally, its Bayesian approach and modular infrastructure are well-suited for integration into Bayesian Global Fits.

There remain some limitations of the code. It currently operates on a series of short-time Fourier transforms, which are less robust to non-stationary noise and data gaps; implementation of the time-frequency basis [41] is a logical next step. Doing so would also enable exploration with BLIP of non-stationary SGWBs like that expected from an EMRI SGWB, and aid in Global Fit integration. Additionally, the noise modelling used in the demonstrations presented here is quite simple compared to our expectations for the LISA instrumental noise; implementation of more advanced noise modelling is warranted, and will be supported by BLIP’s modular framework.

Looking forward, it will be important to consider the case of non-idealized spectra. The Galactic foreground will not have a smooth spectrum, and will include features from both astrophysical populations [e.g., 96] and potential non-Gaussianities from individual subthreshold DWDs [94]. This will complicate spectral separation, especially for other DWD-driven SGWBs like that of the LMC. In this work we were able to perform spectral separation between simple models of the MW and LMC; additional work will be needed to approach the realistic case of spectro-spatial separation between two DWD populations’ SGWBs. This is an exciting direction of future work that will be greatly aided by BLIP’s existing spatial models and straightforward support for new spectral spatial models. The results presented here indicate strong potential for LISA’s ability to characterize important morphological aspects of the MW and LMC spectra, opening a path towards astrophysical comparison between these two local DWD populations.

The challenge of properly modelling the stochastic spectra of the MW and LMC arises from their very nature as astrophysical populations. Simple phenomenological models like those employed in this work, despite accurately describing the overall spectral shape of the MW and LMC, fail to incorporate the detailed effects of bin-to-bin Poisson fluctuations in the number and amplitude of the contingent unresolved DWDs. Moreover, the phenomenological truncation mechanisms are intrinsically unable to account for the departures from true stochasticity as source discreteness becomes relevant [e.g., 99]. On top of this, such models are difficult to interpret from an astrophysical perspective, as their shape can be just as influenced by detector-dependent resolution of individual binaries as they are by the underlying astrophysical population. Taken together, these factors strongly motivate development spectral models for the Galactic foreground (as well as the LMC) which directly arise from treatment of the unresolved astrophysical population. Not only would such models allow for more flexible and accurate characterization of realistic DWD-driven SGWBs, but also for astrophysical population inference from the Galactic foreground and — potentially — the LMC SGWB. In the case of the Galactic foreground, joint population inference with the resolved binaries may even allow for stronger constraints on the Galactic foreground (c.f. a resolved + unresolved population model for next-generation terrestrial detectors [114]), thereby improving almost every LISA science case.

Finally, BLIP’s modular infrastructure provides a strong foundation for transdimensional multi-SGWB modelling, wherein new submodels are proposed and removed, and the data itself can inform the choice of both spectral and spatial modelling. Work is underway towards this end, beginning with the implementation in BLIP of Eryn [62, 63], the transdimensional GPU-accelerated sampler used in the Erebor global fit [64].

Data Availability

BLIP is open source; the primary repository can be found at https://github.com/sharanbngr/blip. The specific version of the code used in this study can be found within the corresponding BLIP release on Zenodo [45]. Simulated data and posterior samples for all analyses are available on Zenodo [44].

Acknowledgments

Portions of this manuscript were adapted from the doctoral dissertation of AWC, Astrophysical Inferences from Multimessenger Ensembles. AWC acknowledges support by NSF grant no. 2125764 and NASA grant 90NSSC19K0318. S.B. acknowledges support from NSF grant PHY-2207945 and from Australian Research Council (ARC) grant CE230100016. The authors are members of the NANOGrav collaboration, which receives support from NSF Physics Frontiers Center award number 1430284 and 2020265. SRT acknowledges support from an NSF CAREER no. 2146016, NSF AST-2007993, and NSF AST-2307719. SRT also acknowledges support from a Chancellor’s Faculty Fellowship from Vanderbilt University. The authors acknowledge the use of computing resources provided by the Minnesota Supercomputing Institute at the University of Minnesota. This work leveraged the resources provided by the Vanderbilt Advanced Computing Center for Research and Education (ACCRE), a collaboratory operated by and for Vanderbilt faculty.

Software

BLIP [this work; 23, 46], numpy [55], scipy [109], ChainConsumer [59], astropy [92, 105, 19], healpy [54], JAX [31], Numpyro [82], and Matplotlib [61].