A pipeline for searching and fitting instrumental glitches in LISA data

Martina Muratore contact: [email protected] \addressii Jonathan Gair \addressii Olaf Hartwig \addressaeihan \addressunihan Michael L. Katz \addres Alexandre Toubiana \addressio \addrese

(May 26, 2025)

Abstract

Instrumental artefacts, such as glitches, can significantly compromise the scientific output of LISA. Our methodology employs advanced Bayesian techniques, including Reversible Jump Markov Chain Monte Carlo and parallel tempering to find and characterize glitches and astrophysical signals. The robustness of the pipeline is demonstrated through its ability to simultaneously handle diverse glitch morphologies and it is validated with a “Spritz”-type data set from the LISA Data Challenge. Our approach enables accurate inference on Massive Black Hole Binaries, while simultaneously characterizing both instrumental artefacts and noise. These results present a significant development in strategies for differentiating between instrumental noise and astrophysical signals, which will ultimately improve the accuracy and reliability of source population analyses with LISA.

I Introduction

The space-based gravitational wave detector, the Laser Interferometer Space Antenna (LISA) mission Amaro-Seoane et al. (2017); Colpi et al. (2024) was adopted as a mission by the European Space Agency in January 2024 and is scheduled for launch in 2035. One of its main objectives is the first direct observation of coalescing massive black holes binaries (MBHBs) with gravitational waves (GWs). Such observations will transform our understanding of the astrophysical processes leading to the formation of massive black holes and driving their coevolution with their host galaxies Gair et al. (2011); Sesana et al. (2011); Klein et al. (2016); Dayal et al. (2019); Barausse et al. (2020a); Toubiana et al. (2021); Chen et al. (2022); Fang and Yang (2023); Spadaro et al. (2025); Langen et al. (2024); Toubiana et al. (2024a). Moreover, these observations will reach signal-to-noise ratios (SNRs) of thousands, allowing for exquisite tests of general relativity Berti et al. (2006, 2016); Chamberlain and Yunes (2017); Barausse et al. (2020b); Maggio et al. (2020); Bhagwat et al. (2022); Corman et al. (2022); Toubiana et al. (2024b); Pitte et al. (2024). Finally, the possible observation of electromagnetic counterparts Mangiagli et al. (2022); Dong-Páez et al. (2023); Izquierdo-Villalba et al. (2024) would allow MBHBs to be used as “bright sirens” for the measurement of cosmological parameters, not only the Hubble constant, but possibly also the equation of state of dark energy, thanks to the wide redshift range of these observations Caprini and Tamanini (2016); Tamanini et al. (2016); Mangiagli et al. (2023). To achieve these science objectives, it is crucial to develop a robust data analysis pipeline for LISA.

Unlike current ground-based GW detectors such as LIGO LIGO Scientific Collaboration et al. (2015), VIRGO Acernese et al. (2015) and KAGRA Somiya (2012); Aso et al. (2013). LISA will be a signal-dominated detector with multiple GW signals overlapping in time and frequency. Indeed, in addition to MBHBs, signals from $\sim 10^{7}$ Galactic binaries are expected to be present in the data stream, of which $\sim 10^{4}$ will be individually resolved, and the rest will form a stochastic noise that will dominate over the instrumental noise around 1 mHz Nelemans et al. (2001); Ruiter et al. (2010); Korol et al. (2017); Lamberts et al. (2019); Korol et al. (2022); Toubiana et al. (2024c). The extreme-mass-ratio inspirals (EMRIs) of stellar-origin compact objects into massive black holes in the centres of galaxies will also generate signals in the LISA data that could last for the whole mission duration. The rate of EMRI events is very uncertain, ranging from a few to tens of thousands Gair et al. (2010); Amaro-Seoane (2018); Babak et al. (2017); Pan and Yang (2021); Seoane et al. (2024), and at the upper end of this range EMRIs could perhaps even form a stochastic foreground Bonetti and Sesana (2020); Pozzoli et al. (2023). Finally, the loudest and most massive of the stellar-origin black hole binary population being observed by ground-based detectors are also expected to be observable during their inspiral phase Sesana (2016); Wong et al. (2018); Gerosa et al. (2019); Toubiana et al. (2020); Buscicchio et al. (2021); Toubiana et al. (2022); Buscicchio et al. (2025) and a variety of cosmological backgrounds could significantly contribute to the total noise (see Auclair et al. (2023) and references therein). Moreover, we will not have a direct and independent measurement of the instrument noise Muratore et al. (2023), meaning that we will have to estimate its properties simultaneously with the parameters of the astrophysical sources. The simultaneous characterization of all the sources of multiple different types that are expected to be present in the LISA data, along with instrumental properties, is usually referred to as the LISA Global Fit Cornish and Crowder (2005); Littenberg and Cornish (2023); Katz et al. (2024); Strub et al. (2024); Deng et al. (2025).

When designing a methodology to characterize instrumental noise, it is crucial to avoid systematic errors arising from noise misinterpretation or mismodeling. The LISA precursor mission, LISA Pathfinder (LPF) Armano et al. (2024a), which shares technology with LISA, was affected by spurious signals of unknown origin, referred to as glitches. Since various types of glitches were detected in LPF Armano et al. (2022a), it is reasonable to expect that, unless they are fully understood and mitigated through hardware requirements, similar glitches may appear in LISA. The LISA hardware, particularly the Gravitational Reference Sensors (GRS), will be largely the same as in LPF, with only minor modifications to accommodate integration into the Movable Optical Sub-assembly (MOSA) and the LISA spacecraft Muratore (2021). Therefore, we anticipate encountering glitches similar to those seen in LPF—especially force impulses—which could impact LISA’s sensitivity.

Glitches are also commonly seen in the data collected by ground-based gravitational wave detectors Zevin et al. (2017). These glitches show a complex range of morphologies, in contrast to the glitches observed in LISA Pathfinder, which were all of a much simpler form Armano et al. (2022b). These differences most likely arise from the difference in the technology of the instruments and in their different operating environments. However, given the differences in the mission environment and duration from LISA and LPF, we should remain open to the possibility that LISA may encounter previously unobserved or more complex types of glitches. Regarding signals observed in ground-based detectors, these are typically very short in duration compared to their typical separation in time. For this reason, segments of data flagged as having glitches are normally excluded from scientific analysis. Some events, for example the first binary neutron star merger, GW170817 Abbott et al. (2017), were clearly identified as being astrophysical in origin, due to good data quality in one of the instruments, while having overlapping glitches in one or more of the other detectors. In those cases, the analysis either gated out the range of times when the glitch was present or fitted out the glitch using a flexible model Pankow et al. (2018). In the case of LISA, signal durations will be typically much longer than the average time between glitches, and so only the latter strategies will be available to us.

In this paper, we present a possible solution to address these types of artefacts—particularly force glitches acting on the test masses—within the framework of the Global Fit development for LISA. Our analysis is performed by using the open-source sampler Eryn taking advantage of its reversible jump (RJ) capabilities Karnesis et al. (2023); Katz and Roberts (2023). We demonstrate the feasibility of performing both searches and parameter estimation to fit multiple glitches of different shape, morphology and SNR along with instrumental noise and an astrophysical source, in this case an MBHB. With this methodology, we aim to contribute to the ongoing efforts to develop robust data analysis techniques for LISA, ensuring its capacity to distinguish between astrophysical signals and instrumental artefacts. We therefore show the following achievements:

1.

Effective search, parameter estimation and fit of multiple glitches while recovering the LISA instrumental noise.
2.

Methodology to distinguish between a glitch and an MBHB signal.
3.

Preservation of parameter estimation accuracy for an MBHB signal in noisy data in the presence of one or multiple glitches, injected at various times—either before, after, or during the time of maximum overlap between the MBHB signal and the glitch.

This paper is organised as follows: in section II we describe the LPF main results in terms of measured and cataloged artefacts; in section III we introduce the models used in this paper for the instrumental noise, glitches, and the MBHB signal; in section IV we provide details on the Bayesian formalism used, including the likelihood and the sampling methodology; in section V, we present an example where we generate our own data containing only glitches and instrumental noise. This demonstrates the algorithm’s ability to identify and fit multiple artefacts with different SNRs at the same time, while distinguishing them from stationary instrumental noise. Then, in section VI, we analyze the impact of glitches with different SNRs on MBHB parameter estimation when left unmodeled and assess the recovery of both glitch and MBHB parameters when performing a joint fit. Finally, in section VII, we apply our methodology to a “Light-Spritz” Data Challenge LISA LDC Team (nd). We address it as “Light-Spritz” because in this analysis we ignore the non-stationarities caused by gaps and Galactic binaries present in the dataset as these are subjects of a follow up study. Additionally, we replace the MBHB signal included in this dataset with one generated using the PhenomHM London et al. (2018) waveform model from the BBHx software package Katz et al. (2023), in order to eliminate systematics arising from differences between the simulated signal and the model used for the fit. We conclude the paper in section IX with a summary of our main results and perspectives for future work.

II Glitches as instrumental artefacts

LPF successfully demonstrated the feasibility of maintaining free-falling test masses with residual acceleration noise below the stringent requirements for LISA Armano et al. (2018). Operating from 2016 to 2017, LPF achieved unprecedented noise levels, however it also showed that measured(-experimental) noise could deviate quite significantly with respect to the derived noise models. In particular, the measured test-mass (TM) acceleration noise at frequencies below $10^{-4}$ is still not fully understood Armano et al. (2018, 2024b). Also, transient acceleration glitches were observed during the mission Sala (2023); Castelli (2020). They spanned a wide amplitude range - transferring impulses from femto-Newton seconds (fN·s) to nano-Newton seconds (nN·s) to the test-masses - and exhibited diverse durations, from seconds to hours. These glitches fell into two main categories: rapid transients in the interferometric readout and longer-lasting forces acting on the test masses. Glitches in the second category arose from spurious impulsive forces of unknown origin acting on the test masses, which occurred approximately once per day, following a Poisson distribution. These impulses exhibited a distinctive time evolution, with a peak followed by an exponential decay over several hours Armano et al. (2022a). A detailed analysis of these phenomena showed that the observed signals could be accurately described using a fitted model Sala (2023), with minimal residual deviations. The strength of these impulsive events can be quantified through their SNRs. Out of the 432 events detected, 98 occurred during ordinary runs, characterized by standard operational conditions (22 or 11 degrees Celsius). The SNR distribution of glitches that occurred during ordinary runs peaked around moderate values of 10. In contrast, more glitches (334) were observed during the cold runs performed in May 2017, during which the instrument was operated at about 0 degree. All glitches observed during the cold runs were short and intense, and had higher SNR values (see Figure 7.12 in Sala (2023)). We refer the reader to Sala (2023) for a more detailed discussion of the possible origin of this phenomenon.

In this paper, we focus primarily on glitches with morphology similar to those observed in the ordinary runs since LISA will operate at standard operational conditions. These exhibit a well-defined and predictable pattern. Indeed, the analysis of LPF glitches leveraged phenomenological modelling, as single or double exponentials with a certain decay time and amplitude Sala (2023). Recently, Baghi et al. (2022) has shown that also shapelet functions could be used to fit LPF glitches and therefore, characterise their properties. These models, even if not directly connected to physical models as the ones derived in Sala (2023); Castelli (2020), laid the groundwork for future efforts to mitigate noise artefacts in LISA that are agnostic to the physical properties of the signals, such as the glitch decay time and impulse(-amplitude), which determine their shape.

To conclude this section, we note that one of the key challenges in GW data analysis is to develop methods that can robustly separate astrophysical signals from instrumental artefacts. To address this, we first need to derive a fiducial glitch model to use in our analysis. Therefore, in the next section, we present the shapelet function used to model the glitches and we relate it to the physical glitch model used to fit the observed LPF glitches. We also describe the waveform model used to describe the MBHB signal and the Time-Delay Interferometry (TDI) technique, which is applied to the raw data as a pre-processing step to reduce laser frequency noise before parameter estimation can be performed.

III Time delay interferometry, glitches, Massive black hole binaries and noise models

III.1 Time-delay interferometry

TDI is used in LISA to suppress laser frequency noise, which would otherwise exceed the GW signal by eight orders of magnitude, due to the un-equality of the interferometer arms. The LISA mission comprises three identical spacecraft arranged in a triangular configuration, separated by 2.5 million kilometres and linked by six active laser connections. While the LISA satellite formation is designed to remain roughly equilateral, gravitational and orbital dynamics cause the arm lengths to fluctuate by approximately $1.5\%$ annually, with relative velocity drifts between spacecraft reaching up to 10 m/s Amaro-Seoane et al. (2017); Colpi et al. (2024). These time-varying and unequal arm lengths introduce laser frequency noise into the measurement channels. TDI achieves laser frequency noise reduction by forming specific linear combinations of time-delayed inter-spacecraft phase measurements Shaddock et al. (2004); Armstrong et al. (1999). Many combinations can be constructed Muratore et al. (2020), but for this work, we use the Michelson TDI combinations denoted as $X,Y$ , and $Z$ ; also because these are the data streams available in “Spritz”. Generally, when we refer to first-generation TDI Shaddock et al. (2004), this means measurement combinations that are constructed to cancel laser frequency noise for a static, but unequal arm length, constellation. In contrast, for second-generation TDI, time-varying and unequal arm lengths are taken into account Shaddock et al. (2003); Vallisneri (2005). In this paper, while we consider both generations, we assume in both cases, equal and non-varying arm lengths. Specifically, the first-generation TDI, is used for generating, and analyzing, frequency data and modelling the noise and glitches in section V. Whereas, the second-generation TDI is used here for modelling the noise, signal, and glitches in the analysis of the “Light--Spritz” Data Challenge LISA LDC Team (nd) in section VII.

Note that by assuming equal arm lengths, uncorrelated noise, and identical optical metrology system (OMS) and TM acceleration noise terms, the combinations $X,Y$ , and $Z$ can be further transformed into three noise-orthogonal channels $A,E$ , and $T$ . This orthogonality simplifies the likelihood computation - in the frequency domain - by not having cross-correlation terms between channels. The channels $A$ , $E$ and $T$ Hartwig and Muratore (2022a); Prince et al. (2002) that will be used here are:

\displaystyle{\rm A}

\displaystyle=\frac{{\rm Z}-{\rm X}}{\sqrt{2}}\;,\hskip 3.0pt{\rm E}=\frac{{% \rm X}-2{\rm Y}+{\rm Z}}{\sqrt{6}}\;,\hskip 3.0pt{\rm T}

\displaystyle=\frac{{\rm X}+2{\rm Y}+{\rm Z}}{\sqrt{3}}.

(1)

In the case of first generation TDI, the channel $X$ is defined as:

	$\displaystyle{\rm X_{1}}$	$\displaystyle=(1-D_{12}D_{21})\times(\eta_{13}+D_{13}\eta_{31})$
		$\displaystyle\qquad-(1-D_{13}D_{31})\times(\eta_{12}+D_{12}\eta_{21})\;,$		(2)

and the same variable in second generation is defined as:

$\displaystyle{\rm X_{2}}$	$\displaystyle=(1-D_{12}D_{21}-D_{12}D_{21}D_{13}D_{31}$
	$\displaystyle\qquad\quad+D_{13}D_{31}D_{12}D_{21}D_{12}D_{21})\times(\eta_{13}% +D_{13}\eta_{31})$
	$\displaystyle\qquad-(1-D_{13}D_{31}-D_{13}D_{31}D_{12}D_{21}$
	$\displaystyle\qquad\quad+D_{12}D_{21}D_{13}D_{31}D_{13}D_{31})\times(\eta_{12}% +D_{12}\eta_{21})$
	.	(3)

The $Y$ and $Z$ variables differ from $X$ by cyclic permutations of the three satellites. The operator $D_{ij}$ is a delay operator, corresponding to a constant time shift by $L_{ij}$ and thus in frequency to multiplication by $\mathcal{F}\{D_{ij}\}=e^{-i2\pi fL_{ij}}$ with $L_{ij}$ representing the light travel time for a beam received at spacecraft $i$ and emitted from spacecraft $j$ , in seconds Heinzel et al. (2011). In this work we will assume the six arm-lengths to be equal and constant and therefore $L_{ij}=L$ . The ${\eta}_{ij}$ are the so-called intermediary variables Hartwig (2021). They represent a synthesized interferometric test-mass to test-mass measurement, here expressed as a fractional frequency shift between the interfering laser beams. The first index $i$ indicates the spacecraft where the measurement is performed at time $t$ , and the second index $j$ indicates the distant spacecraft from which light was emitted at time $t-L_{ij}$ . The intermediary variables ${\eta}_{ij}$ contain the sum of the contributions from GW signals, glitches and instrumental noises, and so the TDI transfer function applies linearly to these different components.

III.2 Shapelet model

In this section, we describe how glitches are modeled and the LISA response to them. We model glitches using the exponential shapelet functions following Baghi et al. (2022) and Bergé et al. (2019). This is a family of exponentially damped functions that are eigen-wavefunctions of the normalized 1D-hydrogen atom. They can be defined in time domain as follows:

\psi_{n}(t)=2\frac{t}{n}e^{-\frac{t}{n}}L^{1}_{n-1}\bigg{(}2\frac{t}{n}\bigg{)% }\Theta(t),

(4)

with $L_{n}^{\alpha}$ denoting the generalized Laguerre polynomials. A glitch perturbation is then modeled by a linear combination of shapelets:

g(t)=\sum_{i=0}^{P}A_{i}\psi_{n_{i}}\bigg{(}\frac{t-\tau_{i}}{\beta_{i}}\bigg{% )},

(5)

with $P+1$ being the number of shapelets. The parameters to be estimated for each shapelet component are the scale parameter $\beta_{i}$ , which represents the characteristic damping time, the starting (or injection) time of the shapelet $\tau_{i}$ , and the shapelet amplitude $A_{i}$ . Note that Eq. (5) is given in units of acceleration, but considering that the phasemeter in LISA outputs phase or frequency, we must convert it to the appropriate units. To achieve this, we utilize the generating function (see Wikipedia contributors (nd) for reference and Appendix A.1) and integrate Eq. (5) -after substituting Eq. (4)- to derive the fractional frequency expression for a single exponential shapelet with $n_{0}=1$ as:

\begin{split}&\frac{\Delta\nu_{g}}{\nu_{0}}(t,\tau_{0},\beta_{0},A_{0})=\\ &\frac{2A_{0}\beta_{0}}{c}\bigg{(}1-\frac{e^{\frac{-t+\tau_{0}}{\beta_{0}}}(t+% \beta_{0}-\tau_{0})}{\beta_{0}}\bigg{)}\Theta\bigg{(}\frac{t-\tau_{0}}{\beta_{% 0}}\bigg{)},\end{split}

(6)

where we have considered that for our case $\alpha=1$ Baghi et al. (2022). The shapelet model zero order depends thus on three parameters $\phi=\{\tau_{0},A_{0},\beta_{0}\}$ . We note that Eq. (6) is the same as the one derived in Baghi et al. (2022)¹¹1Private correspondence with the authors of Baghi et al. (2022) confirmed that there is a typo in Eqs. B2-B4 of their paper and the corrected expressions is consistent with what we present here.. This function is similar to a step function starting from a value of zero and increasing to:

\lim_{t\rightarrow\infty}\frac{\Delta\nu_{g}}{\nu_{0}}(t,\tau_{0},\beta_{0},A_% {0})=\frac{2A_{0}\beta_{0}}{c}.

(7)

The shapelet model with $i=0$ thus reduces to the way glitches were modeled in LPF. Although, since the most significant impact of glitches during the LPF mission was an effective overestimation of the measured LPF acceleration noise, they had been fitted in terms of acceleration (Eq. 7.1 in Sala (2023) and Appendix A.2) using the following model:

h_{g}(t)=\frac{\Delta v}{\tau^{2}}t^{\prime}e^{\frac{-t^{\prime}}{\tau}}\Theta% (t^{\prime}),\qquad\qquad t^{\prime}=t-t_{0}

(8)

which, integrated in time to get the fractional frequency change, gives:

	$\displaystyle\int_{0}^{t}$	$\displaystyle\frac{h_{g}(t^{\prime})}{c}dt^{\prime}$
		$\displaystyle\hskip 14.22636pt=\frac{\Delta v}{c}\left(\frac{e^{\frac{t_{0}-t}% {\tau}}(-t-\tau+t_{0})}{\tau}+1\right)\Theta(t-t_{0}),$		(9)

$\Delta v$ is the total transferred impulse per unit mass, $t_{0}$ is the injection time and $\tau$ is the decay time. Therefore, comparing the shapelet model (Eq. 6) with the glitch model (Eq. III.2), the shapelet amplitude $A_{0}$ can be related to the glitch impulse $\Delta v$ accordingly to $\Delta v=2A_{0}\beta_{0}$ , while $\beta_{0}$ equals the glitch $\tau$ , and $\tau_{0}$ equals the glitch injection time $t_{0}$ . These relationships will be important in sections V.2 and VIII, where we will use them to define the priors on the shapelet for fitting the glitches present in the data. Specifically, the prior for the shapelet amplitude will be chosen based on this relationship such that smaller amplitude values will be allowed compared to the nominal LPF glitches amplitude.

Refer to caption — Figure 1: Impulse versus decay time of LPF glitches at different SNRs computed considering the LISA noise power spectral densities and the three TDI channels $A$ , $E$ and $T$ .

Figure 1 shows the distribution of glitch parameters observed for LPF glitches. It is clear that relatively short glitches have the highest SNR. This is because short glitches have the majority of the power concentrated in the middle of the LISA frequency band, whereas long glitches concentrate the power at very low frequencies, where LISA is less sensitive. We note that SNRs for LISA are higher with respect to LPF (Fig. 1 versus Fig. 8.8 in Sala (2023)) for two reasons: the first is the different acceleration noise transfer functions between LISA and LPF, in particular, LISA shows a factor $4$ enhancement for the $X$ channels in the TM acceleration noise power spectral density (PSD), the second is that we compute the SNR for LISA assuming the three TDI channels ( $A$ , $E$ and $T$ ) such that the SNR squared is doubled with respect to that in a single channel²²2The TDI channel $T$ , or in general any channel null to GWs, being less sensitive to test-mass acceleration noise is also less sensitive to acceleration glitches. Therefore, their contribution to the SNR is negligible.. We show in Figure 2 two examples of shapelets with $i=0$ and $i=9$ .

To conclude, the shapelet expression (Eq. 6) can be written directly in the frequency domain as follows³³3We have tested that using the analytical frequency domain function is sufficiently equivalent to performing a numerical fast Fourier transform of the time domain data.(see Appendix A.1):

\frac{\Delta\tilde{\nu}_{g}}{\nu_{0}}(f,\tau_{0},\beta_{0},A_{0})=\frac{iA_{0}% \beta_{0}e^{-2if\pi\tau_{0}}}{cf\pi(i-2f\pi\beta_{0})^{2}}.

(10)

This expression agrees with the physical model used to fit glitches in LPF data but expressed in fractional frequency units (Eq. 7.3 in Sala (2023) and Appendix A.2).

III.3 TDI transfer function for glitches and shapelets in frequency domain

We can express the single link measurement $\tilde{\eta}_{ij}(\omega)$ , in Eqs. III.1,3 as:

\tilde{\eta}_{ij}(f)=\Delta\tilde{\nu}_{g_{ji}}{(f)}e^{-2i\pi fL}+\Delta\tilde% {\nu}_{g_{ij}}{(f)},

(11)

where the terms $\Delta\tilde{\nu}_{g_{ij}}{(f)}$ indicates the contribution of a glitch (shapelet) acting on the test mass $ij$ and we have assumed constant and equal arm-lengths, $L_{ij}=L$ . Considering the expected glitch distributions from LISA Pathfinder, it is highly unlikely that multiple glitches happen simultaneously. We therefore assume that one of the two terms ( $\Delta\nu_{g_{ij}}{(f)}$ or $\Delta\tilde{\nu}_{g_{ji}}{(f)}$ ) vanishes and that the other one contains a single glitch.

Considering now Eq. III.1, we can model the TDI applied to a glitch (shapelet) as a product between the glitch (shapelet) expression in the frequency domain (Eq. 10) and the transfer functions $\mathcal{T}_{X,Y,Z}^{g}(f)$ which take into account TDI, i.e., $\tilde{g}_{X,Y,Z}(f)=\mathcal{T}_{X,Y,Z}^{g}(f){\Delta\tilde{\nu}_{g}}(f)/\nu_% {0}$ , where $\tilde{g}_{X,Y,Z}(f)$ is the signal of a glitch (shapelet) in the TDI frequency domain data, where we have omitted the constant glitch parameters in the argument of the expression.

In this work we will always apply glitches on $TM_{12}$ . Denoting by $\mathcal{F}$ the Fourier transform operator, and applying it to the delay operators gives $\mathcal{F}{D}=e^{-2i\pi fL}$ . The transfer function for the $X,Y$ and $Z$ first generation TDI variables of a glitch on $TM_{12}$ is:


$\displaystyle\mathcal{T}_{X_{1}}^{g}(f)$	$\displaystyle=-1+e^{-8iL\pi f}$	(12a)
$\displaystyle\mathcal{T}_{Y_{1}}^{g}(f)$	$\displaystyle=4ie^{-4iL\pi f}\sin(2L\pi f)$	(12b)
$\displaystyle\mathcal{T}_{Z_{1}}^{g}(f)$	$\displaystyle=0.$	(12c)

The same type of derivation can be done for second-generation TDI using Eq. 3:


$\displaystyle\mathcal{T}_{X_{2}}^{g}(f)$	$\displaystyle=e^{-16iL\pi f}(-1+e^{8iL\pi f})^{2}$	(13a)
$\displaystyle\mathcal{T}_{Y_{2}}^{g}(f)$	$\displaystyle=8e^{-8iL\pi f}\sin(4L\pi f)\sin(2L\pi f)$	(13b)
$\displaystyle\mathcal{T}_{Z_{2}}^{g}(f)$	$\displaystyle=0.$	(13c)

Note that, in the case of multiple glitches happening at the same time in different test masses their effects sum up linearly. Expressions for the transfer functions for glitches occurring on the other test masses can be found in Appendix A.2.

III.4 Waveform model for Massive Black Hole Binaries

To model the full inspiral-merger-ringdown signal of a MBHB, we use the PhenomHM London et al. (2018) waveform model from the BBHx software package Katz et al. (2023). This is an analytic frequency-domain model for quasi-circular and spin-aligned binaries that includes harmonics beyond the dominant $(2,2)$ mode. The contribution of these harmonics increases with the asymmetry of the system (i.e., unequal mass ratios) and the magnitude of the spins. We expect it to be possible to observe large mass ratios and highly spinning MBHBs with LISA Toubiana et al. (2021); Barausse et al. (2020a). Not including higher harmonics could result in severe biases Pitte et al. (2023). The GW signal of an MBHB, for the spin-aligned case, depends on 11 parameters $\theta=\{M_{t},q,a_{1},a_{2},d_{L},\varphi_{ref},\iota,\lambda_{d},\beta_{d},% \psi,t_{ref}\}$ . They can be split as:

•
Intrinsic Parameters:
- –
  
  $m_{1}$ is the mass of the first black hole and $m_{2}$ of the second, with $m_{2}\leq m_{1}$ .
- –
  
  $a_{1},a_{2}$ : dimensionless spin magnitudes of the two black holes, aligned or anti-aligned with the orbital angular momentum.
•
Extrinsic Parameters:
- –
  
  $d_{L}$ : luminosity distance to the source, scaling the amplitude of the GW signal.
- –
  
  $\varphi_{ref}$ : reference phase at the merger.
- –
  
  $\iota$ : inclination, the angle between the orbital angular momentum and the line of sight, affecting the observed amplitude of higher-order modes.
- –
  
  $\lambda_{d},\beta_{d}$ : ecliptic longitude and latitude, describing the sky location of the source in the (LISA) detector reference frame.
- –
  
  $\psi$ : polarization angle of the waveform, affecting the projection of the signal onto the detector.
- –
  
  $t_{ref}$ : reference time marking a characteristic point in the binary coalescence, such as the merger.

All extrinsic parameters discussed here are defined relative to the LISA constellation reference frame. We also introduce the total mass $M_{t}=m_{1}+m_{2}$ and the mass ratio $q=m_{2}/m_{1}\leq 1$ .

The frequency-domain, source-frame waveform is expressed as Marsat and Baker (2018); Marsat et al. (2021):

\tilde{h}_{lm}(f)=A_{lm}(f)e^{-i\phi_{lm}(f)},

(14)

where $(l,m)$ represent the harmonic mode indices. The TDI variables are obtained as

\tilde{h}_{X,Y,Z}(f)=\sum_{lm}\mathcal{T}^{h}_{X,Y,Z}(f,t_{lm}(f))\frac{\tilde% {h}_{lm}(f)}{d_{L}}.

(15)

where $\mathcal{T}^{h}(f,t_{lm}(f))$ are the frequency-domain transfer functions detailed in Marsat and Baker (2018); Marsat et al. (2021). The transfer functions depend on the extrinsic parameters and vary both in time and frequency due to the LISA constellation’s orbital motion. They are evaluated at time $t_{lm}(f)$ , which is given by the stationary-phase approximation:

t_{lm}(f)=t_{ref}-\frac{1}{2\pi}\frac{d\phi_{lm}(f)}{df}.

(16)

The phase-meter output (LISA observables) will be therefore the sum of GWs signals, $\tilde{h}_{X,Y,Z}(f)$ , and glitches, $\tilde{g}_{X,Y,Z}(f)$ , propagated through TDI.

III.5 Noise models and definition of the signal to noise ratio

In this work, we assume that primary noise sources such as laser frequency noise, spacecraft longitudinal jitters, clock noises and tilt-to-length couplings have been corrected and only leave negligable residuals in the post-processed TDI variables Hartwig and Bayle (2021); Hartig et al. (2022); Chwalla et al. (2020). What remains are the so-called secondary noises. These remaining noises Quang Nam et al. (2023a) fall into two broad categories: residual acceleration noise of each TM and the OMS noise terms limiting the test-mass position readout. Noises inside these categories can have various physical origins, different levels and different transfer functions into the final TDI variables Quang Nam et al. (2023b). For the analysis that will follow in section V, we just consider two main components modeled as uncorrelated noise on each test-mass and uncorrelated noise in each inter-spacecraft optical link, with identical statistical properties for noises of the same kind. We will instead consider a slightly more elaborate noise model for analysing the “Light-Spritz” data set, in section VII.2.

We characterize these noise sources by their PSD and a corresponding transfer function that describes its propagation in the TDI variables. The TM acceleration noise PSD is:

\begin{split}S_{TM}(f)=(N_{TM})^{2}&\left(1+\left(\frac{0.4\times 10^{-3}}{f}% \right)^{2}\right)\\ \times&\left(1+\left(\frac{f}{8\times 10^{-3}}\right)^{4}\right),\end{split}

(17)

with $N_{TM}=$2.4\text{\,}\mathrm{f}\mathrm{m}\mathrm{/}\mathrm{s}\mathrm{{}^{2}}% \mathrm{/}\sqrt{\mathrm{Hz}}$$ . The PSD of the displacement noise from the OMS is:

S_{oms}(f)=(N_{OMS})^{2}\left(1+\left(\frac{2.0\times 10^{-3}}{f}\right)^{4}% \right),

(18)

with $N_{OMS}=$7.9\text{\,}\mathrm{p}\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$$ . This overall OMS noise level is based on the values provided in the LDC tools working group (2022).

Note that the noise parameters, $\gamma=\{N_{TM},N_{OMS}\}$ , are expressed as equivalent acceleration noise in $\mathrm{m}\mathrm{{}^{2}}\mathrm{/}\mathrm{s}^{4}\mathrm{/}\mathrm{Hz}$ or displacement noise in ${\mathrm{m}}^{2}\text{\,}{\mathrm{Hz}}^{-1}$ . For our analysis, we instead use units of fractional frequency deviations, to which we will convert them by multiplying the components with $C_{\mathrm{acc}}^{\mathrm{ffd}}=\left({1}/{(c2\pi f)}\right)^{2}$ or $C_{\mathrm{disp}}^{\mathrm{ffd}}=\left({2\pi f}/{c}\right)^{2}$ , respectively. Here, $c$ is the speed of light in $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ .

For the first generation TDI, the noise in the $A,E$ and $T$ channels are:

$\displaystyle S_{A}(f)=$	$\displaystyle S_{E}(f)=8\sin^{2}(2\pi Lf)\Big{[}$
	$\displaystyle 2C_{\mathrm{acc}}^{\mathrm{ffd}}S_{TM}(3+2\cos(2\pi Lf)+\cos(4% \pi Lf))$	(19)
	$\displaystyle+C_{\mathrm{disp}}^{\mathrm{ffd}}S_{\text{oms}}(2+\cos(2\pi Lf))% \Big{]},$

and:

	$\displaystyle S_{T}(f)=$	$\displaystyle 16C_{\mathrm{disp}}^{\mathrm{ffd}}S_{\text{oms}}(1-\cos(2\pi Lf)% )\sin^{2}(2\pi Lf)$
		$\displaystyle+128C_{\mathrm{acc}}^{\mathrm{ffd}}S_{TM}\sin^{2}(2\pi Lf)\sin^{4% }\left(\pi Lf\right).$		(20)

Where, again, $L$ denotes the LISA arm length in seconds.

III.5.1 Definition of the signal to noise ratio

The SNRs of glitches and MBHB signals are computed based on the three TDI channels using the following formula (Eq. 119 in Tinto and Dhurandhar (2021) and Babak et al. (2021)):

	$\displaystyle SNR=$	$\displaystyle\Big{(}4\Re\Big{[}\int_{0}^{f_{\text{max}}}\|\tilde{Q}(f)\|^{2}\Big% {(}\frac{\|\mathcal{\tilde{T}}_{A}(f)\|^{2}}{S_{A}(f)}+\frac{\|\mathcal{\tilde{T}% }_{E}(f)\|^{2}}{S_{E}(f)}$
		$\displaystyle\quad+\frac{\|\mathcal{\tilde{T}}_{T}(f)\|^{2}}{S_{T}(f)}\Big{)}\,% \mathrm{d}f\Big{]}\Big{)}^{1/2}.$		(21)

Note that the above formula is valid in this case as we assume the three TDI channels to be uncorrelated. Here $\mathcal{\tilde{T}}_{A/E/T}(f)$ are the transfer functions of a GW signal or a glitch into the combinations $A$ , $E$ , and $T$ and $\tilde{Q}(f)$ is the amplitude of the GW signal ( $\tilde{h}_{lm}(f)/d_{L}$ ) or glitch ( $\Delta\tilde{\nu}_{g}(f)/\nu_{0}$ ) in Fourier frequency domain, and $S_{A/E/T}(f)$ are the spectral densities of the noises in the $A$ , $E$ , $T$ combinations respectively. All these quantities have been defined in the previous sections.

IV Bayesian methodology and likelihood formulation

The Bayesian framework provides a powerful approach to the problem we aim to solve by treating the detection and parameter estimation of GW signals, instrumental artefacts and noise as an inference problem. In this context, the observed data, $d(t)$ , is modelled as the sum of contributions from the GW signal, $h(t,\theta)$ , the glitch, $g(t,\phi)$ , and the noise, $n(t,\gamma)$ :

d(t,\theta,\phi,\lambda)=h(t,\theta)+g(t,\phi)+n(t,\gamma),

(22)

where $\theta,\phi$ and $\gamma$ represent the parameters of the GW signal, the glitch and the noise model, respectively. Denoting by $\mu=(\theta,\phi,\gamma)$ the set of all parameters, we write $p(\mu|d,M)$ the posterior distribution of the parameters. We note $M$ the assumed model, which includes $M_{1}$ , our chosen shapelet model for the glitches, $M_{2}$ , our model for the MBHB signal, and $M_{3}$ , our model for the instrumental noise. The posterior distribution $p(\mu|d,M$ ) is computed using Bayes’ theorem Sivia and Skilling (2006); Karnesis et al. (2023):

p(\mu|d,M)=\frac{p(d|\mu,M)\,p(\mu|M)}{p(d|M)}.

(23)

where:

•

$p(d|\mu,M)$ is the likelihood, quantifying the probability of observing data $d$ under the models $M_{1}$ , $M_{2}$ and $M_{3}$ with parameters $\theta,\phi$ , and $\gamma$ ;
•

$p(\mu|M)$ is the prior distribution for the GW, glitch and noise parameters under the models $M_{1}$ , $M_{2}$ and $M_{3}$ : $p(\theta|M_{1})$ , $p(\phi|M_{2})$ and $p(\gamma|M_{3})$ ;
•

$p(d|M)$ is the evidence, which acts as a normalization constant for parameter estimation, but can also be used to compare different models $M$ .

For the problem at hand, we will use the Whittle likelihood, for which the log-likelihood in the frequency domain is:

	$\displaystyle\log\mathcal{L}=-\frac{1}{2}\sum_{f}\Bigg{[}$	$\displaystyle 4\Delta f\frac{\left\|\tilde{d}(f)-\tilde{h}(f,\theta)-\tilde{g}(% f,\phi)\right\|^{2}}{S(f,\gamma)}$
		$\displaystyle+\log S(f,\gamma)\Bigg{]},$		(24)

where $\tilde{d}(f),\tilde{h}(f,\theta)$ and $\tilde{g}(f,\phi)$ are the Fourier transforms of the same quantities appearing in Eq. 22; $S(f,\gamma)$ is the one-sided noise PSD, parameterized by $\gamma$ ; and $\Delta f$ is the frequency resolution which is equal to $1/T_{\rm obs}=1/(Ndt)$ , with $T_{\rm obs}$ the observation time, $N$ the length of the time series and $dt$ the inverse of the sampling rate.

The posterior can be sampled using advanced Bayesian inference techniques such as Markov Chain Monte Carlo (MCMC) Metropolis et al. (1953); Hastings (1970); Hitchcock (2003); Gilks et al. (1995). In this analysis, we wish to fit simultaneously an unknown number of glitches. Therefore, we use a Reversible Jump MCMC (RJ-MCMC) Green (1995) algorithm, a generalization of the Metropolis–Hastings algorithm that allows to vary the dimensionsality of the parameter space while sampling over the parameters. We use the RJ-MCMC implementation of Eryn Karnesis et al. (2023). Several walkers are run in parallel and exchange information periodically in order to increase the sampling efficiency. We also use parallel tempering Swendsen and Wang (1986); Hukushima and Nemoto (1996); Vousden et al. (2015) to improve the exploration of multi-modalities in the posterior distribution. During the run, a new point $\mu^{\prime}$ is proposed from a point $\mu$ according to a proposal function, $q(\mu^{\prime}|\mu)$ , and is accepted with probability:

\alpha=\min\left(1,\frac{p(\mu^{\prime}|d,M)\,q(\mu|\mu^{\prime})}{p(\mu|d,M)% \,q(\mu^{\prime}|\mu)}\right).

(25)

In our analysis, we use different proposals for the parameters of the MBHB ( $\theta$ ), of the intrumental noise ( $\gamma$ ) and of the glitches ( $\phi$ ). They are updated successively, using the Gibbs sampling implementation of Eryn Karnesis et al. (2023). In the Eryn ‘language’, branches represent different types of models used to fit the data. Each branch contains leaves, which correspond to individual instances of each model. In our case, leaves represent individual glitches, massive black hole binaries, or noise components, which form different branches Karnesis et al. (2023).

IV.1 Moves

In an RJ-MCMC sampler, there are two types of moves:

•

In-model moves: changes are made to the parameters, keeping the number of leaves in each branch fixed.
•

Out-of-model moves: the dimensionality is varied by adding or removing leaves, in one or several branches. For instance, in our case, adding or removing glitches.

In this work, only the number of glitches is varied, so that the latter type of moves applies only to their parameters. Detailed explanations of the various moves available in Eryn, as well as the most commonly used ones in the literature, can be found in Katz ; Karnesis et al. (2023). Here, we focus on the in-model moves specifically developed for glitch searches, MBHB parameter estimation, and noise parameter estimation. For the out-of-model move for glitches, we propose new parameters from the prior distribution during the search phase. Instead, during the final parameter estimation phase, as explained in the upcoming section, we propose sampling from a Gaussian distribution, whose covariance matrix is constructed from an initial estimate of the parameters of the detected glitches during the search phase.

IV.1.1 Glitch search and parameter estimation (in-model moves)

Group stretch move

In the stretch move new samples are generated by leveraging the distribution of points in the ensemble, using the formula:

X_{k}=X_{j}+z(X_{k}-X_{j}),

(26)

where $X_{k}$ is the current sample, $X_{j}$ is another point in the ensemble, and $z$ is a scaling factor drawn from a predefined distribution. In RJ-MCMC, however, the variable dimensionality of models complicates the direct use of ensemble-based proposals. To deal with the variable dimensionality of models the Group Stretch Move is used Karnesis et al. (2023) which builds upon the affine-invariant stretch proposal Goodman and Weare (2010). This group introduces ‘friends’, a static subset of reference points that are defined after a first burn of the sampler and as soon as the group move is used. These friends simulate the posterior distribution for a given model and remain fixed throughout multiple iterations. For instance, in our case, we group the estimated glitches into ‘friend groups’ according to their mean time value. The goal is to set $X_{j}$ properly for each posterior mode such that $X_{j}$ and $X_{k}$ are always found on the same mode of the posterior. This approach improves efficiency by ensuring that proposed points correspond to the posterior distribution of the currently sampled glitch, rather than a different one.
The stretch move is used 50 $\%$ of the time, and in the implementation of this move, the set of ‘friends’ associated with each glitch is tracked, to improve the efficiency of the RJ-MCMC to add and remove glitches Karnesis et al. (2023); Katz and Roberts (2023).

Selected Covariance Group Move

The ‘Selected Covariance Group Move’ replaces the standard Gaussian move. In the context of nested models, the posterior distribution may display distinct modes, corresponding to different glitches in the data. This complexity makes it challenging to efficiently tune covariance matrices for each mode. The core idea behind this move is to group glitches based on their mean injection time values and compute corresponding distributions for glitch parameters using a Gaussian Mixture Model (GMM). This approach leverages covariance matrices specific to each glitch group. When glitches are clearly distinguishable, the group proposal becomes highly efficient, aided by the branch supplemental object in Eryn Katz et al. (2023) to maintain the association between glitches and their respective covariance matrices⁴⁴4The branch supplemental object tracks the friends’ information associated with each glitch. This is particularly needed in the presence of RJ-MCMC and parallel tempering, where we have swaps in temperature or birth/death of new glitches.. Note that we leverage the “fix friends” function from Eryn when a proposed birth of a new glitch is accepted. This function assigns a covariance matrix (or friends) to the new glitch from the precomputed friends group.

The computation(-update) of the covariance matrices for this move is done as follows:

1.

Cluster specification: every 500 iterations, the fraction of samples with zero leaves is computed with respect to the total number of leaves Karnesis et al. (2023). If this fraction is less than 0.5, then a GMM is initialised with several modes equal to the maximum allowed number of glitch leaves, allowing it to model complex distributions as a combination of multiple Gaussian distributions. We assume that each component has a distinct, full-rank covariance matrix. So each distinct glitch-cluster estimated from the samples corresponds to a potential glitch, with the total number of possible clusters constrained by the maximum number of glitches allowed by the problem’s dimensionality.
2.

Covariance Calculation: based on the identified number of clusters, a GMM is fitted to the data. The GMM provides the mean and covariance for each cluster, capturing the local structure of the posterior around each glitch⁵⁵5In a GMM, each Gaussian component does not necessarily correspond to an independent cluster. If a posterior distribution is not very Gaussian in shape, more than one component may be needed to approximate it. Therefore, identifying a GMM component is not the same as identifying a true cluster represented by the posterior of a single glitch. However, for the problem analyzed here, we have verified that the glitch posteriors are sufficiently Gaussian and well-separated, so our assumptions remain valid. Dividing a true cluster into more sub-clusters than are present in the data leads to proposal steps that are smaller than necessary. This can result in a high acceptance rate but a low effective sample size, meaning the sampler requires more steps to converge. To improve sampling efficiency, this proposal method should be refined in future work. One development we have already begun testing is using x-means instead of k-means to automatically determine the optimal number of GMM components. . These parameters are then utilized in the in-model move, enabling it to propose moves that are well-aligned with the posterior’s structure. This approach ensures the sampler’s efficiency, particularly when the shape of the posterior varies significantly between glitches, warranting the use of distinct covariance matrices for each glitch. By iteratively estimating the covariance matrices from the samples via a GMM, and adapting the sampler to the posterior structure, this methodology establishes a robust framework to be used in case of glitch searches. Note that, adaptation is applied for the first 50000 iterations and then stopped to preserve detailed balance. The covariance group move is used 50 $\%$ of the time.

IV.1.2 Glitch search and parameter estimation (out-of model moves)

During the parameter estimation phase, we incorporate a move proposal mechanism based on detected glitches. Specifically, $70\%$ of the time, the glitch parameters are drawn from a mixture of independent Gaussian distributions describing the shapelet parameter $\tau_{0}$ , $ln(A_{0})$ and $\beta_{0}$ . Each parameter is sampled from a multivariate normal distribution with the corresponding mean and covariance of the previously detected shapelets. We propose from the prior distribution the other $30\%$ of the time.

IV.1.3 Noise parameter estimation (in-model moves)

We consider both the stretch and Gaussian proposal moves Katz . Initially, the Gaussian move uses a diagonal covariance matrix with small entries (following a suggestion from the Eryn tutorial). Then, the covariance matrix is re-estimated using the noise samples every 500 iterations and the adaptation is done for the first 50000 iterations to preserve detailed balance. The stretch move is used 85 $\%$ of the time while the Gaussian move 15 $\%$ .

IV.1.4 Massive black hole binaries parameter estimation (in-model moves)

Following Ref. Katz (2022), we use a stretch proposal described above, which is used for approximately $80\%$ of proposals, and a sky-mode hopping proposal. We also add a Gaussian move proposal used $5\%$ of the time with a small covariance matrix. This matrix is then recursively estimated from the MBHB samples with the same cadence used to estimate the noise covariance matrix. We notice that this additional move helps the convergence of the algorithm in possible cases where $X_{k}$ and $X_{j}$ (Eq. 26) are very close in value as they will spread out slowly otherwise if we only use the stretch move. Instead, the covariance Gaussian move will speed this up.

Regarding the sky-mode proposal, when sampling in the LISA reference frame Marsat et al. (2021), as we do here, there exist eight distinct sky modes in terms of the extrinsic parameters. These consist of four longitudinal modes at:

\{\lambda_{d}+(0,1,2,3)\pi/2,\psi+(0,1,2,3)\pi/2\}

and two latitudinal modes:

\{\pm\beta_{d},\pm\cos\iota,\pm\cos\psi\}.

Three varieties of this proposal are used:

•

One is to move the walker to a mode drawn randomly with replacement from all eight sky modes. This is used for approximately $4\%$ of proposals.
•

Proposals to change just the longitudinal or latitudinal mode are used approximately $3\%$ and $8\%$ of the time, respectively.

The relative percentages of the sky-mode hopping proposals reflect that, at higher frequencies, the eight modes reduce to just the two latitudinal modes at the proper ecliptic longitude. In this situation, longitudinal mode-hopping proposals are not efficient.

V Testing the reversible jump with simulated data

To test our pipeline, we generate mock data by injecting glitches—extracted from the LPF catalogue Armano et al. (2022b) and shown in Fig. 1—into a noise-only dataset. For this test, we do not include any astrophysical signals, as our goal is to investigate the ability of the RJ-MCMC to identify glitches of varying shapes and SNRs, and to distinguish them from noise in order to establish an approximate threshold below which such glitches are consistent with instrumental noise. In Section VI, we also quantify the impact of three fiducial glitches—with low, medium, and high SNRs of 21, 72, and 1544, respectively—on the parameter estimation of an MBHB.

V.1 Data generation

We generate 30 days of LISA data where seven glitches are injected according to the amplitudes and damping times extracted from the LPF catalogue Armano et al. (2022b), and are randomly distributed over these 30 days. In LPF the glitches were observed to have a Poisson distribution with about one glitch per test mass per day. We use a slightly lower rate here to simplify the computational burden as we do not believe this will affect our results qualitatively. Moreover, we allow also for glitches injected very close in time to show the ability of the methodology to distinguish between such glitches. Table 1 reports the parameters of the injected glitches as well as the corresponding number indicating the position of the glitch in the LPF catalogue (the run index is in accordance with LPF numbering).
We model the injected glitches using the single exponential model reported in section III Eq. III.2, but we use shapelets with $n_{0}=1$ to fit for them.

We simulate LISA noise directly in the frequency domain following Muratore et al. (2024). Assuming that the real-time series $n(t)$ follows a stationary, zero-mean, Gaussian and ergodic random process, the Fourier transform of the noise, $\tilde{n}(f)$ is normally distributed with zero-mean and variance given by the one-sided power spectral density:


$\displaystyle<\tilde{n}^{*}(f^{\prime})\tilde{n}(f)>$	$\displaystyle=\frac{1}{2}S_{n}(f)\delta(f-f^{\prime}),$	(27a)
$\displaystyle<\tilde{n}(f^{\prime})\tilde{n}(f)>$	$\displaystyle=0,$	(27b)

where $<>$ denotes the expectation value over the data-generating process. For a finite grid of frequencies, defining $n_{k}=\tilde{n}(f_{k})$ , Eq. (27a) can be written as

\displaystyle<\tilde{n}^{*}_{k}\tilde{n}_{j}>

\displaystyle=\frac{T}{2}S_{n}(f_{k})\delta_{jk}.

(28)

Therefore the variance of the real and imaginary parts of the noise is given by (where the real and imaginary part are uncorrelated):

<\Re[\tilde{n}_{k}]^{2}>=<\Im[\tilde{n}_{k}]^{2}>=\frac{T}{4}S_{n}(f_{k})=% \sigma_{k}^{2}\,.

(29)

The noise in the Fourier domain is then generated as:

\tilde{n}_{k}=\sigma_{k}\cdot\left(X_{\text{real}}+iX_{\text{imag}}\right),

(30)

where $X_{\text{real}}\sim\mathcal{N}(0,1)$ and $X_{\text{imag}}\sim\mathcal{N}(0,1)$ are independent normally distributed random variables with mean 0 and variance 1.

For this example the noise PSD is given by the sum of the TM acceleration and OMS noises which follow the models described in section III.5. The assumed value for the TM acceleration noise is $N_{TM}^{2}=\left(2.4\times 10^{-15}~{}\textrm{m}\,\textrm{s}^{-2}\textrm{Hz}^{% -1/2}\right)^{2}$ and for the OMS noise $N_{OMS}^{2}=\left(7.9~{}\textrm{pm}/\sqrt{\textrm{Hz}}\right)^{2}$ .

Figure 3 shows the Amplitude Spectral Density (ASD) in the TDI channels resulting from the injected glitches. The channels $A$ and $E$ show an excess of power with respect to the reference noise model in nominal conditions, shown in black. This excess of power is mild in the null channel $T$ , as anticipated in Muratore (2021). The reason is that for null channels, acceleration-glitches have a mild effect due to the OMS and TM acceleration noise transfer functions of these channels. In particular, the latter shows a suppression at low frequencies with respect to the OMS such that the OMS noise dominates for these channels at all frequencies Muratore et al. (2023).

Glitch	$t_{0}$ $[s]$	$\ln(\Delta v/[ms^{-2}])$	$\tau$ $[s]$	SNR
#3	$2.97\times 10^{5}$	$-25.5$	$184.47$	$72$
#118	$7.50\times 10^{5}$	$-25.5$	$4.23$	$230$
#90	$1.25\times 10^{6}$	$-24.4$	$43.21$	$479$
#72	$1.75\times 10^{6}$	$-24.5$	$1.86$	$637$
#179	$2.20\times 10^{6}$	$-23.9$	$4.77$	$1172$
#248	$2.25\times 10^{6}$	$-23.6$	$13.96$	$1544$
#41	$2.50\times 10^{6}$	$-28.3$	$5.57$	$21$

Table 1: Parameters

t_{0}

\ln(\Delta v)

\tau

and SNR for each injected glitch. For each glitch we also indicate the corresponding number (

\#

) from the LPF catalogue.

V.2 Data analysis

We analyze a complete 30 days dataset in a single run. Using the shapelet model with $n_{0}=1$ , we fit for glitches allowing up to 10 glitches in the data, also including the possibility of no glitches. The noise shape is considered fixed, while the two amplitude parameters of the OMS and TM noise, $N_{OMS}$ and $N_{TM}$ , vary within the uniform prior range specified in Table 2. We also use uniform priors for the time of injection ( $\tau_{0}$ ) and the decay time ( $\beta_{0}$ ) of the shapelet, whereas for the amplitude ( $A_{0}$ ) of the shapelet we consider a uniform prior in natural log amplitude. In this way, we avoid sampling parameters with vastly different numerical scales—such as amplitudes on the order of $10^{-12}$ and injection times around $10^{4}\sim 10^{6}$ . This choice tends to give more weight to low-SNR glitches. The goal here is to ensure that the sampler spends enough time in this region of parameter space, where most glitches occur. Figure 1 illustrates this relationship by showing the impulse, duration, and SNRs of the glitches in the LPF catalog. The analysis takes about 2 days using 70 walkers and 10 temperatures, running on a single A100 GPU. We expect that this timing can be improved in the future.

Parameter	Lower Bound	Upper Bound
$\tau_{0}$ [s]	0	2592000
$ln(A_{0}/[ms^{-2}])$	-35	-20
$\beta_{0}$ [s]	1	1e5
$N_{TM}$ [ $ms^{-2}Hz^{-1/2}$ ]	(7e-12)	(8e-12)
$N_{OMS}$ [ $mHz^{-1/2}$ ]	(2e-15)	(3e-15)

Table 2: Table of priors used for shapelet and noise parameters.

Figure 4 reports the posterior for the parameters of the shapelets used to fit the injected glitches.

The parameters of the injected glitches (Table 1) — converted to shapelet parameters according to $\tau_{0}=t_{0}$ , $A_{0}=\frac{\Delta v}{2\beta_{0}}$ and $\beta_{0}=\tau$ — are indicated within the Figure.

It can be observed in Fig. 4 that the parameters of the glitch with the lowest SNR (21) cannot be recovered, and that the posteriors of the two glitches that are close in time overlap with each other such that in the 1D plot we observe a fifth bigger Gaussian for the time of injection $\tau_{0}$ . The failure to recover the low-SNR glitch is likely due to the fact that the search did not detect it, possibly because proposed points rarely fall within the correct region of parameter space. This issue is probably linked to the choice of a uniform prior on $\beta_{0}$ (along with its lower amplitude), which turns out to be too broad relative to the $\tau_{0}$ values for some of the injected glitches —despite this choice being consistent with the general LPF decay time distribution shown in Fig. 1. To verify this hypothesis, we repeated the analysis using a logarithmic prior on $\beta_{0}$ and we observed that the low-SNR glitch is better constrained⁶⁶6Additionally, to exclude the possibility that the low-SNR glitch is consistent with noise, we have performed another test in which we fixed the number of glitches to seven, initialized the sampler at the true parameters and allowed the sampler to adjust slightly. Once the parameters stabilized, we restarted the run with RJ turned on, allowing the sampler to remove glitches. The sampler did not remove any of the glitches once they were found, indicating that—at least under the specific settings used in this example—this low SNR glitch is not consistent with noise.. However, it is important to underline that the recovery of the two noise parameters is nonetheless unbiased, as shown in Fig. 5. Therefore, the recovered PSD is well-constrained (Fig. 6) and closely aligns with the PSD model used to simulate the data, as illustrated in Fig. 6. In the lower plot of Fig. 4, we also observe a larger number of fictitious low-amplitude glitches (“ghost glitches”) arising from random fluctuations in the stationary component of the noise. These glitches have SNRs $\lesssim 24$ . We also plot the glitches from LPF with $SNR<15$ and $15<SNR<30$ . Notably, the low-SNR glitch with $SNR=21$ , which was not fitted, lies within the region of the latter. This suggests that glitches with parameters in the same region as those glitches are consistent with noise, given our choice of prior.

Finally, Fig. 7 shows the number of estimated glitches in the data at zero temperature. Even though seven glitches were injected, the sampler prefers to have eight glitches, which means adding two “ghost glitches”, since one of the injected ones is not recovered. The key difference between the “ghost glitches” and the missing one is that the amplitudes of the former are so low that the associated sampled value of $\beta_{0}$ and $\tau_{0}$ have no impact on the likelihood, and so, if such signals have significant prior support, they will readily be added and removed by the RJ step of the sampler.
The amplitude of the missing glitch is small enough to evade recovery, but not so low that the associated $\tau_{0}$ and $\beta_{0}$ have no impact on the likelihood, thus glitches with a similar amplitude are less likely to be accepted. Therefore, only values of $\tau_{0}$ and $\beta_{0}$ in the right portion of the parameter space give high enough likelihood. To conclude, the addition of “ghost glitches” is the result of allowing for low-amplitude glitches in the prior, which can be used to fit random excursions in the instrumental noise. As these fictitious glitches do not absorb much of the power in the data, this does not impact noise estimation, but nonetheless there are techniques to mitigate this effect. One approach, used for example in Littenberg et al. (2020) in the context of including white dwarf binaries in the LISA Global Fit, is to use a stricter prior on the SNR of the white dwarf binaries that are proposed to be added.

To estimate the lowest SNR glitch that can be distinguished from noise, we considered different scenarios in which the $\#41$ glitch, the low SNR one that was not recovered in our analysis, was replaced by other glitches of increasing SNRs, starting from SNR 20 to 65. We concluded that, for the noise model used in this study, glitches with an SNR of a few tens cannot be distinguished from noise. In the next section, we show that such glitches do not introduce significant biases in the parameter estimation of MBHB. Biases start to appear from SNRs of $\sim 60$ . This observation could also be used to define a prior for glitches based on their SNRs to help mitigate the confusion noise from the low-SNR sources discussed in the preceding paragraph.

VI Effect of glitches on parameter estimation of a Massive Black Holes Binary

We consider now three different cases where we inject three glitches with SNRs of 21, 72 and 1544 respectively. The first glitch corresponds to $\#41$ , the second to $\#3$ and the third to $\#248$ of the LPF catalogue. The parameters of these glitches are reported in the Table 1. We inject one glitch at a time, substituting it with another in each subsequent example. We simulate for all the three cases a fiducial MBHB signal at redshift 4, with SNR 1152, in noisy data. The choice of the MBHB is driven by Fig. 3.5 of Colpi et al. (2024) which shows the typical SNRs and redshift of MBHBs sources expected in LISA as well as Fig. 3 in Amaro-Seoane et al. (2017).

The MBHB priors used for this analysis are given in Table 3. Uniform distributions are used on parameters $\{ln(M_{t}),q,a_{1},a_{2},d_{L},\varphi_{ref},\lambda_{d},\psi,t_{ref}\}$ . The prior on total mass is natural log-uniform. For the inclination, we use a uniform prior in $\cos\iota$ between -1 and 1, while the prior on the ecliptic latitude is uniform in $\sin\beta_{d}$ .

Parameter	True	Lower Bound	Upper Bound
$ln(M_{t}/M_{\odot})$	14.48	ln( $10^{5}$ )	ln( $10^{8}$ )
$q$	0.4628	0.01	0.999
$a_{1}$	0.7474	-1	1
$a_{2}$	0.8388	-1	1
$d_{L}[Gpc]$	37	0.01	$1000$
$\varphi_{ref}$ $[rad]$	1.194	0	2 $\pi$
$\cos\iota$	0.5	-1	1
$\lambda_{d}$ $[rad]$	1.288	0	2 $\pi$
$\sin\beta_{d}$	-0.298	-1	1
$\psi$ $[rad]$	$\pi/6$	0	$\pi$
$t_{ref}[s]$	$2.628\times 10^{6}$	$t_{ref}$ - 100	$t_{ref}$ + 100

Table 3: True values of the MBHB parameters together with the prior range used for parameter estimations. We assume a flat prior on the specified range. Values of the masses are given in the detector frame.

We choose the injection time for each of the three glitches to maximize the impact of the glitch on the MBHB parameter estimation, following the procedure described in Spadaro et al. (2023). The idea is to maximize the overlap between the glitch and the MBHB signal, defined as:

\mathcal{O}(h_{1},h_{2})=\frac{(h_{1}|h_{2})}{\sqrt{(h_{1}|h_{1})(h_{2}|h_{2})% }}.

(31)

Here, $h_{1}$ stands for the MBHB signal and $h_{2}$ for the glitch. Defining $\Delta t$ the difference between the MBHB merger time and the glitch injection time, we found that the best match for the glitch with SNR 21 is at $|\Delta t|=6.80361$ minutes before the merger, for the glitch with SNR 72 it is at $|\Delta t|=25.2338$ minutes before the merger, and for the glitch with SNR 1544 at $|\Delta t|=11.5665$ minutes before the merger. Figure 8 shows the variation of the overlap versus time for the three glitches, where we have fixed the amplitudes and decay time of the glitches while varying the injection time. The dotted vertical lines indicate the injection times, with respect to the merger, at which the overlap is maximized for each of the three glitches. We show in Fig. 9 the three glitches injected at the maximum overlap time and the fiducial MBHB.

Note that the frequency of the MBHB evolves over time, and the maximum impact of a glitch is expected when its frequency overlaps with the instantaneous frequency of the MBHB signal. Acceleration glitches typically occur at frequencies on the order of $10^{-3}$ Hz, which aligns more closely with the early inspiral phase. However, if the MBHB merger occurred at lower frequencies—as in the case studied in Spadaro et al. (2023)—the glitch would instead have the strongest overlap during the merger phase.

The posteriors shown in Fig. 10 and Fig. 11 illustrate the effect of the presence of these different glitches if we attempt to recover the parameters of this MBHB without simultaneously fitting for both the artefact and the GW source. For comparison, the posteriors of the MBHB parameters obtained from data containing no glitches are reported in both plots, and the true values of the parameters are indicated by dashed black line. In the case of the SNR 21 glitch, shown in Fig 10, the glitch has less impact in the intrinsic parameters with respect to the impact we see in the extrinsic ones, as the shift in the posterior is less than the width of the posterior. The bias starts to be significant for the glitch with SNR 72 for several of the parameters, e.g. the total mass, the reference phase at merger and the ecliptic longitude and latitude. For the SNR 1544 glitch shown in Fig. 11, the two posteriors are completely disjoint, with all parameters significantly biased. To help the reader visualise the effect of not fitting such a glitch, we report only the posteriors for the intrinsic parameters. It is important to specify that the noise model has been fixed for these analyses and we expect that variations in the noise model would allow to partially absorb the impact of these glitches making the MBHB posterior more consistent with the injected values.

Let us recall that the injection times of the glitches were chosen specifically to maximize their impact on the parameter estimation of the MBHB. Their effect is expected to be smaller if they occur during the early inspiral. This point is further discussed in Appendix B, where we present posteriors illustrating how the presence of the SNR 1544 glitch affects parameter estimation for the same MBHB, depending on its injection time. However, since the coalescence is the loudest part of MBHB signals and largely drives parameter estimation, failing to captures glitches occuring near to it would severely jeopardize the science objectives of LISA. For instance, biases in sky localization and distance—such as those shown in Fig. 11 when not accounting for the glitch—could prevent the detection of an electromagnetic counterpart and lead to misidentification of the host galaxy, thereby compromising cosmological measurements based on bright or dark siren methods Gair et al. (2023); Caprini and Tamanini (2016); Tamanini et al. (2016); Mangiagli et al. (2023). This issue also extends to applications in fundamental physics, which rely on detecting extremely small deviations from general relativity. In this context, even small biases can lead to false detections and must be mitigated as much as possible Kwok et al. (2022); Gupta et al. (2025). Our analysis shows that, to avoid such pitfalls, glitches with SNRs of several tens or higher will need to be included in the model fitted to the data. This threshold also corresponds to the point at which glitches first become identifiable in the data, as discussed in the previous paragraph.

In Fig. 12, we now show the posterior on the parameters of the MBHB signal obtained when simultaneously fitting for the glitch, and compare this to the posterior obtained in the absence of the glitch. We observe that by fitting the glitch, the MBHB posteriors are no longer biased, and the posteriors for the (only-)MBHB case overlap with those where the glitch has been accounted for, indicating that the best strategy is to fit for both. For completeness, we also show posteriors for the fitted glitch parameters in Fig. 13.

To conclude the section, we also investigated the effect of injecting glitches ( $\#41$ ) and ( $\#3$ ) into MBHB with a SNR of 320, in order to reproduce the example studied in Spadaro et al. (2023). The posteriors plot are reported in the Appendix C. For the short-duration glitch with SNR 21 ( $\#41$ ), we find results consistent with those reported there: the inferred parameters remain largely unbiased. In contrast, glitch ( $\#3$ ), with an SNR of 72, though also short in duration, introduces a more significant bias.

VII “(Light-)Spritz” data set

In this section we describe the “Spritz” data set and the strategy we have used to analyze it. We describe the mock data in section VII.1 and the noise model used for this analysis in section VII.2, since it differs from the noise that was used to generate and analyse the data in the previous section VI. Finally, in Section VIII, we present the methodology used for the initial search and parameter estimation of glitches, noise, and MBHBs, followed by the subsequent analysis of the full dataset to refine the posterior estimates of the noise, gravitational wave signal, and artifacts.

VII.1 Introduction

Full details of the “Spritz” data-sets can be found in LDC (2022). The “Spritz” data release comprised three main data sets, of which we focus on the first one, which contains a single loud MBHB and three glitches as visible in Fig. 14 which shows the dataset in the time-frequency domain, and Fig. 15 that shows it in the frequency domain. The injected glitches in “Spritz” are referred to as double decay exponential glitches Sala (2023). All the data sets are realized with TDI $X,Y$ and $Z$ second generation assuming Keplerian orbits.

The dataset also features more realistic instrumental noise (see section VII.2), a background of Galactic binaries and data gaps. Auxiliary data included in the dataset allows individual components to be perfectly subtracted, allowing analysts to focus on selected aspects of the dataset. Since we focus on the detection of glitches and their impact on MBHB parameter estimation, we removed the gaps and Galactic binaries, and for this reason, we renamed the modified data set as “Light-Spritz”.
Furthermore, the version of BBHx we are using for this paper (version 0.0.0) does not support the response model used to generate the MBHB in “Spritz”. The “Spritz” data was generated using the PhenomD waveform model Husa et al. (2016); Khan et al. (2016), with the LISA response evaluated for Keplerian orbits (time varying arm-lengths), but our version of BBHx does not include time varying arm-lengths response. Therefore, to rule out any systematic errors from this waveform mismatch Dhani et al. (2024), we replace the original “Spritz” injected signal with one with the same parameters (see Table 3 in sec. VI), but generated using the PhenomHM waveform (instead of PhenomD) with an equal arm response. We plan to include a more realistic treatment of the LISA constellation to compute the response, in a follow-up study that will analyse the full data set.

Regarding glitch modelling, Ref. Muratore (2021) showed that the mismatch between equal- and unequal/time-varying-arm models does not significantly impact glitch parameter reconstruction. For the noise, however, using equal, unequal, or time-varying arm-length models affects the high-frequency range, where the zeros of the TDI appear. To avoid related issues, we restrict our analysis to frequencies up to $2\text{\times}{10}^{-2}\text{\,}\mathrm{Hz}$ Hartwig and Muratore (2022b), which does not alter the outcome of our analysis. Additionally, the null channel $T$ cannot be well approximated using the equal-arm-length model Muratore et al. (2022); Adams and Cornish (2010). For these reasons, we perform the analysis using only the two sensitive channels, $A$ and $E$ , employing second-generation TDI under the assumption of equal arm lengths, for which sensitivity to discrepancies in arm-length modelling is negligible. For consistency, the MBHB signal is also computed using second-generation TDI variables, assuming constant and equal arm lengths.

As discussed in Sec. IV, we analyse the data in the frequency domain. However, converting the time-domain data provided in “Light-Spritz” to the frequency domain by using a fast-Fourier transform (FFT) without any pre-processing causes significant artefacts in the form of spectral leakage Hartwig (2021), linked to significant high-frequency content that leaks into the low frequency part of the “Light-Spritz” data set. We mitigate this effect by applying a low-pass filter to the data before taking the FFT. Another approach would be to apply a window function to the data before taking the FFT. This has the downside of introducing correlations between adjacent frequency bins, potentially biasing the parameter estimation using our diagonal likelihood model. Instead, by filtering the data, we obtain a diagonal likelihood (see Appendix E). Concretely, we use a first-order Butterworth filter with a cut-off frequency of $1\text{\times}{10}^{-3}\text{\,}\mathrm{Hz}$ using the scipy.signal.butter function and then use the scipy.signal.filtfilt function to apply the filter twice in opposite time directions, such that it has zero impact on the phase of the underlying data (i.e., it does not cause any time-shift of the data). Moreover, in order to not bias our analysis, we apply a frequency domain model of the filter response to the noise, glitches and GW signal templates used in the likelihood as well. Due to the double application of the filter, the Fourier transforms appearing in the likelihood are scaled by $|H(f)|^{2}$ , whereas PSDs are scaled by $|H(f)|^{4}$ , where $H(f)$ is the complex filter transfer function:

	$\displaystyle\log\mathcal{L}=$
	$\displaystyle-\frac{1}{2}\sum_{f}\Bigg{[}4\Delta f\,\frac{\left\|\tilde{d}(f)-% \tilde{h}(f,\theta)\|H(f)\|^{2}-\tilde{g}(f,\phi)\|H(f)\|^{2}\right\|^{2}}{S(f,% \lambda)\|H(f)\|^{4}}$
	$\displaystyle\qquad\qquad+\log\left\{S(f,\lambda)\|H(f)\|^{4}\right\}\Bigg{]}.$		(32)

VII.2 Noise model

To study – and model – just the noise in the “Light-Spritz” dataset, we further removed the glitch and MBHB components described in section VII.1. We then estimated the spectra using the scipy.signal.welch method with the following parameters: window=(’kaiser’, 10), nperseg=10713, noverlap=5356, detrend=’linear’. The LDC manual for “Spritz” LDC (2020) suggests that laser, backlink, acceleration and readout noises were enabled in the simulation. However, we find that the corresponding noise model provided as part of the LDC software suite (working group (2022), v.1.2.4) does not perfectly describe the “Spritz” simulated data.

To accurately fit the noise of the data set, we cross-checked against the default values for the noise sources included in the instrument simulation at the time Bayle et al. (2022) and consequently included the following noise components, expressed as PSDs:

•

Inter-spacecraft interferometer OMS noise,
$S_{\mathrm{OMS}}^{\mathrm{ISI}}=\left(N_{OMS}^{ISI}\right)^{2}$ (white),
with $N_{OMS}^{ISI}=$6.35\text{\,}\mathrm{p}\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$$
•

Reference interferometer OMS noise,
$S_{\mathrm{OMS}}^{\mathrm{RFI}}=\left(N_{OMS}^{RSI}\right)^{2}$ (white),
with $N_{OMS}^{RSI}=$3.32\text{\,}\mathrm{p}\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$$
•

Test-mass interferometer OMS noise,
$S_{\mathrm{OMS}}^{\mathrm{TMI}}=\left(N_{OMS}^{TMI}\right)^{2}$ (white),
with $N_{OMS}^{TMI}=$1.42\text{\,}\mathrm{p}\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$$
•

TM acceleration noise,
$S_{\mathrm{TM}}=\left(N_{TM}\right)^{2}\cdot\left(1+\left(\frac{0.4\times 10^{% -3}}{f}\right)^{2}\right)$ ,
with $N_{TM}=$2.4\text{\,}\mathrm{f}\mathrm{m}\mathrm{/}\mathrm{s}\mathrm{{}^{2}}% \mathrm{/}\sqrt{\mathrm{Hz}}$$
•

RFI and TMI backlink noise,
$S_{\mathrm{BL}}^{\mathrm{RFI/TMI}}=\left(N_{BL}^{\mathrm{RFI/TMI}}\right)^{2}% \cdot\left(1+\left(\frac{2\times 10^{-3}}{f}\right)^{4}\right)$ ,
with $N_{BL}^{\mathrm{RFI/TMI}}=$3\text{\,}\mathrm{p}\mathrm{m}\mathrm{/}\sqrt{% \mathrm{Hz}}$$ .

For the analysis, we define the following combinations of the noise parameters to estimate the effective noise levels:

\gamma=\begin{cases}N_{1}\equiv\left[(N_{\mathrm{OMS}}^{\mathrm{ISI}})^{2}+(N_% {\mathrm{OMS}}^{\mathrm{RFI}})^{2}+(N_{\mathrm{BL}}^{\mathrm{RFI/TMI}})^{2}% \right]^{1/2},\\ N_{2}\equiv\left[(N_{\mathrm{OMS}}^{\mathrm{TMI}})^{2}+(N_{\mathrm{BL}}^{% \mathrm{RFI/TMI}})^{2}\right]^{1/2},\\ N_{\mathrm{TM}},\end{cases}

(33)

with $N_{1}=3.32\times 10^{-12}$ , $N_{2}=7.78\times 10^{-12}$ , and $N_{\mathrm{TM}}=2.4\times 10^{-15}$ . As in section III.5, we first express the noises as equivalent displacement noise in ${\mathrm{m}}^{2}\text{\,}{\mathrm{Hz}}^{-1}$ or acceleration noise in $\mathrm{m}\mathrm{{}^{2}}\mathrm{/}\mathrm{s}^{4}\mathrm{/}\mathrm{Hz}$ . We then convert them to units of fractional frequency fluctuations as used in the LDC by applying the factors $C_{\mathrm{disp}}^{\mathrm{ffd}}=\left({2\pi f}/{c}\right)^{2}$ and $C_{\mathrm{acc}}^{\mathrm{ffd}}=\left({c2\pi f}\right)^{-2}$ , where $c$ is again the speed of light.

Fig. 16 shows the ratio between the estimated spectra (the PSD computed from the “Spritz” only-noise data), the model included with the LDC tools, and the model we used. We see artefacts around the zeros of the TDI transfer function for both plots, owing to the use of an equal arm approximation in the model we are using to fit the noise. For our model, we estimated the armlength to be used in the equal arm model as an average over all six time-varying armlengths, evaluating to $L/c\approx 8.323s$ . Whereas we used the default value set in the LDC tools for the LDC model. This causes slightly increased artefacts for the LDC model around the zeros. We also observe that the LDC model slightly underestimates the noise in the upper half of the frequency band.

To conclude, some of these noises have different TDI transfer functions, which are derived in Quang Nam et al. (2023b) for TDI $X,Y,Z$ as well as $A,E,T$ . For our work, we only need the transfer functions for $A$ and $E$ , which turn out to be identical under our assumptions. We summarize them here for convenience. They include a common TDI factor for second generation TDI:

C_{xx}=16\sin^{2}(2\pi fL)\sin^{2}(4\pi fL).

(34)

Using this, we have

$\displaystyle\|H_{\text{ISI/RFI (OMS)}}^{AA}\|^{2}$	$\displaystyle=2C_{xx}\left(2+\cos(2\pi fL)\right)$	(35)
$\displaystyle\|H_{\text{TMI (OMS)}}^{AA}$	$\displaystyle\|^{2}=$
	$\displaystyle C_{xx}\left(3+2\cos(2\pi fL)+\cos(4\pi fL)\right)$	(36)
$\displaystyle\|H_{\text{TM}}^{AA}\|^{2}$	$\displaystyle=4\cdot\|H_{\text{TMI (OMS)}}^{AA}\|^{2},$	(37)

such that the total noise can be described by

\begin{split}S_{AA}^{\mathrm{total}}&=S_{EE}^{\mathrm{total}}=C_{\mathrm{disp}% }^{\mathrm{ffd}}\Big{[}|H_{\mathrm{TMI(OMS)}}^{AA}|^{2}\left(S_{\mathrm{OMS}}^% {\mathrm{TMI}}+S_{\mathrm{BL}}^{\mathrm{TMI}}\right)\\ &+|H_{\mathrm{ISI/RFI(OMS)}}^{AA}|^{2}\left(S_{\mathrm{OMS}}^{\mathrm{ISI}}+S_% {\mathrm{OMS}}^{\mathrm{RFI}}+S_{\mathrm{BL}}^{\mathrm{RFI}}\right)\Big{]}\\ &+C_{\mathrm{acc}}^{\mathrm{ffd}}|H_{\mathrm{TM}}^{AA}|^{2}S_{\mathrm{TM}}.% \end{split}

(38)

VIII Analysis of the data set

The search-phase algorithm is visualized in Fig. 17 and explained in the subsection VIII.1, while the parameter estimation one is visualized in Fig. 18 and explained in section VIII.2.

VIII.1 Search phase

We divide the time series into segments of approximately one day, with a $20\%$ overlap between adjacent segments. We therefore segment the “Light-Spritz” dataset into 29 stretches of data⁷⁷7We start the counting of the stretches of data from number 0. and we regenerate a time grid spanning from $t_{0}=$0\text{\,}\mathrm{s}$$ to $t_{end}=$2\,675\,890\text{\,}\mathrm{s}$$ . We consider broad prior ranges for the glitch parameters, allowing the distribution of $A_{0}$ and $\beta_{0}$ to extend beyond what is observed in LPF (Fig. 1), aiming to remain as agnostic as possible. The priors and the corresponding true values on the glitch and noise parameters for the search phase are reported in Table 4. The true value of the amplitude and decay times for the glitch parameters are taken from the “Spritz” LDC manual LISA LDC Team (nd), whereas the injection times are not explicitly given. In particular the amplitude $\Delta v$ is $2.20\text{\,}\mathrm{pm}\text{\,}{\mathrm{s}}^{-1}$ , while the two decay times are $\tau_{1}=$10\text{\,}\mathrm{s}$$ and $\tau_{2}=$11\text{\,}\mathrm{s}$$ . Our search is split in two phases.

Parameter	True Value	Lower Bound	Upper Bound
$\tau_{0}$ [s]	—	t-start	t-end
$\ln(A_{0}/[m\,s^{-2}])$	-27.58	-35	-20
$\beta_{0}$ [s]	10.5	1	$1\text{\times}{10}^{5}$
$N_{1}$ [ $\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$ ]	$7.77\text{\times}{10}^{-12}$	$7\text{\times}{10}^{-12}$	$8\text{\times}{10}^{-12}$
$N_{2}$ [ $\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$ ]	$3.32\text{\times}{10}^{-12}$	$2.5\text{\times}{10}^{-12}$	$3.5\text{\times}{10}^{-12}$
$N_{TM}$ [ $\mathrm{m}\mathrm{/}\mathrm{s}\mathrm{{}^{2}}\mathrm{/}\sqrt{\mathrm{Hz}}$ ]	$2.4\text{\times}{10}^{-15}$	$2\text{\times}{10}^{-15}$	$3\text{\times}{10}^{-15}$

Table 4: Table of priors for the glitch and noise parameters along with the true values. The variables ‘t-start’ and ‘t-end’ change depending on the stretch considered. Note that the following relationships hold between a double decay exponential glitch and the shapelet model used to fit it:

A_{0}=\Delta v\frac{1}{1+\tau_{2}/\tau_{1}}

and

\beta_{0}=\frac{\tau_{1}+\tau_{2}}{2}

VIII.1.1 First phase

First, we fit only for both noise and glitches in each segment separately. We use the in-model and out-of-model moves described in sec. IV.1, to explore the parameter space. The number of glitches in each segment is estimated between 0 and 1 (i.e. the number of leaves in the glitch branch is 0 or 1). The search includes a stopping criterion such that if after 1000 iterations, more than $70\%$ of the samples have no glitch leaves, we conclude that no glitches are present in the segment and proceed to analyze the next one. Otherwise the segment is flagged as containing a possible glitch signal.

Table 5 reports the results of the first phase of the search pipeline, showing in which stretches of data glitches have been found and their SNRs. Note that the SNR reported in the table is estimated through the simultaneous fitting of glitches and noise. Specifically, to calculate the SNR, we compute the mean of the posterior distribution of SNR values for each data stretch during the search phase. This mean is obtained by drawing 3000 random samples from the posterior distribution of the glitch and noise parameters. However, the instrumental noise is not sufficiently constrained; therefore, these SNRs have to be considered as a rough estimate. Indeed, we observed that only a day of data is not enough to properly estimate and constrain all the noise parameters. This is particularly true in stretches of data where the late inspiral part of the MBHB is present (N26-N28). Therefore, we need to properly constrain the noise when processing the full month of data. This is further discussed in the Appendix D.

To perform this initial search, we used 25 walkers and 10 temperatures, and the run for each segment takes about 20 minutes.

N	Time [s]	% of 0	% of 1	Glitch
		Leaves	Leaves	SNR
0	0.0 - 107995	86	13	0.60
1	72000 - 197995	66	36	1.27
2	162000 - 287995	87	13	0.39
3	252000 - 377995	86	14	0.47
4	342000 - 467995	85	15	0.50
5	432000 - 557995	83	16	0.54
\rowcolorred!30 6	522000 - 647995	10	90	294
7	612000 - 737995	58	42	1.28
8	702000 - 827995	87	13	0.42
9	792000 - 917995	87	13	0.43
10	882000 - 1007995	56	44	1.42
11	972000 - 1097995	79	21	1.46
12	1062000 - 1187995	85	15	0.51
13	1152000 - 1277995	84	15	0.52
14	1242000 - 1367995	86	14	0.45
15	1332000 - 1457995	87	12	0.40
16	1422000 - 1547995	87	13	0.45
17	1512000 - 1637995	80	20	0.65
18	1602000 - 1727995	84	16	0.65
19	1692000 - 1817995	86	14	0.48
20	1782000 - 1907995	86	14	0.47
21	1872000 - 1997995	85	15	3.34
22	1962000 - 2087995	85	15	0.44
\rowcolorred!30 23	2052000 - 2177995	3	97	269.14
\rowcolorred!30 24	2142000 - 2267995	6	93	265.68
25	2232000 - 2357995	81	18	0.77
\rowcolorcyan!3026	2322000 - 2447995	8	92	3.74
\rowcolorcyan!3027	2412000 - 2537995	3	97	5.85
\rowcolorred!30 28	2502000 - 2627995	1	98	188.99
\rowcoloryellow!50 29	2592000 - 2675890	1	99	520.26

Table 5: Table showing the number of leaves in each model and the glitch SNR for each stretch of data analyzed.

N

labels the segment considered. The “% of 0 or 1 leaf” columns indicate the fraction of samples with 0 or 1 leaf, respectively. The orange lines indicate stretches of data where glitches were identified, the cyan lines mark regions where the inspiral phase of the MBHB becomes significant, and the yellow line highlights the stretch where the merger of the MBHB is present.

We observe five stretches with high SNR signals. Of these, four have been identified as glitches (orange line) and one as a MBHB (yellow line). All these stretches have at least $90\%$ of samples with only one leaf. Additionally, two other stretches (cyan line) have more than $92\%$ of samples with one leaf, although their estimated signal SNRs are much lower. This occurs because the inspiral phase of the MBHB waveform starts to contribute a non-negligible SNR, triggering the algorithm with evidence of a low-SNR signal. In the remaining stretches, the SNR values result from random noise fluctuation.

In the first phase of the search, the glitch model serves as a trigger to detect the presence of a signal in each stretch, whether it corresponds to a glitch or an MBHB signal. This approach is effective because the MBHB in “Light-Spritz” has a high SNR, leading to behavior similar to matched filtering, where the overlap between our template bank and the GW signal increases with the SNR (the characteristic frequency of the glitch and the instantaneous GW frequency of the MBHB). However, in cases with lower SNR MBHBs, the search may not be sensitive enough to detect such signals. To resolve ambiguities between the presence of a glitch or an MBHB signal in each stretch, the second phase of the search evaluates the Bayes factor comparing the glitch and MBHB hypotheses.

VIII.1.2 Second phase

The first phase of our search aims to identify the presence of a signal in each stretch of data. In the second phase, we determine whether the signal is better described as a glitch or a MBHB. To this end, we perform a dedicated MBHB search to assess the presence of a GW signal in stretches where a signal was detected or where the estimated SNR is sufficiently high (N6, N23, N24, N26, N27, N28, and N29). We then compute the Bayes factor comparing the MBHB and glitch hypotheses⁸⁸8Note that this analysis could alternatively be performed using the RJ-MCMC framework, which would allow for automatic computation of Bayes factors between glitches and MBHBs by leveraging the Erebor Global Fit infrastructure Katz et al. (2024). However, incorporating this additional complexity is not necessary for the current analysis.. For this, we employ both the stepping-stone method Maturana-Russel et al. (2019) and thermodynamic integration Karnesis et al. (2023). This information is subsequently used to initialize the walkers for the parameter estimation phase. We then compare the Bayes factors between the glitch and MBHB hypotheses for each of these stretches. To this end, we use 30 temperatures and run the sampler with a fixed-dimensional parameter space. We fit only for the signal, assuming the noise to be at a fiducial level, following the strategy outlined in Katz et al. (2024). This approach is robust given the high SNR of both glitches and MBHBs, and no ambiguous Bayes factor values were observed that would suggest the simultaneous presence of both components. In potential edge cases, the existing methodology in the Erebor Global Fit framework Katz et al. (2024) offers a viable and complementary strategy for MBHB detection.

We show in Table 6 three representative examples of the log Bayes factor for stretches N6, N24, and N29, where a signal (either a glitch or an MBHB) was identified. For brevity, we report only the results obtained using the stepping-stone method. In the Appendix F are reported the results obtained using the thermodynamic integration for completeness. The Bayes factors across all considered stretches were found to be highly decisive, strongly favouring one hypothesis over the other, either the presence of a glitch or that of an MBHB. For this reason, we present only a subset of representative results here, as the conclusions were clear-cut in all cases. Specifically, for stretches N6 and N24, the presence of a glitch is strongly favored, while for stretch N29, the MBHB hypothesis is preferred.

	N6	N24	N29
$\log\mathcal{B}_{\mathrm{MBHB}/\mathrm{glitch}}$	$\sim-6.1\times 10^{4}$	$\sim-3.9\times 10^{4}$	$\sim 1.5\times 10^{7}$
$d\log\mathcal{B}_{\mathrm{MBHB}/\mathrm{glitch}}$	$\sim 0.61$	$\sim 0.22$	$\sim 0.44$

Table 6: Log Bayes factors comparing the MBHB and glitch models for datasets N6, N24, and N29 using the stepping-stone method.

d\log\mathcal{B}_{\mathrm{MBHB}/\mathrm{glitch}}

gives an estimate on the error on

\log\mathcal{B}_{\mathrm{MBHB}/\mathrm{glitch}}

During the search phase, along with computing the Bayes factor we also obtain a rough estimate of the parameters of the MBHB. Doing parameter estimation for the MBHB using only a day of data highlights several key considerations regarding frequency domain analysis and the use of the discrete Fourier transform. One primary concern is spectral leakage, particularly evident when the pre-merger/inspiral phase of the signal is analyzed, which necessitates the use of windowing techniques to mitigate its effects Davies et al. (2024). Indeed, when analyzing frequency-domain waveforms over short time segments, particular caution is required to ensure accuracy. However, in the context of this research, we did not window the data as it was not strictly necessary since we were not interested in characterizing the pre-merger part of the signals⁹⁹9As illustrated in Appendix E, filtering the data favours no correlation between adjacent frequency bins.. These issues reflect a fundamental challenge in performing low-latency analyses of MBHB signals Cornish (2022), but once again we do not need to be so careful in the search phase, as the final characterization of the signals uses the full data set.

To conclude on a more technical side, to be able to perform the analysis of the MBHB with BBHx but restricting it to a single day, we need to account for the time window considered when generating the MBHB template with BBHx. This is achieved by shifting the phase with respect to the starting time of our window, and by defining a start and an end time of the observation to cut the waveform, such that any frequency components of the signal that are emitted outside of this range are excluded.

VIII.2 Parameter estimation phase

For the final parameter estimation phase we consider the full month of data. We include an MBHB template in our likelihood and consider the same prior as the one reported in Table 3, but restrict the prior on the time of coalescence between $t_{ref}-24s<t_{ref}<t_{ref}+24s$ according to the MBHB posteriors found during the search phase. We also consider a slightly restricted prior range for the glitch parameters based on the findings from the search phase as shown in Table 7. In particular, we restrict the times of injection and exclude combinations of $\beta_{0}$ and $A_{0}$ that do not reflect the glitch distribution found in LPF. If such glitches were present, they would have already been identified during the search phase given the high SNR of those glitches.

We then set the minimum number of glitches for the reversible jump to match the number of high SNR glitches identified in the search phase, and we set the maximum number of glitches to be double the minimum to allow for glitches that could have been missed during the search.

We consider the same in-model and out-of-model moves used during the search phase for the glitches. The in-model moves get updated every 20 iterations up to a total of 50000 iterations. Instead, for the in-model Gaussian move proposal of the MBHB we estimate the covariance matrix using samples from the MBHB parameter posteriors of the search phase. We also incorporate the covariance matrices computed from the three estimated glitches -during the search phase- to define an additional reversible jump move proposal (see section IV.1.2) which proposes the birth or death of a glitch based on those identified during the search. To be on the safe side, we still allow $40\%$ of proposed glitch parameters to be drawn from the full prior range, ensuring that we do not miss any potential glitches.

Finally, we initialise the MCMC walkers using samples found during the search phase, and we estimate the parameters of the noise, signal, and glitches simultaneously following a Gibbs sampling approach Karnesis et al. (2023).

The posteriors for the three glitches are shown in Figs. 19,20,21, for the MBHB are shown in Figs. 22 and for the noise parameters in Fig. 23. The MBHB found has an SNR of 1522 while the glitches have SNRs of around 300. For all the posteriors, we only plot the last 500 samples of each walker, as we consider the previous samples, where we perform recursive updates of the in-model moves, as a burn-in period.

Parameter	Lower	Upper
	Bound	Bound
$\tau_{0}$ [s]	432000	2627995
$ln(A_{0}/[ms^{-2}])$	-32	-20
$\beta_{0}$ [s]	1	1e4
$N_{1}$ [ $\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$ ]	$7\text{\times}{10}^{-12}$	$8\text{\times}{10}^{-12}$
$N_{2}$ [ $\mathrm{m}\mathrm{/}\sqrt{\mathrm{Hz}}$ ]	$2.5\text{\times}{10}^{-12}$	$3.5\text{\times}{10}^{-12}$
$N_{TM}$ [ $\mathrm{m}\mathrm{/}\mathrm{s}\mathrm{{}^{2}}\mathrm{/}\sqrt{\mathrm{Hz}}$ ]	$2.2\text{\times}{10}^{-15}$	$2.7\text{\times}{10}^{-15}$

Table 7: Table of priors for the glitch and noise parameters. The true values are given in Table 4.

IX Conclusion and future perspective

In this work, we have demonstrated a potential Global Fit approach for LISA to estimate MBHB parameters in the presence of noise artefacts by implementing a search pipeline to find and fit for glitches along with the noise properties of the instrument and the parameters of an MBHB, using realistic simulated LISA data. Our methodology provides a general framework for detecting and modelling glitches, ensuring that such instrumental artefacts do not compromise LISA’s scientific output.

We have developed and validated our pipeline in two steps. First, we generated synthetic data where multiple glitches of varying SNRs, including those closely spaced in time, were injected to test the effectiveness of our reversible jump approach in handling diverse glitch morphologies. Our injected glitches were randomly drawn from the LPF glitch distribution to provide a realistic test scenario. We established an approximate threshold for distinguishing glitches from instrumental noise and investigated the impact of low, medium and high SNR glitches on the parameter estimation of a fiducial MBHB with SNR of 1152. We found that unfitted glitches with SNR $\gtrsim$ 60 introduce non-negligible biases in MBHB parameter recovery, especially for low-SNR MBHB signals. For instance, for an unfitted glitch with SNR 21 we observed a slight shift in the parameter posterior relative to the injected values, but this shift was small compared to the statistical uncertainties, whereas for a glitch with SNR 72, the true values were excluded at high confidence from the posterior support for several parameters, such as the total mass or the sky location. For the extreme case of an SNR 1544 glitch, biases were very large and the posterior was completely disjoined to the case were no glitch is present. Moreover, we noticed that the relative SNR between the glitch and the MBHB signal, as well as the glitch’s temporal proximity to the merger, are critical factors in determining the extent of bias. However, these are not the sole contributors. For instance, when comparing the posterior distributions of MBHB systems with SNRs 1152 and 320, the glitch number 41 in the LPF catalogue appears to cause a more pronounced bias in the recovery of parameters for the low-SNR system hinting that also the relationship between the glitch frequency and the instananeous frequency of the MBHB plays a role. When fitting for the glitch(-es) (single glitch as well as multiple glitches), we corrected for these biases and estimated the parameters of the MBHB correctly, with almost the same accuracy when no glitch was present. This study stresses the importance of accounting for glitches to fully exploit the tremendous scientific potential of MBHB observations with LISA—whether for fundamental physics Berti et al. (2006, 2016); Chamberlain and Yunes (2017); Barausse et al. (2020b); Maggio et al. (2020); Bhagwat et al. (2022); Corman et al. (2022); Toubiana et al. (2024b); Pitte et al. (2024), astrophysics Gair et al. (2011); Sesana et al. (2011); Klein et al. (2016); Dayal et al. (2019); Barausse et al. (2020a); Toubiana et al. (2021); Chen et al. (2022); Fang and Yang (2023); Spadaro et al. (2025); Langen et al. (2024); Toubiana et al. (2024a), or cosmology Caprini and Tamanini (2016); Tamanini et al. (2016); Mangiagli et al. (2023).

Second, we applied our search algorithm to a simplified version of the “Spritz” dataset, which was obtained by removing the Galactic binaries and gaps from the original dataset. Our results show that we can successfully implement a search and a parameter estimation pipeline to detect and model glitches present in the data while recovering both the instrumental noise and the MBHB parameters without bias. Additionally, we provided a detailed noise model for this dataset after identifying inconsistencies between the official LDC “Spritz” documentation LDC (2022) and the actual procedure used to generate data for this mock data set. This information will be valuable for future researchers analyzing the same dataset.

A key aspect that warrants further investigation is the accuracy and efficiency of the GMM clustering approach used to build proposal distributions for glitches in the parameter estimation phase, particularly in determining the number of components. This is particularly important in view of including the search and parameter estimation for glitches, developed in this work, in the framework of the Erebor Global Fit Katz et al. (2024), in which multiple MBHBs and other types of GW source are present. In the current implementation, the number of clusters is fixed to match the maximum number of leaves (components) in the model, without attempting to infer the true number of underlying glitch populations. While this approach offers flexibility, and has demonstrated to work very well in our framework, it may introduce inefficiencies when the number of GMM components significantly exceeds the actual number of distinct glitches. More specifically, in scenarios where the GMM overestimates the number of glitch clusters, the model may fragment individual glitch posteriors into sub-posteriors, which can hinder efficient sampling. As a result, MCMC chains may struggle to traverse the full posterior landscape, especially across different glitch modes, even if local acceptance rates remain high. We have begun exploring alternative strategies to dynamically estimate the appropriate number of clusters, potentially improving the exchange of samples across posterior modes and enhancing the overall sampling efficiency. This might play a more significant role when more GW sources or more complicated noise components, with respect to the “Light-Spritz” data set, are considered.

Nonetheless, the glitch search, determination, and parameter estimation components developed in this work will form the foundation of the glitch handling module within Erebor. From Erebor, we will incorporate the already developed modules for MBHB searches and parameter estimation as well as the search and parameter estimation for Galactic binaries, noise and foreground fitting. Over time, we expect to refine both the glitch and MBHB modules based on the insights gained here, particularly in trans-dimensional sampling strategies. We will also update the noise fitting by leveraging the more complete noise models reported in this work.

To conclude, future work will also extend this approach to account for instrumental non-stationarities, such as data gaps and the presence of the Galactic foreground so that data with the more realistic properties of the “Spritz” data set can be handled. While our study serves as a guideline for glitch mitigation in LISA, it does not encompass all possible MBHB populations, and our reference glitch distribution is based on LPF data. However, the shapelet-based approach should readily adapt to different glitch morphologies, with respect to the signal exponential ones addressed in this paper, present in the data in particular by leveraging higher order shapelets, but its efficacy needs to be properly quantified. However, since LISA will use similar technology to LPF, at least for the GRS, the acceleration glitches are unlikely to differ substantially, so the results obtained here should represent the expected performance on actual LISA data.

All code used in this work is open-source and can be accessed at https://github.com/martinaAEI/artifacts.

Acknowledgements.

M.M. gratefully acknowledges the support of the German Space Agency, DLR. The work is supported by the Federal Ministry for Economic Affairs and Climate Action based on a resolution of the German Bundestag (Project Ref. No. FKZ 50 OQ 2301). The authors thank the following people (alphabetically ordered) for useful discussion in improving the work: Quentin Baghi, Jean-Baptiste Bayle, Ollie Burke, Eleonora Castelli, Raffi Enficiaud, Hector Hestelles, Niklas Houba, Nikos Karnesis, Lorenzo Pompili, Lorenzo Sala, Elise Sangër, Alessandro Santini, Lorenzo Speri, Alice Spadaro, Mauro Pieroni and Sebastian Völkel.

Appendix A Glitch model

A.1 Derivation of shapelet model in time and frequency domains

Shapelets can be defined in the time domain as follows:

\psi_{n}(t)=2\frac{t}{n}e^{-\frac{t}{n}}L^{1}_{n-1}\bigg{(}2\frac{t}{n}\bigg{)% }\Theta(t),

(39)

with $L_{n}^{\alpha}$ denoting a generalized Laguerre polynomial. A glitch perturbation is then modeled by a linear combination of shapelets:

g(t)=\sum_{i=0}^{P}A_{i}\psi_{n_{i}}\bigg{(}\frac{t-\tau_{i}}{\beta_{i}}\bigg{% )},

(40)

with $P+1$ being the number of shapelets. Substituting Eq. (39) into Eq. (40) we have:

	$\displaystyle g(t)$	$\displaystyle=\sum_{i=0}^{P}2A_{i}\frac{1}{n_{i}}\frac{t-\tau_{i}}{\beta_{i}}e% ^{-\frac{1}{n_{i}}\frac{t-\tau_{i}}{\beta_{i}}}$
		$\displaystyle\quad\times L^{1}_{n_{i}-1}\bigg{(}2\frac{1}{n_{i}}\frac{t-\tau_{% i}}{\beta_{i}}\bigg{)}\Theta\bigg{(}\frac{t-\tau_{i}}{\beta_{i}}\bigg{)},$		(41)

note that Eq. (A.1) refers to a glitch in acceleration but we need to derive the correspondent fractional frequency shift, which means integrating $g(t)$ and dividing by the velocity of light $c$ .

To perform the integral, we can consider the generating function:

\sum_{n=0}^{\infty}v^{n}L_{n}^{\alpha}(x)=\frac{1}{(1-v)^{\alpha+1}}e^{-\frac{% xv}{1-v}}.

(42)

To align it with Eq. (39), we adjust the index of the series as follows:

\sum_{n=1}^{\infty}v^{n}L_{n-1}^{\alpha}(x)=\frac{v}{(1-v)^{\alpha+1}}e^{\frac% {-xv}{1-v}}.

(43)

and then define

	$\displaystyle G(t,\tau,n,\beta,A,v)$
	$\displaystyle\hskip 14.22636pt=\frac{2A}{c}\int_{0}^{t}\left(\frac{t^{\prime}-% \tau}{n\beta}\right)\exp\left[{-\frac{t^{\prime}-\tau}{n\beta}}\right]$
	$\displaystyle\hskip 28.45274pt\times\frac{v}{(1-v)^{\alpha+1}}\exp\left[{-2% \frac{t^{\prime}-\tau}{n\beta}\frac{v}{(1-v)}}\right]$
	$\displaystyle\hskip 56.9055pt\times\Theta\left(\frac{t^{\prime}-\tau}{\beta}% \right)\,{\rm d}t^{\prime}.$		(44)

This expression enables us to obtain different shapelets orders by rewriting the relationship between the left and right components of Eq. (44) as follows:

	$\displaystyle\frac{\Delta\nu_{g}}{\nu_{0}}(t,\vec{\tau},\vec{n},\vec{\beta},% \vec{A})\equiv\int_{0}^{t}\frac{g(t^{\prime})}{c}\,dt^{\prime}$
	$\displaystyle\hskip 28.45274pt=\sum_{i=0}^{P}\frac{1}{n_{i}!}\left[\frac{{\rm d% }^{n_{i}}}{{\rm d}v^{n_{i}}}G(t,\tau_{i},n_{i},\beta_{i},A_{i},v)\right]_{\|v=0}.$		(45)

Eq. (44) can be explicitely written as:

	$\displaystyle G(t,\tau,n,\beta,A,v)=\frac{2vA}{c(1-v)^{\alpha+1}}n\beta$
	$\displaystyle\hskip 56.9055pt\times\int_{-\frac{\tau}{n\beta}}^{\frac{t-\tau}{% n\beta}}ue^{-u}e^{\frac{-2uv}{1-v}}\Theta(nu)\,{\rm d}u$
	$\displaystyle\hskip 28.45274pt=-\frac{2vnA\beta(1-v)^{-\alpha}}{c(1+v)^{2}}% \times\Theta\bigg{(}\frac{t-\tau}{\beta}\bigg{)}$
	$\displaystyle\times\bigg{(}-1+v+\frac{e^{\frac{(1+v)(t-\tau)}{n(-1+v)\beta}}(t% +tv+n\beta-nv\beta-(1+v)\tau)}{n\beta}\bigg{)}.$		(46)

In this way, the generating function is used to perform the integral. Equation 46 needs to be computed for each value of $A_{i},\tau_{i}$ and $\beta_{i}$ that is needed such that the sum over the shapelets can be performed.

We can now derive the Fourier transform of $\frac{\Delta v_{g}}{\nu_{0}}(t,\vec{\tau},\vec{n},\vec{\beta},\vec{A})$ . This reads:

	$\displaystyle\frac{\Delta\tilde{\nu}_{g}}{\nu_{0}}(f,\vec{\tau},\vec{n},\vec{% \beta},\vec{A})=$		(47)
	$\displaystyle\hskip 28.45274pt\sum_{i=0}^{P}\frac{1}{n_{i}!}\left[\frac{{\rm d% }^{n_{i}}}{{\rm d}v^{n_{i}}}\tilde{G}(f,\tau_{i},n_{i},\beta_{i},A_{i},v)% \right]_{\|v=0}.$		(48)

Considering Eq. 46, we obtain:

	$\displaystyle\tilde{G}(f,\tau,n,$	$\displaystyle\beta,A,v)=$
		$\displaystyle\frac{iAn(1-v)^{1-\alpha}v\beta e^{-2if\pi\tau}}{cf\pi(i(1+v)+2fn% \pi(-1+v)\beta)^{2}},$		(49)

and therefore,

	$\displaystyle\frac{\Delta\tilde{\nu}_{g}}{\nu_{0}}(f,\vec{\tau},\vec{n},\vec{% \beta},\vec{A})=\noindent$		(50)
	$\displaystyle\sum_{i=0}^{P}\frac{1}{n_{i}!}\left[\frac{{\rm d}^{n_{i}}}{{\rm d% }v^{n_{i}}}\bigg{(}\frac{iAn(1-v)^{1-\alpha}v\beta e^{-2if\pi\tau}}{cf\pi(i(1+% v)+2fn\pi(-1+v)\beta)^{2}}\bigg{)}\right]_{\|v=0}.$		(51)

The zero order shapelet ( $i=0$ ) in the frequency domain is

\frac{\Delta\tilde{\nu}_{g}}{\nu_{0}}(f,\tau_{0},\beta_{0},A_{0})=\frac{iA_{0}% \beta_{0}e^{2if\pi\tau_{0}}}{cf\pi(i-2f\pi\beta_{0})^{2}}.

(52)

A.2 Glitch model in time and frequency domain and respective TDI transfer functions

The effect of a typical LPF glitch is evaluated starting from a model for its differential acceleration, $\Delta g(t)$ , Sala (2023); Muratore (2021):

h_{g}(t)=\frac{\Delta v}{\tau^{2}}t{{}^{\prime}}e^{-\frac{t{{}^{\prime}}}{\tau% }}\Theta(t^{\prime}),\quad\quad t{{}^{\prime}}=t-t_{0}

(53)

where $t_{0}$ is the occurrence time, $\Theta(t)$ is the Heaviside step function, $\Delta v$ is the total transferred impulse per unit mass, and the $\tau$ parameter determines the duration of the event. We can also express the fractional frequency shift from the glitch as:

\frac{\Delta\nu_{g}(t)}{\nu_{0}}=\frac{\Delta v}{c}\Bigg{(}1-\frac{e^{\frac{-t% +T}{\tau}}(t-T+\tau)}{\tau}\Bigg{)}\Theta(t-t_{0}).

(54)

From a frequency-domain perspective, $\Delta v$ is also the zero-frequency limit of the Fourier transform of the $h_{g}(t)$ template. This is relevant, since it implies that it has a strong low-frequency component. The Fourier transform of the template is:

\frac{\Delta\nu_{g}(f)}{\nu_{0}}=-\frac{\Delta ve^{-i2\pi ft_{0}}}{ci2\pi f% \left(-i+2\pi f\tau\right)^{2}},

(55)

where we need to divide by $c\omega$ to express it as fractional frequency shift. The expression is similar to the shapelet when $n_{0}=1$ (Eq. 10) with $\Delta v=2A_{0}\beta_{0}$ and $t_{0}=\tau_{0}$ and $\tau=\beta_{0}$ .

Considering now Eq. III.1 in section III.3, we can model the TDI response to a glitch as a product between the glitch expression in the frequency domain (Eq. 55) and the transfer function $\mathcal{T}_{X,Y,Z}^{g}(f)$ which takes into account TDI (see section III.3). We provide below the expressions for the first-generation TDI transfer functions, noting $\mathcal{F}{D_{ij}}=e^{-i(2\pi f)L}$ :

In case of a glitch happening on $TM_{21}$ we have:


$\displaystyle\mathcal{T}_{X_{1}}^{g}(f)$	$\displaystyle=-2e^{-3iL(2\pi f)}(-1+e^{2iL(2\pi f)})$	(56a)
$\displaystyle\mathcal{T}_{Y_{1}}^{g}(f)$	$\displaystyle=1-e^{-4iL(2\pi f)}$	(56b)
$\displaystyle\mathcal{T}_{Z_{1}}^{g}(f)$	$\displaystyle=0$	(56c)

In case the same glitch is happening on $TM_{13}$


$\displaystyle\mathcal{T}_{X_{1}}^{g}(f)$	$\displaystyle=1-e^{-4iL(2\pi f)}$	(57a)
$\displaystyle\mathcal{T}_{Y_{1}}^{g}(f)$	$\displaystyle=0$	(57b)
$\displaystyle\mathcal{T}_{Z_{1}}^{g}(f)$	$\displaystyle=-2e^{-3iL(2\pi f)}(-1+e^{2iL(2\pi f)})$	(57c)

or on $TM_{31}$


$\displaystyle\mathcal{T}_{X_{1}}^{g}(f)$	$\displaystyle=4ie^{-2iL(2\pi f)}\sin(L(2\pi f))$	(58a)
$\displaystyle\mathcal{T}_{Y_{1}}^{g}(f)$	$\displaystyle=0$	(58b)
$\displaystyle\mathcal{T}_{Z_{1}}^{g}(f)$	$\displaystyle=-1+e^{-4iL(2\pi f)}$	(58c)

or on $TM_{23}$


$\displaystyle\mathcal{T}_{X_{1}}^{g}(f)$	$\displaystyle=0$	(59a)
$\displaystyle\mathcal{T}_{Y_{1}}^{g}(f)$	$\displaystyle=-1+e^{-4iL(2\pi f)}$	(59b)
$\displaystyle\mathcal{T}_{Z_{1}}^{g}(f)$	$\displaystyle=4ie^{-2iL(2\pi f)}\sin(L(2\pi f))$	(59c)

and finally in case of $TM_{23}$ :


$\displaystyle\mathcal{T}_{X_{1}}^{g}(f)$	$\displaystyle=0$	(60a)
$\displaystyle\mathcal{T}_{Y_{1}}^{g}(f)$	$\displaystyle=-2e^{-3iL(2\pi f)}(-1+e^{2iL(2\pi f)})$	(60b)
$\displaystyle\mathcal{T}_{Z_{1}}^{g}(f)$	$\displaystyle=1-e^{-4iL(2\pi f)}$	(60c)

For the second generation TDI a similar computation can be done using Eq. 3 in section III.3.

Appendix B Posteriors of SNR 1152 MBHB with high SNR glitch of 1544 at different injection times with respect to its merger time

The impact of a glitch in the recovery of a MBHB is maximised when the glitch occurs closer to the merger; this, in practice, means that the same type of glitch would not have as great an impact if it happened at a different time. To show this, we also report below the impact on parameter estimation of the fiducial MBHB with SNR 1152 and redshift 4 – considered also in the analysis of the “Light-Spritz” data set – with the 1544 SNR glitch injected at different times with respect to maximum overlap computed in section VI: five hours before the merger and five minutes after the merger (Fig. 24). In all cases, noisy data have been used, which introduce some random statistical fluctuations that are visible in the posteriors.

Looking at Fig. 24, where the glitch has been injected five hours before the merger time, we can see that not fitting for the glitch causes less bias with respect to the glitch occurring five minutes after the merger. We expect that glitches happening earlier (or later with respect to the merger time) would have an even smaller effect. If the glitch occurs closer to or at the merger, not fitting for it clearly biases the parameters of the MBHB, such as the masses, spins or sky positions. This behaviour stems from the fact that the merger is the loudest part of the signal, and drives the parameter estimation of the MBHB, so glitches that corrupt that part of the signal have a stronger effect.

Appendix C Posteriors of SNR 320 MBHB with glitches $N41$ and $N3$ injected at their maximum overlap

In Fig. 10 are reported the posteriors of an SNR 320 MBHB when the two glitches corresponding to glitch $\#41$ and $\#3$ of the LPF catalogue are injected. The parameters of the MBHB are reported in Table 8. The glitches were injected at their maximum overlap, that is $0.19[h]$ after the merger for glitch $\#41$ , and $0.11[h]$ after the merger for glitch $\#3$ , as can be seen from the Figure 25. It is possible to observe that sky parameters have multimodality. Indeed, there are up to eight modes in the sky that are exactly degenerate for a fixed LISA configuration and non-evolving sources. However, this degeneracy is broken by LISA’s motion and changes in frequency. High-mass and short-duration sources, like this case, are more likely to exhibit such degeneracies. The degeneracy primarily involves sky location and polarization ( $\beta_{d}$ , $\psi$ , $\lambda_{d}$ and $\cos(\iota)$ ) Marsat et al. (2021); Katz (2022), as can also be seen in the posteriors plot showing the extrinsic parameters (right plot Fig. 26). The left plot in Fig. 26 shows the posteriors for the intrinsic parameters.

Parameter	True Value
$ln(M_{t}/M_{\odot})$	17.91
$q$	0.33
$a_{1}$	0.3
$a_{2}$	0.4
$d_{L}[Gpc]$	47.6
$\varphi_{ref}$ $[rad]$	1.2
$\cos\iota$	0.29
$\lambda$ $[rad]$	2.0
$\sin\beta$	0.82
$\psi$ $[rad]$	1.6
$t_{ref}[h]$	30

Table 8: Massive Black Hole Binary true values. Values of the masses are given in the detector frame.

Appendix D Discussion on posteriors for noise estimated during the search phase

Fig. 27 shows the posterior distributions for the three noise parameters ( $N_{1}$ , $N_{2}$ , and $N_{\mathrm{TM}}$ ), estimated during the search for segments N0, N2, N5, N8, N10, N13, N15, N20, and N26, along with the true values reported in section VII.2. The TM acceleration noise is not well characterized, showing a bias away from the true value. This bias increases for data chunks that include the late inspiral phase of the MBHB (segments N26–N28), as this occurs in the frequency range around $10^{-3}$ Hz where TM noise dominates, as also visible in Fig. 28.
The $N_{2}$ component spans the entire prior range, indicating poor constraint. In contrast, $N_{1}$ is better estimated, as it dominates at higher frequencies. Even with just about a day of data, its value can be meaningfully constrained. On the other hand, the TMI noise remains poorly characterized, as its contribution lies in a frequency range where it is overshadowed by the stronger OMS noise, as discussed in section VII.2. Otherwise, in other words, the OMS noise that enters in $N_{2}$ is much smaller than the OMS noise that contributes to $N_{1}$ . Longer data segments are needed to improve the inference of TMI noise, as demonstrated in section VI, where the fit is performed on the full “Light-Spritz” dataset. Nonetheless, this initial estimate is useful as an indicator of a reasonable prior range to be used during the parameter estimation phase, where the full data set is analysed.

Appendix E Orthogonality of FFT of “Spritz”-like data.

In section VII.1, we propose to low-pass filter the “Light-Spritz” data before taking the FFT to better approximate a diagonal likelihood. To verify this approach, we produce a toy-model with simplified “Spritz”-like noise in the time domain and investigate the covariance matrix of the FFT. We use the same sampling rate as in “Spritz”, $f_{s}=$0.2\text{\,}\mathrm{Hz}$$ , and produce 500 realizations of a noise data stream of 4000 samples. The noise in this toy-model is the sum of the two components,

n_{\mathrm{toy}}=n_{\mathrm{TM}}+n_{\mathrm{OMS}},

(61)

where the noise PSDs of $n_{\mathrm{TM}}$ and $n_{\mathrm{OMS}}$ are given by $C_{\mathrm{acc}}^{\mathrm{ffd}}S_{\mathrm{TM}}$ and $C_{\mathrm{disp}}^{\mathrm{ffd}}S_{\mathrm{oms}}^{\mathrm{ISI}}$ , as described in section VII.2.

We approximate these noise shapes by adding white noise time series, shaped using repeated applications of cumsum and gradient operations. We then inject this noise into the variable $\eta_{12}$ of the 2nd generation TDI Michelson combination (see Eq. (3)), evaluated with equals arms and computed using the open source software pytdi Staab et al. (2023). Expressed in terms of delay operators, this corresponds to

n_{\mathrm{toy}}^{\mathrm{TDI}}=(1-D^{2})(1-D^{4})n_{\mathrm{toy}}.

(62)

We then consider three scenarios:

1.

Directly compute the FFT without any pre-processing
2.

Compute the FFT after applying a window function (hann window)
3.

Compute the FFT after applying a low-pass filter, as described in section VII.1.

After computing the FFT series, we estimate the covariance matrix by computing

C_{ij}=E[\mathrm{FFT}_{i}\mathrm{FFT}_{j}^{*}],

(63)

where $\mathrm{FFT}_{i}$ is the i’th frequency bin of the FFT of $n_{\mathrm{toy}}^{\mathrm{TDI}}$ , considering only positive frequencies. The expectation value is approximated by averaging over our 500 realizations.

The square root of the diagonal elements of this matrix are displayed in Fig. 29, rescaled by $\sqrt{\frac{2}{f_{s}\sum{w[i]^{2}}}}$ to approximate an ASD Heinzel et al. (2002). Here, $w[i]$ are the window function coefficients applied to the i’th sample.

The direct FFT is strongly affected by spectral leakage and deviates from the expected ASD at low frequencies, whereas the windowed FFT and filtered FFT approximate their respecitve model well. Note that the low-pass filter is accounted for in the modelled noise shape for the latter.

To study correlations between frequency bins with the three methods, we normalize $C_{ij}$ by the values on the diagonal and plot

\tilde{C}_{ij}=|C_{ij}/\sqrt{C_{ii}C_{jj}}|.

(64)

The results are displayed in Fig. 30. For all three scenarios, we display the matrix for the full frequency range on the left as well as a zoomed in version of just the first 100 frequency bins on the right. The top row of Fig. 30 is the direct FFT. It shows strong long-range correlations around the frequencies corresponding to the zeros of the TDI transfer function as well as at DC, but otherwise a narrow diagonal. The middle row is the windowed FFT. The window strongly suppresses the long-range correlations, but introduces strong correlations in a narrow band along the diagonal. Finally, the filtered FFT shows even stronger correlations at high frequencies, but is perfectly diagonal in the band up to $20\text{\,}\mathrm{mHz}$ used for our analysis. This validates the use of a diagonal likelihood model in Fig. VII.1.

Appendix F Bayes factor computation

In Fig. 31, we report three examples of the stretches N6, N24 and N29 as an illustration of the methodology to compute the Bayes factor with thermodynamic integration. We use 70 temperatures to perform the thermodynamic integration method, where $\beta=\frac{1}{kT}$ and $T$ represents the temperatures. As expected Karnesis et al. (2023), the stepping-stone algorithm produces more accurate estimates of the marginal likelihood than methods that use thermodynamic integration Xie et al. (2011); this is why in the computation shown here, the errors on the Bayes factor are higher. The y-axis shows the average log-likelihood at each temperature for the MBHB (glitch) model, and the x-axis shows the inverse of the temperature. This curve is integrated to calculate the log evidence, $\log Z$ . The difference between the log evidences of the MBHB model versus the glitch model gives the log Bayes factor, $\log\mathcal{B}_{\mathrm{MBHB}/\mathrm{glitch}}$ .

References

Amaro-Seoane et al. (2017) P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, J. Camp, C. Caprini, V. Cardoso, M. Colpi, J. Conklin, N. Cornish, C. Cutler, K. Danzmann, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, J. Gair, L. G. Bote, D. Giardini, F. Gibert, C. Grimani, H. Halloin, G. Heinzel, T. Hertog, M. Hewitson, K. Holley-Bockelmann, D. Hollington, M. Hueller, H. Inchauspe, P. Jetzer, N. Karnesis, C. Killow, A. Klein, B. Klipstein, N. Korsakova, S. L. Larson, J. Livas, I. Lloro, N. Man, D. Mance, J. Martino, I. Mateos, K. McKenzie, S. T. McWilliams, C. Miller, G. Mueller, G. Nardini, G. Nelemans, M. Nofrarias, A. Petiteau, P. Pivato, E. Plagnol, E. Porter, J. Reiche, D. Robertson, N. Robertson, E. Rossi, G. Russano, B. Schutz, A. Sesana, D. Shoemaker, J. Slutsky, C. F. Sopuerta, T. Sumner, N. Tamanini, I. Thorpe, M. Troebs, M. Vallisneri, A. Vecchio, D. Vetrugno, S. Vitale, M. Volonteri, G. Wanner, H. Ward, P. Wass, W. Weber, J. Ziemer, and P. Zweifel, “Laser interferometer space antenna,” (2017), arXiv:1702.00786 [astro-ph.IM] .
Colpi et al. (2024) M. Colpi, K. Danzmann, M. Hewitson, K. Holley-Bockelmann, P. Jetzer, G. Nelemans, A. Petiteau, D. Shoemaker, C. Sopuerta, R. Stebbins, N. Tanvir, H. Ward, W. J. Weber, I. Thorpe, A. Daurskikh, A. Deep, I. F. Núñez, C. G. Marirrodriga, M. Gehler, J.-P. Halain, O. Jennrich, U. Lammers, J. Larrañaga, M. Lieser, N. Lützgendorf, W. Martens, L. Mondin, A. P. Niño, P. Amaro-Seoane, M. A. Sedda, P. Auclair, S. Babak, Q. Baghi, V. Baibhav, T. Baker, J.-B. Bayle, C. Berry, E. Berti, G. Boileau, M. Bonetti, R. Brito, R. Buscicchio, G. Calcagni, P. R. Capelo, C. Caprini, A. Caputo, E. Castelli, H.-Y. Chen, X. Chen, A. Chua, G. Davies, A. Derdzinski, V. F. Domcke, D. Doneva, I. Dvorkin, J. M. Ezquiaga, J. Gair, Z. Haiman, I. Harry, O. Hartwig, A. Hees, A. Heffernan, S. Husa, D. Izquierdo, N. Karnesis, A. Klein, V. Korol, N. Korsakova, T. Kupfer, D. Laghi, A. Lamberts, S. Larson, M. L. Jeune, M. Lewicki, T. Littenberg, E. Madge, A. Mangiagli, S. Marsat, I. M. Vilchez, A. Maselli, J. Mathews, M. van de Meent, M. Muratore, G. Nardini, P. Pani, M. Peloso, M. Pieroni, A. Pound, H. Quelquejay-Leclere, A. Ricciardone, E. M. Rossi, A. Sartirana, E. Savalle, L. Sberna, A. Sesana, D. Shoemaker, J. Slutsky, T. Sotiriou, L. Speri, M. Staab, D. Steer, N. Tamanini, G. Tasinato, J. Torrado, A. Torres-Orjuela, A. Toubiana, M. Vallisneri, A. Vecchio, M. Volonteri, K. Yagi, and L. Zwick, “Lisa definition study report,” (2024), arXiv:2402.07571 [astro-ph.CO] .
Gair et al. (2011) J. R. Gair, A. Sesana, E. Berti, and M. Volonteri, Class. Quant. Grav. 28, 094018 (2011), arXiv:1009.6172 [gr-qc] .
Sesana et al. (2011) A. Sesana, J. Gair, E. Berti, and M. Volonteri, Phys. Rev. D 83, 044036 (2011), arXiv:1011.5893 [astro-ph.CO] .
Klein et al. (2016) A. Klein et al., Phys. Rev. D 93, 024003 (2016), arXiv:1511.05581 [gr-qc] .
Dayal et al. (2019) P. Dayal, E. M. Rossi, B. Shiralilou, O. Piana, T. R. Choudhury, and M. Volonteri, Mon. Not. Roy. Astron. Soc. 486, 2336 (2019), arXiv:1810.11033 [astro-ph.GA] .
Barausse et al. (2020a) E. Barausse, I. Dvorkin, M. Tremmel, M. Volonteri, and M. Bonetti, Astrophys. J. 904, 16 (2020a), arXiv:2006.03065 [astro-ph.GA] .
Toubiana et al. (2021) A. Toubiana, K. W. K. Wong, S. Babak, E. Barausse, E. Berti, J. R. Gair, S. Marsat, and S. R. Taylor, Phys. Rev. D 104, 083027 (2021), arXiv:2106.13819 [gr-qc] .
Chen et al. (2022) H.-Y. Chen, A. Ricarte, and F. Pacucci, (2022), arXiv:2202.04764 [astro-ph.GA] .
Fang and Yang (2023) Y. Fang and H. Yang, Mon. Not. Roy. Astron. Soc. 523, 5120 (2023), arXiv:2209.14509 [gr-qc] .
Spadaro et al. (2025) A. Spadaro, R. Buscicchio, D. Izquierdo-Villalba, D. Gerosa, A. Klein, and G. Pratten, Phys. Rev. D 111, 023004 (2025), arXiv:2409.13011 [astro-ph.HE] .
Langen et al. (2024) V. Langen, N. Tamanini, S. Marsat, and E. Bortolas, (2024), arXiv:2409.06527 [astro-ph.CO] .
Toubiana et al. (2024a) A. Toubiana, L. Sberna, M. Volonteri, E. Barausse, S. Babak, R. Enficiaud, D. I. Villalba, J. R. Gair, J. E. Greene, and H. Quelquejay Leclere, (2024a), arXiv:2410.17916 [astro-ph.GA] .
Berti et al. (2006) E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
Berti et al. (2016) E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Phys. Rev. Lett. 117, 101102 (2016), arXiv:1605.09286 [gr-qc] .
Chamberlain and Yunes (2017) K. Chamberlain and N. Yunes, Phys. Rev. D 96, 084039 (2017), arXiv:1704.08268 [gr-qc] .
Barausse et al. (2020b) E. Barausse et al., Gen. Rel. Grav. 52, 81 (2020b), arXiv:2001.09793 [gr-qc] .
Maggio et al. (2020) E. Maggio, L. Buoninfante, A. Mazumdar, and P. Pani, Phys. Rev. D 102, 064053 (2020), arXiv:2006.14628 [gr-qc] .
Bhagwat et al. (2022) S. Bhagwat, C. Pacilio, E. Barausse, and P. Pani, Phys. Rev. D 105, 124063 (2022), arXiv:2201.00023 [gr-qc] .
Corman et al. (2022) M. Corman, A. Ghosh, C. Escamilla-Rivera, M. A. Hendry, S. Marsat, and N. Tamanini, Phys. Rev. D 105, 064061 (2022), arXiv:2109.08748 [gr-qc] .
Toubiana et al. (2024b) A. Toubiana, L. Pompili, A. Buonanno, J. R. Gair, and M. L. Katz, Phys. Rev. D 109, 104019 (2024b), arXiv:2307.15086 [gr-qc] .
Pitte et al. (2024) C. Pitte, Q. Baghi, M. Besançon, and A. Petiteau, Phys. Rev. D 110, 104003 (2024), arXiv:2406.14552 [gr-qc] .
Mangiagli et al. (2022) A. Mangiagli, C. Caprini, M. Volonteri, S. Marsat, S. Vergani, N. Tamanini, and H. Inchauspé, Phys. Rev. D 106, 103017 (2022), arXiv:2207.10678 [astro-ph.HE] .
Dong-Páez et al. (2023) C. A. Dong-Páez, M. Volonteri, R. S. Beckmann, Y. Dubois, A. Mangiagli, M. Trebitsch, S. Vergani, and N. Webb, Astron. Astrophys. 676, A2 (2023), arXiv:2303.09569 [astro-ph.HE] .
Izquierdo-Villalba et al. (2024) D. Izquierdo-Villalba, A. Sesana, M. Colpi, D. Spinoso, M. Bonetti, S. Bonoli, and R. Valiante, Astron. Astrophys. 686, A183 (2024), arXiv:2401.10983 [astro-ph.GA] .
Caprini and Tamanini (2016) C. Caprini and N. Tamanini, JCAP 10, 006 (2016), arXiv:1607.08755 [astro-ph.CO] .
Tamanini et al. (2016) N. Tamanini, C. Caprini, E. Barausse, A. Sesana, A. Klein, and A. Petiteau, JCAP 04, 002 (2016), arXiv:1601.07112 [astro-ph.CO] .
Mangiagli et al. (2023) A. Mangiagli, C. Caprini, S. Marsat, L. Speri, R. R. Caldwell, and N. Tamanini, (2023), arXiv:2312.04632 [astro-ph.CO] .
LIGO Scientific Collaboration et al. (2015) LIGO Scientific Collaboration, J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. Adya, C. Affeldt, N. Aggarwal, O. D. Aguiar, A. Ain, P. Ajith, A. Alemic, B. Allen, D. Amariutei, S. B. Anderson, W. G. Anderson, K. Arai, M. C. Araya, C. Arceneaux, J. S. Areeda, G. Ashton, S. Ast, S. M. Aston, P. Aufmuth, C. Aulbert, B. E. Aylott, S. Babak, P. T. Baker, S. W. Ballmer, J. C. Barayoga, M. Barbet, S. Barclay, B. C. Barish, D. Barker, B. Barr, L. Barsotti, J. Bartlett, M. A. Barton, I. Bartos, R. Bassiri, J. C. Batch, C. Baune, B. Behnke, A. S. Bell, C. Bell, M. Benacquista, J. Bergman, G. Bergmann, C. P. L. Berry, J. Betzwieser, S. Bhagwat, R. Bhandare, I. A. Bilenko, G. Billingsley, J. Birch, S. Biscans, C. Biwer, J. K. Blackburn, L. Blackburn, C. D. Blair, D. Blair, O. Bock, T. P. Bodiya, P. Bojtos, C. Bond, R. Bork, M. Born, S. Bose, P. R. Brady, V. B. Braginsky, J. E. Brau, D. O. Bridges, M. Brinkmann, A. F. Brooks, D. A. Brown, D. D. Brown, N. M. Brown, S. Buchman, A. Buikema, A. Buonanno, L. Cadonati, J. Calderón Bustillo, J. B. Camp, K. C. Cannon, J. Cao, C. D. Capano, S. Caride, S. Caudill, M. Cavaglià, C. Cepeda, R. Chakraborty, T. Chalermsongsak, S. J. Chamberlin, S. Chao, P. Charlton, Y. Chen, H. S. Cho, M. Cho, J. H. Chow, N. Christensen, Q. Chu, S. Chung, G. Ciani, F. Clara, J. A. Clark, C. Collette, L. Cominsky, M. Constancio, Jr., D. Cook, T. R. Corbitt, N. Cornish, A. Corsi, C. A. Costa, M. W. Coughlin, S. Countryman, P. Couvares, D. M. Coward, M. J. Cowart, D. C. Coyne, R. Coyne, K. Craig, J. D. E. Creighton, T. D. Creighton, J. Cripe, S. G. Crowder, A. Cumming, L. Cunningham, C. Cutler, K. Dahl, T. Dal Canton, M. Damjanic, S. L. Danilishin, K. Danzmann, L. Dartez, I. Dave, H. Daveloza, G. S. Davies, E. J. Daw, D. DeBra, W. Del Pozzo, T. Denker, T. Dent, V. Dergachev, R. T. DeRosa, R. DeSalvo, S. Dhurandhar, M. D´ıaz, I. Di Palma, G. Dojcinoski, E. Dominguez, F. Donovan, K. L. Dooley, S. Doravari, R. Douglas, T. P. Downes, J. C. Driggers, Z. Du, S. Dwyer, T. Eberle, T. Edo, M. Edwards, M. Edwards, A. Effler, H. B. Eggenstein, P. Ehrens, J. Eichholz, S. S. Eikenberry, R. Essick, T. Etzel, M. Evans, T. Evans, M. Factourovich, S. Fairhurst, X. Fan, Q. Fang, B. Farr, W. M. Farr, M. Favata, M. Fays, H. Fehrmann, M. M. Fejer, D. Feldbaum, E. C. Ferreira, R. P. Fisher, Z. Frei, A. Freise, R. Frey, T. T. Fricke, P. Fritschel, V. V. Frolov, S. Fuentes-Tapia, P. Fulda, M. Fyffe, and J. R. Gair, Classical and Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
Acernese et al. (2015) F. Acernese et al. (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
Somiya (2012) K. Somiya (KAGRA), Class. Quant. Grav. 29, 124007 (2012), arXiv:1111.7185 [gr-qc] .
Aso et al. (2013) Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto (KAGRA), Phys. Rev. D 88, 043007 (2013), arXiv:1306.6747 [gr-qc] .
Nelemans et al. (2001) G. Nelemans, L. R. Yungelson, and S. F. Portegies Zwart, Astron. Astrophys. 375, 890 (2001), arXiv:astro-ph/0105221 .
Ruiter et al. (2010) A. J. Ruiter, K. Belczynski, M. Benacquista, S. L. Larson, and G. Williams, Astrophys. J. 717, 1006 (2010), arXiv:0705.3272 [astro-ph] .
Korol et al. (2017) V. Korol, E. M. Rossi, P. J. Groot, G. Nelemans, S. Toonen, and A. G. A. Brown, Mon. Not. Roy. Astron. Soc. 470, 1894 (2017), arXiv:1703.02555 [astro-ph.HE] .
Lamberts et al. (2019) A. Lamberts, S. Blunt, T. B. Littenberg, S. Garrison-Kimmel, T. Kupfer, and R. E. Sanderson, Mon. Not. Roy. Astron. Soc. 490, 5888 (2019), arXiv:1907.00014 [astro-ph.HE] .
Korol et al. (2022) V. Korol, N. Hallakoun, S. Toonen, and N. Karnesis, Mon. Not. Roy. Astron. Soc. 511, 5936 (2022), arXiv:2109.10972 [astro-ph.HE] .
Toubiana et al. (2024c) A. Toubiana, N. Karnesis, A. Lamberts, and M. C. Miller, Astron. Astrophys. 692, A165 (2024c), arXiv:2403.16867 [astro-ph.SR] .
Gair et al. (2010) J. R. Gair, C. Tang, and M. Volonteri, Phys. Rev. D 81, 104014 (2010), arXiv:1004.1921 [astro-ph.GA] .
Amaro-Seoane (2018) P. Amaro-Seoane, Living Rev. Rel. 21, 4 (2018), arXiv:1205.5240 [astro-ph.CO] .
Babak et al. (2017) S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D 95, 103012 (2017), arXiv:1703.09722 [gr-qc] .
Pan and Yang (2021) Z. Pan and H. Yang, Phys. Rev. D 103, 103018 (2021), arXiv:2101.09146 [astro-ph.HE] .
Seoane et al. (2024) P. A. Seoane, Y. Lin, and K. Tzanavaris, Phys. Rev. D 110, 064011 (2024), arXiv:2403.12134 [astro-ph.CO] .
Bonetti and Sesana (2020) M. Bonetti and A. Sesana, Phys. Rev. D 102, 103023 (2020), arXiv:2007.14403 [astro-ph.GA] .
Pozzoli et al. (2023) F. Pozzoli, S. Babak, A. Sesana, M. Bonetti, and N. Karnesis, Phys. Rev. D 108, 103039 (2023), [Erratum: Phys.Rev.D 110, 049903 (2024)], arXiv:2302.07043 [astro-ph.GA] .
Sesana (2016) A. Sesana, Phys. Rev. Lett. 116, 231102 (2016), arXiv:1602.06951 [gr-qc] .
Wong et al. (2018) K. W. K. Wong, E. D. Kovetz, C. Cutler, and E. Berti, Phys. Rev. Lett. 121, 251102 (2018), arXiv:1808.08247 [astro-ph.HE] .
Gerosa et al. (2019) D. Gerosa, S. Ma, K. W. K. Wong, E. Berti, R. O’Shaughnessy, Y. Chen, and K. Belczynski, Phys. Rev. D 99, 103004 (2019), arXiv:1902.00021 [astro-ph.HE] .
Toubiana et al. (2020) A. Toubiana, S. Marsat, S. Babak, J. Baker, and T. Dal Canton, Phys. Rev. D 102, 124037 (2020), arXiv:2007.08544 [gr-qc] .
Buscicchio et al. (2021) R. Buscicchio, A. Klein, E. Roebber, C. J. Moore, D. Gerosa, E. Finch, and A. Vecchio, Phys. Rev. D 104, 044065 (2021), arXiv:2106.05259 [astro-ph.HE] .
Toubiana et al. (2022) A. Toubiana, S. Babak, S. Marsat, and S. Ossokine, Phys. Rev. D 106, 104034 (2022), arXiv:2206.12439 [gr-qc] .
Buscicchio et al. (2025) R. Buscicchio, J. Torrado, C. Caprini, G. Nardini, N. Karnesis, M. Pieroni, and A. Sesana, JCAP 01, 084 (2025), arXiv:2410.18171 [astro-ph.HE] .
Auclair et al. (2023) P. Auclair et al. (LISA Cosmology Working Group), Living Rev. Rel. 26, 5 (2023), arXiv:2204.05434 [astro-ph.CO] .
Muratore et al. (2023) M. Muratore, O. Hartwig, D. Vetrugno, S. Vitale, and W. J. Weber, Phys. Rev. D 107, 082004 (2023).
Cornish and Crowder (2005) N. J. Cornish and J. Crowder, Phys. Rev. D 72, 043005 (2005), arXiv:gr-qc/0506059 .
Littenberg and Cornish (2023) T. B. Littenberg and N. J. Cornish, Phys. Rev. D 107, 063004 (2023), arXiv:2301.03673 [gr-qc] .
Katz et al. (2024) M. L. Katz, N. Karnesis, N. Korsakova, J. R. Gair, and N. Stergioulas, “An efficient gpu-accelerated multi-source global fit pipeline for lisa data analysis,” (2024), arXiv:2405.04690 [gr-qc] .
Strub et al. (2024) S. H. Strub, L. Ferraioli, C. Schmelzbach, S. C. Stähler, and D. Giardini, Phys. Rev. D 110, 024005 (2024), arXiv:2403.15318 [gr-qc] .
Deng et al. (2025) S. Deng, S. Babak, M. Le Jeune, S. Marsat, E. Plagnol, and A. Sartirana, (2025), arXiv:2501.10277 [gr-qc] .
Armano et al. (2024a) M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi, E. Castelli, A. Cavalleri, A. Cesarini, V. Chiavegato, A. M. Cruise, D. Dal Bosco, K. Danzmann, M. De Deus Silva, I. Diepholz, G. Dixon, R. Dolesi, L. Ferraioli, V. Ferroni, E. D. Fitzsimons, M. Freschi, L. Gesa, D. Giardini, F. Gibert, R. Giusteri, C. Grimani, J. Grzymisch, I. Harrison, M. S. Hartig, G. Heinzel, M. Hewitson, D. Hollington, D. Hoyland, M. Hueller, H. Inchauspé, O. Jennrich, P. Jetzer, B. Johlander, N. Karnesis, B. Kaune, N. Korsakova, C. J. Killow, J. A. Lobo, J. P. López-Zaragoza, R. Maarschalkerweerd, D. Mance, V. Martín, L. Martin-Polo, F. Martin-Porqueras, J. Martino, P. W. McNamara, J. Mendes, L. Mendes, N. Meshksar, M. Nofrarias, S. Paczkowski, M. Perreur-Lloyd, A. Petiteau, E. Plagnol, J. Ramos-Castro, J. Reiche, F. Rivas, D. I. Robertson, G. Russano, L. Sala, J. Slutsky, C. F. Sopuerta, T. Sumner, D. Texier, J. I. Thorpe, D. Vetrugno, S. Vitale, G. Wanner, H. Ward, P. Wass, W. J. Weber, L. Wissel, A. Wittchen, C. Zanoni, and P. Zweifel (LISA Pathfinder Collaboration), Phys. Rev. D 110, 042004 (2024a).
Armano et al. (2022a) M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi, E. Castelli, A. Cavalleri, A. Cesarini, V. Chiavegato, A. Cruise, D. Dal Bosco, K. Danzmann, M. De Deus Silva, I. Diepholz, G. Dixon, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, M. Freschi, L. Gesa, D. Giardini, F. Gibert, R. Giusteri, C. Grimani, J. Grzymisch, I. Harrison, M. Hartig, G. Heinzel, M. Hewitson, D. Hollington, D. Hoyland, M. Hueller, H. Inchauspé, O. Jennrich, P. Jetzer, B. Johlander, N. Karnesis, B. Kaune, N. Korsakova, C. Killow, J. Lobo, J. López-Zaragoza, R. Maarschalkerweerd, D. Mance, V. Martín, L. Martin-Polo, F. Martin-Porqueras, J. Martino, P. McNamara, J. Mendes, L. Mendes, N. Meshksar, M. Nofrarias, S. Paczkowski, M. Perreur-Lloyd, A. Petiteau, E. Plagnol, J. Ramos-Castro, J. Reiche, F. Rivas, D. Robertson, G. Russano, L. Sala, P. Sarra, J. Slutsky, C. Sopuerta, T. Sumner, D. Texier, J. Thorpe, D. Vetrugno, S. Vitale, G. Wanner, H. Ward, P. Wass, W. Weber, L. Wissel, A. Wittchen, C. Zanoni, and P. Zweifel, Physical Review D 106 (2022a), 10.1103/physrevd.106.062001.
Muratore (2021) M. Muratore, “Time delay interferometry for LISA science and instrument characterization,” (2021).
Zevin et al. (2017) M. Zevin et al., Class. Quant. Grav. 34, 064003 (2017), arXiv:1611.04596 [gr-qc] .
Armano et al. (2022b) M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi, E. Castelli, A. Cavalleri, A. Cesarini, V. Chiavegato, A. M. Cruise, D. Dal Bosco, K. Danzmann, M. De Deus Silva, I. Diepholz, G. Dixon, R. Dolesi, L. Ferraioli, V. Ferroni, E. D. Fitzsimons, M. Freschi, L. Gesa, D. Giardini, F. Gibert, R. Giusteri, C. Grimani, J. Grzymisch, I. Harrison, M. S. Hartig, G. Heinzel, M. Hewitson, D. Hollington, D. Hoyland, M. Hueller, H. Inchauspé, O. Jennrich, P. Jetzer, B. Johlander, N. Karnesis, B. Kaune, N. Korsakova, C. J. Killow, J. A. Lobo, J. P. López-Zaragoza, R. Maarschalkerweerd, D. Mance, V. Martín, L. Martin-Polo, F. Martin-Porqueras, J. Martino, P. W. McNamara, J. Mendes, L. Mendes, N. Meshksar, M. Nofrarias, S. Paczkowski, M. Perreur-Lloyd, A. Petiteau, E. Plagnol, J. Ramos-Castro, J. Reiche, F. Rivas, D. I. Robertson, G. Russano, L. Sala, P. Sarra, J. Slutsky, C. F. Sopuerta, T. Sumner, D. Texier, J. I. Thorpe, D. Vetrugno, S. Vitale, G. Wanner, H. Ward, P. Wass, W. J. Weber, L. Wissel, A. Wittchen, C. Zanoni, and P. Zweifel (LISA Pathfinder Collaboration), Phys. Rev. D 106, 062001 (2022b).
Abbott et al. (2017) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
Pankow et al. (2018) C. Pankow et al., Phys. Rev. D 98, 084016 (2018), arXiv:1808.03619 [gr-qc] .
Karnesis et al. (2023) N. Karnesis, M. L. Katz, N. Korsakova, J. R. Gair, and N. Stergioulas, Monthly Notices of the Royal Astronomical Society 526, 4814–4830 (2023).
Katz and Roberts (2023) M. Katz and J. Roberts, “mikekatz04/bbhx: New release!” (2023).
LISA LDC Team (nd) LISA LDC Team, “Lisa data challenges: Challenge 2b,” (n.d.), accessed: 2025-01-22.
London et al. (2018) L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa, X. Jiménez-Forteza, C. Kalaghatgi, F. Ohme, and F. Pannarale, Phys. Rev. Lett. 120, 161102 (2018), arXiv:1708.00404 [gr-qc] .
Katz et al. (2023) M. Katz, N. Karnesis, and N. Korsakova, “mikekatz04/eryn: first full release,” (2023).
Armano et al. (2018) M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi, E. Castelli, A. Cavalleri, A. Cesarini, A. M. Cruise, K. Danzmann, M. de Deus Silva, I. Diepholz, G. Dixon, R. Dolesi, L. Ferraioli, V. Ferroni, E. D. Fitzsimons, M. Freschi, L. Gesa, F. Gibert, D. Giardini, R. Giusteri, C. Grimani, J. Grzymisch, I. Harrison, G. Heinzel, M. Hewitson, D. Hollington, D. Hoyland, M. Hueller, H. Inchauspé, O. Jennrich, P. Jetzer, N. Karnesis, B. Kaune, N. Korsakova, C. J. Killow, J. A. Lobo, I. Lloro, L. Liu, J. P. López-Zaragoza, R. Maarschalkerweerd, D. Mance, N. Meshksar, V. Martín, L. Martin-Polo, J. Martino, F. Martin-Porqueras, I. Mateos, P. W. McNamara, J. Mendes, L. Mendes, M. Nofrarias, S. Paczkowski, M. Perreur-Lloyd, A. Petiteau, P. Pivato, E. Plagnol, J. Ramos-Castro, J. Reiche, D. I. Robertson, F. Rivas, G. Russano, J. Slutsky, C. F. Sopuerta, T. Sumner, D. Texier, J. I. Thorpe, D. Vetrugno, S. Vitale, G. Wanner, H. Ward, P. J. Wass, W. J. Weber, L. Wissel, A. Wittchen, and P. Zweifel, Phys. Rev. Lett. 120, 061101 (2018).
Armano et al. (2024b) M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi, E. Castelli, A. Cavalleri, A. Cesarini, V. Chiavegato, A. Cruise, D. Dal Bosco, K. Danzmann, M. De Deus Silva, I. Diepholz, G. Dixon, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, M. Freschi, L. Gesa, D. Giardini, F. Gibert, R. Giusteri, C. Grimani, J. Grzymisch, I. Harrison, M. Hartig, G. Heinzel, M. Hewitson, D. Hollington, D. Hoyland, M. Hueller, H. Inchauspé, O. Jennrich, P. Jetzer, B. Johlander, N. Karnesis, B. Kaune, N. Korsakova, C. Killow, J. Lobo, J. López-Zaragoza, R. Maarschalkerweerd, D. Mance, V. Martín, L. Martin-Polo, F. Martin-Porqueras, J. Martino, P. McNamara, J. Mendes, L. Mendes, N. Meshksar, M. Nofrarias, S. Paczkowski, M. Perreur-Lloyd, A. Petiteau, E. Plagnol, J. Ramos-Castro, J. Reiche, F. Rivas, D. Robertson, G. Russano, L. Sala, J. Slutsky, C. Sopuerta, T. Sumner, D. Texier, J. Thorpe, D. Vetrugno, S. Vitale, G. Wanner, H. Ward, P. Wass, W. Weber, L. Wissel, A. Wittchen, C. Zanoni, and P. Zweifel, Physical Review D 110 (2024b), 10.1103/physrevd.110.042004.
Sala (2023) L. Sala, “Residual test mass acceleration in lisa pathfinder: in-depth statistical analysis and physical sources,” (2023).
Castelli (2020) E. Castelli, “LISA Pathfinder noise performance results: disturbances in the sub-mHz frequency band and projection to LISA,” (2020).
Baghi et al. (2022) Q. Baghi, N. Korsakova, J. Slutsky, E. Castelli, N. Karnesis, and J.-B. Bayle, Phys. Rev. D 105, 042002 (2022), arXiv:2112.07490 [gr-qc] .
Shaddock et al. (2004) D. A. Shaddock, B. Ware, R. E. Spero, and M. Vallisneri, Phys. Rev. D 70, 081101 (2004).
Armstrong et al. (1999) J. W. Armstrong, F. B. Estabrook, and M. Tinto, The Astrophysical Journal 527, 814 (1999).
Muratore et al. (2020) M. Muratore, D. Vetrugno, and S. Vitale, Class. Quant. Grav. 37, 185019 (2020), arXiv:2001.11221 [astro-ph.IM] .
Shaddock et al. (2003) D. A. Shaddock, M. Tinto, F. B. Estabrook, and J. W. Armstrong, Phys. Rev. D 68, 061303 (2003).
Vallisneri (2005) M. Vallisneri, Phys. Rev. D 72, 042003 (2005).
Hartwig and Muratore (2022a) O. Hartwig and M. Muratore, Phys. Rev. D 105, 062006 (2022a).
Prince et al. (2002) T. A. Prince, M. Tinto, S. L. Larson, and J. W. Armstrong, Phys. Rev. D 66, 122002 (2002).
Heinzel et al. (2011) G. Heinzel, J. J. Esteban, S. Barke, M. Otto, Y. Wang, A. F. Garcia, and K. Danzmann, Classical and Quantum Gravity 28, 094008 (2011).
Hartwig (2021) O. Hartwig, “Instrumental modelling and noise reduction algorithms for the Laser Interferometer Space Antenna,” (2021).
Bergé et al. (2019) J. Bergé, R. Massey, Q. Baghi, and P. Touboul, Mon. Not. Roy. Astron. Soc. 486, 544 (2019), arXiv:1903.05837 [astro-ph.IM] .
Wikipedia contributors (nd) Wikipedia contributors, “Laguerre polynomials,” (n.d.), accessed: 2025-01-22.
Pitte et al. (2023) C. Pitte, Q. Baghi, S. Marsat, M. Besançon, and A. Petiteau, Phys. Rev. D 108, 044053 (2023), arXiv:2304.03142 [gr-qc] .
Marsat and Baker (2018) S. Marsat and J. G. Baker, “Fourier-domain modulations and delays of gravitational-wave signals,” (2018), arXiv:1806.10734 [gr-qc] .
Marsat et al. (2021) S. Marsat, J. G. Baker, and T. D. Canton, Physical Review D 103 (2021), 10.1103/physrevd.103.083011.
Hartwig and Bayle (2021) O. Hartwig and J.-B. Bayle, Phys. Rev. D 103, 123027 (2021).
Hartig et al. (2022) M.-S. Hartig, S. Schuster, and G. Wanner, Journal of Optics 24, 065601 (2022).
Chwalla et al. (2020) M. Chwalla, K. Danzmann, M. D. Álvarez, J. E. Delgado, G. Fernández Barranco, E. Fitzsimons, O. Gerberding, G. Heinzel, C. Killow, M. Lieser, M. Perreur-Lloyd, D. Robertson, J. Rohr, S. Schuster, T. Schwarze, M. Tröbs, G. Wanner, and H. Ward, Phys. Rev. Appl. 14, 014030 (2020).
Quang Nam et al. (2023a) D. Quang Nam, J. Martino, Y. Lemière, A. Petiteau, J.-B. Bayle, O. Hartwig, and M. Staab, Phys. Rev. D 108, 082004 (2023a).
Quang Nam et al. (2023b) D. Quang Nam, J. Martino, Y. Lemière, A. Petiteau, J.-B. Bayle, O. Hartwig, and M. Staab, Physical Review D 108 (2023b), 10.1103/physrevd.108.082004.
working group (2022) L. D. C. working group, “Lisa data challenge software,” (2022).
Tinto and Dhurandhar (2021) M. Tinto and S. Dhurandhar, Living Rev Relativ 24 (2021), 10.1007/s41114-020-00029-6.
Babak et al. (2021) S. Babak, M. Hewitson, and A. Petiteau, “Lisa sensitivity and snr calculations,” (2021), arXiv:2108.01167 [astro-ph.IM] .
Sivia and Skilling (2006) D. Sivia and J. Skilling, Data Analysis: A Bayesian Tutorial, 2nd ed. (Oxford University Press, 2006).
Metropolis et al. (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953).
Hastings (1970) W. K. Hastings, Biometrika 57, 97 (1970).
Hitchcock (2003) C. Hitchcock, Contemporary Debates in Philosophy of Science (Blackwell Publishing, 2003).
Gilks et al. (1995) W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman & Hall/CRC, 1995).
Green (1995) P. J. Green, Biometrika 82, 711 (1995).
Swendsen and Wang (1986) R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986).
Hukushima and Nemoto (1996) K. Hukushima and K. Nemoto, Journal of the Physical Society of Japan 65, 1604 (1996), https://doi.org/10.1143/JPSJ.65.1604 .
Vousden et al. (2015) W. D. Vousden, W. M. Farr, and I. Mandel, Monthly Notices of the Royal Astronomical Society 455, 1919 (2015).
(108) M. Katz, “Eryn documentation,” https://mikekatz04.github.io/Eryn/html/index.html, accessed: 2025-05-23.
Goodman and Weare (2010) J. Goodman and J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).
Katz (2022) M. L. Katz, Phys. Rev. D 105, 044055 (2022).
Muratore et al. (2024) M. Muratore, J. Gair, and L. Speri, Phys. Rev. D 109, 042001 (2024).
Littenberg et al. (2020) T. B. Littenberg, N. J. Cornish, K. Lackeos, and T. Robson, Physical Review D 101 (2020), 10.1103/physrevd.101.123021.
Spadaro et al. (2023) A. Spadaro, R. Buscicchio, D. Vetrugno, A. Klein, D. Gerosa, S. Vitale, R. Dolesi, W. J. Weber, and M. Colpi, Physical Review D 108 (2023), 10.1103/physrevd.108.123029.
Gair et al. (2023) J. R. Gair et al., Astron. J. 166, 22 (2023), arXiv:2212.08694 [gr-qc] .
Kwok et al. (2022) J. Y. L. Kwok, R. K. L. Lo, A. J. Weinstein, and T. G. F. Li, Phys. Rev. D 105, 024066 (2022), arXiv:2109.07642 [gr-qc] .
Gupta et al. (2025) A. Gupta, K. G. Arun, E. Barausse, L. Bernard, E. Berti, S. A. Bhat, A. Buonanno, V. Cardoso, S. Y. Cheung, T. A. Clarke, S. Datta, A. Dhani, J. M. Ezquiaga, I. Gupta, N. Guttman, T. Hinderer, Q. Hu, J. Janquart, N. Johnson-McDaniel, R. Kashyap, N. V. Krishnendu, P. D. Lasky, A. Lundgren, E. Maggio, P. Mahapatra, A. Maselli, P. Narayan, A. B. Nielsen, L. K. Nuttall, P. Pani, L. Passenger, E. Payne, L. Pompili, L. Reali, P. Saini, A. Samajdar, S. Tiwari, H. Tong, C. Van Den Broeck, K. Yagi, H. Yang, N. Yunes, and B. S. Sathyaprakash, SciPost Physics Community Reports (2025), 10.21468/scipostphyscommrep.5.
LDC (2022) “Lisa data challenge: Spritz,” (2022).
Husa et al. (2016) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. J. Forteza, and A. Bohé, Physical Review D 93 (2016), 10.1103/physrevd.93.044006.
Khan et al. (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. J. Forteza, and A. Bohé, Physical Review D 93 (2016), 10.1103/physrevd.93.044007.
Dhani et al. (2024) A. Dhani, S. Völkel, A. Buonanno, H. Estelles, J. Gair, H. P. Pfeiffer, L. Pompili, and A. Toubiana, “Systematic biases in estimating the properties of black holes due to inaccurate gravitational-wave models,” (2024), arXiv:2404.05811 [gr-qc] .
Hartwig and Muratore (2022b) O. Hartwig and M. Muratore, Physical Review D 105 (2022b), 10.1103/physrevd.105.062006.
Muratore et al. (2022) M. Muratore, D. Vetrugno, S. Vitale, and O. Hartwig, Phys. Rev. D 105, 023009 (2022).
Adams and Cornish (2010) M. R. Adams and N. J. Cornish, Phys. Rev. D 82, 022002 (2010).
LDC (2020) “Lisa data challenge: Sangria,” (2020).
Bayle et al. (2022) J.-B. Bayle, O. Hartwig, A. Petiteau, and M. Lilley, “Lisanode,” (2022).
Maturana-Russel et al. (2019) P. Maturana-Russel, R. Meyer, J. Veitch, and N. Christensen, Physical Review D 99 (2019), 10.1103/physrevd.99.084006.
Davies et al. (2024) G. C. Davies, I. Harry, M. J. Williams, D. Bandopadhyay, L. Barack, J.-B. Bayle, C. Hoy, A. Klein, H. Middleton, C. J. Moore, L. Nuttall, G. Pratten, A. Vecchio, and G. Woan, “Premerger observation and characterization of massive black hole binaries,” (2024), arXiv:2411.07020 [hep-ex] .
Cornish (2022) N. J. Cornish, Phys. Rev. D 105, 044007 (2022), arXiv:2110.06238 [gr-qc] .
Staab et al. (2023) M. Staab, J.-B. Bayle, and O. Hartwig, “Pytdi,” (2023).
Heinzel et al. (2002) G. Heinzel, A. Rüdiger, and R. Schilling, Max Plank Inst 12 (2002).
Xie et al. (2011) W. Xie, P. O. Lewis, Y. Fan, L. Kuo, and M.-H. Chen, Systematic Biology 60, 150 (2011).