Noisy pipelines II: Non-parametric Noise Inference for LISA

Introduction

Similar to existing ground-based gravitational-wave (GW) detectors, the raw data from future space-based detectors, such as the Laser Interferometer Space Antenna (LISA, Audley et al. 2017; Colpi et al. 2024), are expected to be dominated by fundamental, technical¹ and unknown noise sources. Once the dominant technical noise sources are suppressed, making LISA source (or, in part of its band, confusion) limited, data analysis methods, such as searching for signals (see e.g. Bandopadhyay & Moore 2023; Weaving et al. 2024; Bandopadhyay & Moore 2025; Cabourn Davies et al. 2025; Strub et al. 2026; Cabourn Davies & Harry 2026) and inferring their properties (see e.g. Vecchio 2004; Berti et al. 2005; Lang & Hughes 2006; Arun et al. 2009; Cornish 2022; Katz 2022; Littenberg & Cornish 2023; Klein et al. 2022; Pratten et al. 2023; Hoy & Nuttall 2024; Hoy et al. 2024), rely on being able to statistically differentiate the data from the expected noise distribution (Finn 1992; Veitch & Vecchio 2008; Veitch et al. 2015; Thrane & Talbot 2019). Having a well understood noise model is crucial to ensure accurate astrophysical parameter estimation, concrete discoveries of new physics, and to maximise the scientific return of these space-based missions.

Time delayed interferometry (TDI, Tinto & Armstrong 1999) is a post-processing technique to suppress one of LISA’s most significant noise sources: laser frequency noise. This technical noise source is orders of magnitude larger than any expected GW signal, making GW detection challenging. TDI is a technique that combines the data from LISA’s three individual spacecraft to form virtual interferometric observables. Numerous observables have been proposed in the past with different underlying assumptions, but the second-generation TDI observables: $X$, $Y$ and $Z$, allow for linearly changing arm lengths, making them the most applicable for realistic LISA scenarios (Armstrong et al. 1999; Estabrook et al. 2000; Vallisneri 2005; Tinto & Dhurandhar 2021; Tinto et al. 2023).

Even after TDI has been performed, estimating LISA’s noise distribution is more challenging than its ground-based counterpart. This is because unlike the existing LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2015), and KAGRA (Akutsu et al. 2021) GW detectors, the $X$, $Y$ and $Z$ TDI observables are correlated in their noise properties. As such, the cross spectral densities (CSDs) can no longer be assumed zero, and must be modelled alongside the power spectral densities (PSDs). A standard fix is to linearly combine the $X$, $Y$ and $Z$ observables, forming a set of approximately uncorrelated TDI observables for analyses. However, these $A$, $E$ and $T$ channels are only uncorrelated under idealised assumptions: equal arm lengths and identical noise levels in every test mass and optical link (Prince et al. 2002). It is therefore important to develop methods beyond those used in ground-based GW detectors to estimate LISA’s correlated noise distribution.

Several techniques have already been proposed to estimate LISA’s noise distribution. These can broadly be divided into modelled, non-parametric and Welch-based approaches. Katz et al. (2025); Deng et al. (2025) modelled the instrumental noise as the sum of the two fundamental LISA noise sources – the Test Mass Acceleration (TMA) and Optical Metrology System (OMS) noise – each with a fixed power-law and a free amplitude, together with a five-parameter phenomenological model for the confusion foreground from unresolved Galactic binaries. In contrast, Littenberg & Cornish (2023) adapted the non-parametric BayesLine algorithm (Littenberg & Cornish 2015) to fit the LISA noise with an Akima spline, using a reference noise spectrum that encodes important features in the LISA response. Baghi et al. (2023); Burke et al. (2026a) both used flexible splines that were transformed to the TDI observables through given transfer functions, and Aimen et al. (2026) proposed using penalised B-splines (P-splines) to model the noise from a parametric guess. Vajpeyi et al. (2026) later extended this work to model the noise with P-splines assuming no knowledge of the LISA detector. Finally, Strub et al. (2024) uses a combination of Welch PSDs (Welch 1967) and a spline to fit the LISA noise.

Since flexible non-parametric models are regularly used throughout GW astronomy, we present a non-parametric noise model for LISA that requires no prior information about the LISA observatory. We show that our model can fit the full noise covariance matrix – both the PSDs and CSDs for correlated data – for the $X$, $Y$ and $Z$ TDI observables. Owing to its flexibility, we envision that this model can be used as part of the global fit, as well as in combination with other modelled approaches to identify possible deviations from the fundamental noise sources arising from unknown effects. Further investigations are needed to identify how the model performs on e.g. different data durations (see e.g. Vajpeyi et al. 2026), including gaps in the data (see e.g. Burke et al. 2026b).

This paper is structured as follows: in § 2 we discuss the expected noise sources for the LISA observatory and in § 3 we discuss different techniques for modelling the LISA’s noise distribution, including the new non-parametric noise model described in this work. In § 4 we review the Bayesian methods used to infer the noise model properties and in § 5 we discuss our results on simulated LISA data. Finally we conclude in § 6.

LISA noise sources

LISA will consist of three spacecraft, each containing a free falling test mass. Lasers will be relayed between each of the spacecraft, and the distance between free falling test masses will be monitored over time.

LISA’s fundamental noises are well known, and their frequency dependence well documented: TMA noise limits the low-frequency sensitivity and follows a powerlaw with displacement PSD $\propto f^{-4}$; OMS noise limits the high-frequency sensitivity and follows a broken powerlaw with a displacement PSD $\sim1 + (2\mathrm{mHz}/f)^4$ (Robson et al. 2019). TMA noise arises from stray forces acting on the free-falling test masses inside each spacecraft, such as random, position independent forces and couplings to the noisy motion of the spacecraft shield (Bortoluzzi et al. 2004b). The TMA noise becomes dominant below $\sim1\mathrm{mHz}$, as demonstrated by LISA Pathfinder (Armano et al. 2018). The high frequency sensitivity is controlled by the OMS noise, which measures the distance of the free-falling test masses with respect to the reference test mass (Bortoluzzi et al. 2004a). The OMS noise is fundamentally limited by photon shot noise. However, for LISA it dominates above a few $\mathrm{mHz}$, where the TMA noise becomes subdominant (see e.g. Bayle & Hartwig 2023). Other technical noise sources, such as tilt-to-length coupling noise, clock noise, laser modulation and sideband noise will also affect LISA (Hartwig et al. 2022; Chwalla et al. 2020; Hartwig & Bayle 2021), but for this work, we only consider OMS and TMA secondary noise sources.

The major noise contribution of the LISA observatory is laser frequency noise. This will be several orders of magnitude larger than any change in arm length due to expected GW. The laser frequency noise can be suppressed by performing an additional post-processing step known as TDI. TDI is a technique that combines data from the three spacecraft into virtual interferometric observables, where the laser frequency noise is suppressed by several orders of magnitude. The second-generation TDI observables, $X$, $Y$ and $Z$, allowed for linearly changing arm lengths, making them the most applicable for realistic LISA scenarios. For this reason, we only consider the $X$, $Y$ and $Z$ variates in this work, but this brings about additional challenges as they are correlated in their noise properties.

LISA noise models

Numerous techniques have been proposed to model the noise for the LISA observatory (Littenberg & Cornish 2023; Baghi et al. 2023; Strub et al. 2024; Katz et al. 2025; Deng et al. 2025; Aimen et al. 2026; Burke et al. 2026a). A Welch estimate is the most common technique in astronomy (see e.g. Allen et al. 2012), where the timeseries data are divided into smaller, overlapping segments and the individual periodograms are averaged to reduce the variance and noise (Welch 1967). However, a Welch estimate is non-ideal for LISA as the noise will be non-stationary with many overlapping signals. Unlike in ground-based GW detectors, LISA will not have a signal-free segment of data to estimate the noise distribution off-source, and LISA’s arm lengths will breathe by a few percent a year. This means that the nulls in the transfer functions in particular will drift with time (Tinto et al. 2004; Tinto & Dhurandhar 2014; Bayle & Hartwig 2023; Burke et al. 2026a).

Various modelled approaches have been proposed in the past. In Burke et al. (2026a), we proposed a technique that took advantage of the known geometry of the LISA observatory, specifically the light travel times between individual spacecraft. This technique modeled the expected funda mental LISA noises with splines and used the given transfer functions to project the single-link measurements into TDI observables. This assumes that there are no other noise sources beyond the OMS and TMA noise; the splines have the freedom to fit additional noise sources, but the inferred TMA and OMS noise measurements will be biased. Once LISA flies, this underlying assumption may not be valid, especially for long duration data where glitches, and potentially instrumental malfunctions are more likely. For example, although ground-based GW detectors have been built and studied in detail for decades, there still exists unknown noise sources that affect the sensitivity of the instrument especially at low frequencies, see e.g. Capote et al. (2025). Similarly, in more realistic LISA scenarios, the secondary noise spectra become more complex, significantly increasing the the dimensionality of the problem, and number of splines to sample over (Quang Nam et al. 2023; Hartwig et al. 2023). For this reason, there is motivation to develop an unmodeled approach that relaxes these assumptions.

Owing to their flexibility, non-parametric approaches that estimate the noise spectrum are routinely used for ground-based GW astronomy. For example, all GW detections to date have used non-parametric approaches to estimate the noise spectrum around the time of the event (Abac et al. 2026a,b). The BayesLine (Littenberg & Cornish 2015) algorithm is typically used, which independently fits the noise from each detector by using a two-component non-parametric model composed of a cubic spline to fit the broadband noise spectrum and a linear sum of Lorentzians to fit the narrowband features, such as calibration lines. It uses a reversible jump Markov-Chain Monte-Carlo (RJMCMC) algorithm to determine the number of control points in the spline and the number of Lorentzians to best fit the data. Non-parametric approaches have not been regularly proposed for LISA global fit applications. An exception is GLASS (Littenberg & Cornish 2023). GLASS uses a simplified version of the BayesLine algorithm, which estimates the noise with an Akima spline (Akima 1970). It does not use a cubic spline as it was found to be prone to oscillations between control points, which was amplified from the large dynamic range and steeply changing spectrum at low frequencies for LISA. Unlike for ground-based GW applications, GLASS does not use additional model components to fit features beyond the base spline as a reference noise spectrum that encodes e.g. the TDI transfer function nulls is used.

The GW likelihood

Accurately describing the noise is crucial to find GW signals and avoid biased parameter estimates with misrepresented uncertainties. Similar to inferring the properties of e.g. a massive binary black hole in LISA (see e.g. Katz 2022; Hoy & Nuttall 2024), a noise model can be fit to the observed data through a process known as Bayesian inference. Although non-Bayesian methods exist, for example Welch’s method or maximum likelihood estimates (Burke et al. 2026a), Bayesian methods are preferred as they provide the full posterior distribution of the parameters rather than a single point estimate, allowing uncertainties and correlations to be quantified consistently.

In general, Bayesian inference is the process of estimating the parameters $\theta$ from some data $d$ given a model $m$. The parameters are represented by the posterior probability distribution, $P(\theta|d,m)$, and calculated through Bayes theorem,

where $\mathcal{L}(d|\theta,m)$ is the likelihood, $\Pi(\theta|m)$ is the prior and $\zeta= \int\mathcal{L}(d|\theta,m)\Pi(\theta|m)\,d\theta$ is the evidence. Assuming the noise is Gaussian and stationary, the GW log likelihood takes the Whittle form (Wainstein & Zubakov 1970; Veitch & Vecchio 2008; Veitch et al. 2015),

where $\langle a|b\rangle$ is the inner product between two frequency series (Finn 1992),

and $S(f)$ is the known one-sided PSD of the noise and $\tilde{a}^{*}$ is the complex conjugate of $\tilde{a}$. In practice, the data are sampled at $N$ discrete frequencies, and the inner product is simplified to,

When the PSD is unknown, it can also be estimated through Bayesian inference (see e.g. Littenberg & Cornish 2015; Cornish & Littenberg 2015). For example, assuming a noise model $\mathcal{M}$ with parameters $\lambda$, the posterior distribution for the PSD can be inferred alongside the parameters of the observed GW signal through Bayes theorem,

Although not necessary, we have assumed that the joint prior is independent and can be factorised into individual signal and noise components. The discrete likelihood in Eq. (2) becomes,

and the discrete inner product in Eq. (9) is,

The noise covariance matrix

Although $S$ is often treated as a function of frequency, it is more precisely (up to normalisation) the noise covariance matrix in the frequency domain (see e.g. Edy et al. 2021; Burke et al. 2025, 2026b). For stationary noise, the different frequency bins are uncorrelated, and the matrix is diagonal and real valued,

where $N$ is the number of frequency bins and $S_N(\lambda)$ is the PSD at $N$th frequency bin. In this form, Eq. (4) is a matrix multiplication and can be written as,

where $\tilde{\mathbf{a}} = (\tilde{a}_1, \ldots, \tilde{a}_N)^{\mathrm{T}}$ and similarly for $\tilde{\mathbf{b}}$.

When the noise is non-stationary, such that different frequency bins are correlated with one another, the covariance matrix has off-diagonal terms and is complex valued (with the complex terms encoding information about time, see e.g. Cornish 2026; Burke et al. 2025),

where $S_{kk^{\prime}}$ is the covariance between the $k$th and $k^{\prime}$th frequency bins.

Up until now, we have assumed a single GW detector. For a network of $N_{\mathrm{d}}$ detectors with correlated noise, each element of Eq. (10) becomes an $N_{\mathrm{d}} \times N_{\mathrm{d}}$ matrix. In general,

where $S_{kk^{\prime}}^{IJ}$ is the covariance between the $k$th frequency bin of detector $I$ and the $k^{\prime}$th frequency bin of detector $J$, with $S_{k^{\prime}k}^{JI} = (S_{kk^{\prime}}^{IJ})^{*}$. For LISA, $N_{\mathrm{d}} = 3$, representing the 3 TDI observables. From now on, we will refer to the network of detectors $N_{\mathrm{d}}$, as the three TDI channels since we are focusing on the LISA observatory in this work.

If the noise is assumed to be jointly stationary across the three TDI channels, different frequency bins are uncorrelated and only the blocks with $k = k^{\prime}$ are non-zero, with correlations between channels at the same frequency retained. $S_{kk}^{II}$ is then the PSD of the individual detector $I$, and the off-diagonal elements $S_{kk}^{IJ}$ are the CSDs between detectors $I$ and $J$. In this work, we generate non-stationary data, with correlated noise properties across the different TDI channels. However, we assume that the noise is locally non-stationary in our analyses, simplifying the noise covariance matrix. The validity of this assumption was explored in Burke et al. (2026a). We leave the full non-stationary extension for future work.

Stochastic sampling

Although the GW likelihood is well known, it is difficult to analytically compute the posterior distribution. As such, stochastic sampling algorithms, such as nested sampling (Skilling 2004, 2006) and Markov-Chain Monte-Carlo (Metropolis & Ulam 1949), are often used to draw samples from the unknown posterior distribution; although see (Pankow et al. 2015; Lange et al. 2018; Gabbard et al. 2022; Green & Gair 2021; Green et al. 2020; Dax et al. 2021; Tiwari et al. 2023; Fairhurst et al. 2023; Dax et al. 2025; Raymond et al. 2025; Roulet et al. 2024; Mushkin et al. 2025; Williams et al. 2025; Roulet et al. 2026) for alternative approaches. Unfortunately, classical stochastic samplers scale poorly with dimension. This can make it difficult to sufficiently explore high-dimensional distributions, especially with narrow, multimodal or strongly correlated posterior distributions. As a result, sampling over the full correlated likelihood in Eq. (6) is difficult without rigorous sampler settings; for example, assuming each term in the noise covariance matrix is modelled with 20 parameters, the full likelihood involves sampling over a 120 dimensional parameter space. One way to improve sampler convergence is to reduce the dimensionality of the likelihood. This can be achieved by considering a blocked approach.

To reduce the computational challenge of stochastically sampling a large parameter space, we can approximate Eq. (6) by employing a pairwise likelihood, in which the noise covariance matrix is split into several $2 \times2$ blocks. For three TDI channels, $X$, $Y$ and $Z$, the $3 \times3$ noise covariance matrix can be replaced by 3, $2 \times2$ matrices,

which can be solved in parallel. The computational cost can be further reduced by first inferring the PSDs of the individual channels, $S^{XX}$, $S^{YY}$ and $S^{ZZ}$, independently and in parallel. The CSD for each pair $(I, J)$ can then inferred with the PSDs $S^{II}$ and $S^{JJ}$ fixed at their best fit values. For three TDI channels, 3 analyses can be launched in parallel to model the PSDs, and once complete, a further 3 analyses can be launched in parallel to model the CSDs.

This blocked approach significantly reduces the size and degeneracies in the parameter space, as each analysis is now sampling over a e.g. 20 dimensional parameter space. This significantly improves sampler convergence and reduces the wall-time of analyses. This is at the cost of only providing a single estimate for the PSDs of each detector, and therefore the uncertainties on these measurements are neglected when estimating the CSDs. Also, by using this blocked approach, information about e.g. $S^{XY}$ is neglected when considering $S^{XZ}$, and positive definiteness of the full $3 \times3$ matrix is not enforced.

To check when the blocked approximation breaks down, we explicitly compared the blocked approach against the full likelihood for different combinations. We found that, as expected, when at most one pair of channels is coherent, the blocked likelihood matches the full likelihood exactly. When two or more pairs are coherent, as is the case for the $X$, $Y$ and $Z$ channels, the blocked likelihood diverges from the full likelihood. By fitting each cross-spectrum independently, the blocked likelihood neglects the coupling between the coherencies, and the discrepancy grows as the coherency increases. Nevertheless, we believe the blocked approach has merit for modelling the noise in LISA. For example, the blocked approximation holds true and will significantly reduce the computational cost when analysing the $A$, $E$ and $T$ channels (which is common, see e.g. (Katz 2022; Buscicchio et al. 2021; Pratten et al. 2023; Weaving et al. 2024; Hoy & Nuttall 2024; Cabourn Davies et al. 2025; Katz et al. 2025; Sharma et al. 2025)), which are quasi-orthogonal in second generation TDI and therefore have small coherences.

The PSD model

Similar to ground-based GW applications, we propose using a two-component model to fit the PSDs of the LISA observatory: one component to fit the broadband noise spectrum and another to fit high frequency regions where the noise formally goes to zero due to the TDI transfer functions. The latter is crucial due to the geometry of the LISA detector: the noise drops to exactly zero for signals with wavelength that fit an integer number of cycles within the detector arm.

We model the broadband noise spectrum with an Akima spline, $g(f)$, and the high frequency features with a linear combination of inverse Lorentzians. The Akima spline is controlled by $N + N_g$ knots, where $N$ is the number of knots within the frequency band of interest and $N_g$ are a limited number of ghost knots to control edge effects.

In second generation TDI, the nulls in the TDI transfer function occur due to a $\sin^2(2\pi fL/c)\sin^2(4\pi fL/c)$ term in the transfer function (Armstrong et al. 1999; Estabrook et al. 2000; Vallisneri 2005; Tinto & Dhurandhar 2021; Tinto et al. 2023), where $L$ is the arm length of the observatory. At odd multiples of 30 mHz only the $\sin^2(4\pi fL/c)$ term vanishes, whereas at multiples of 60 mHz both terms vanish, producing a deeper and sharper feature; a fourth order zero for the 60 mHz null compared to a second order zero at the 30 mHz and 90 mHz nulls. Owing to this, we use two different inverse Lorentzian models for the different structures: a locally quadratic model for the 30 mHz and 90 mHz nulls, and a locally quartic model for the 60 mHz null. In all cases, the inverse Lorentzian, $L(f)$, is parametrised with a central frequency $f_0$, a quality factor $\gamma$, and a minimum value at the resonant frequency $S_{\min}$,

where $\mathcal{N} = 2$ for the 30 mHz and 90 mHz nulls, and $\mathcal{N} = 4$ for the 60 mHz null.

We show an example of our two-component model in Figure 1. The full PSD is the sum of the base Akima spline and $\mathcal{N}$ inverse Lorentzian fits where $\mathcal{N} = 3$ in this work due to our chosen high frequency cut-off. To smoothly transition between the base Akima spline and the inverse Lorentzians, we apply a smooth Sine-squared taper, enabling continuity at the connection frequencies $f_{\mathrm{low}}$ and $f_{\mathrm{high}}$.

The CSD model

For the CSDs, we assume the cross spectrum is a combination of real and imaginary components. Rather than fitting $S^{IJ}_{kk'}$ directly, as in the PSD case, we model $S^{IJ}_{kk'} = \rho_{IJ}\sqrt{S^{II}_{kk'}S^{JJ}_{kk'}}$ where $\rho_{IJ}$ is defined as the complex coherency between detectors $I$ and $J$. Under this parameterisation, $\rho_{IJ}$ quantifies how strongly the two detectors are linked, and ranges between 0, uncorrelated, to 1, perfectly correlated. For the model to be valid, $\lvert\rho_{IJ} \rvert< 1$.

To account for both the real and imaginary components while enforcing $\lvert\rho_{IJ}\rvert< 1$, we define $\rho_{IJ} = \kappa r$, where $r = g_{\mathrm{re}} + i g_{\mathrm{im}}$ and,

In this form, we are free to model the real and imaginary components of $r$ independently with individual Akima splines.

In Burke et al. (2026a), we demonstrated that unbiased and accurate noise inference can be obtained for the CSDs when modelling only the real component; neglecting the imaginary components widen the knot uncertainties by at most 3% even at the most asymmetric LISA constellation in a two year period. We nevertheless model both the real and imaginary components in this work to stress test the CSD model.

Results

Analysis setup

To verify the non-parametric models described in § 3, and to test the different likelihood approaches described in § 4, we analysed 30 days of data with realistic LISA orbits, split into 7 day segments. We inferred the noise model through Bayesian inference between the frequencies $10^{-4}\,\mathrm{Hz} < f < 10^{-1}\,\mathrm{Hz}$. The simulation parameters for this dataset are nearly identical to those listed in Table I of Burke et al. (2026a). We generated a noise-only dataset generated with lisainstrument (Bayle et al. 2023; Bayle & Hartwig 2023) using ESA trailing orbits described in lisaorbits (Bayle et al. 2024), and converted to the TDI observables using extended precision via pytdi (Staab et al. 2022). The data was recast using the mojito package (Bayle 2026) and post-processed using MojitoProcessor (Burke 2026). The only difference with Burke et al. (2026a), is that we elected not to downsample the data to a cadence of $\Delta t = 5\,\mathrm{s}$, as we wanted to explicitly show the efficacy of our PSD and CSD models for the sharper null at 60 mHz, which would otherwise be removed as a result of downsampling (Burke et al. 2026a). We analyse 7 days of data so we can safely assume stationarity, and hence validate our assumption that the likelihood is described by a Whittle function. In 1 day of data, LISA’s arm length is expected to vary by at most 0.035% (Hartig et al. 2025), meaning $< 0.25\%$ in 7 days of data.

To sample the posterior distribution, we used Markov-Chain Monte-Carlo (MCMC) methods. We specifically used the Eryn (Karnesis et al. 2023) sampler since it allows for Reversible-Jump parallel-tempered MCMC, where the dimensionality of the model can change during the sampling (Reversible-Jump), and multiple tempered chains can be launched in parallel to avoid local maxima (parallel-tempering). This is useful for e.g. dynamically choosing the number of knots in the base Akima splines of our models. We performed a small convergence study on the 1 day of stationary Gaussian noise coloured by the SciRDv1 noise PSD (Babak et al. 2021) to identify sufficient sampler settings for robust results in a large prior volume. We settled on 800 walkers, 8 parallel temperatures, 2000 steps, and 5000 samples discarded as burn-in. These settings are comparable to those used in previous works with the Eryn sampler (Puecher et al. 2026). We found that these settings reduced the variance in median inferred PSD, as well as reduced the inferred 90% credible interval when using the blocked likelihood. We leave a detailed sampler convergence study for the full and blocked likelihoods for future work. Unless otherwise stated, we seeded the MCMC walkers with a random point chosen from the prior.

For both the PSD and CSD models described in § 3, we placed 8 knots uniformly in log frequency within the window: $0.05\ \mathrm{mHz} < f < 100\ \mathrm{mHz}$, and 6 ghosts knots, 3 either side of the window placed between $0.001\ \mathrm{mHz} < f < 0.05\ \mathrm{mHz}$ and $100\ \mathrm{mHz} < f < 1\ \mathrm{Hz}$ respectively to control edge effects. The number of knots were chosen based on a preliminary investigation, where an RJMCMC analysis was performed on the broadband noise spectrum consisting of only the analytic TMA and OMS noise. When allowing the number of knots to vary, we found a small preference for 8 knots. Future work will allow the number of knots to vary during each MCMC analysis. In both the PSD and CSD model, we assumed a log-uniform distribution for the knot amplitudes between $10^{-47}\ \mathrm{Hz}^{-1} < A < 10^{-37}\ \mathrm{Hz}^{-1}$. Finally, for the 3 inverse Lorentzians used in the PSD model, we assumed a uniform prior for the central frequency $f_0$, and uniform distributions for the lower and upper connection frequencies $f_{\mathrm{low}}$ and $f_{\mathrm{high}}$ to the base spline. Owing to the different functional forms of the inverse Lorentzians, we employed different priors for the quality factors: we assumed a uniform distribution between 10 and 30 for the 30 and 90 mHz nulls and a uniform distribution between 2 and 30 for the 60 mHz null. In our analyses, we do not allow for the “lack of TDI null” case, where the quality factor is 0 by definition. We purposely use wide priors for the knot amplitudes to remain as agnostic as possible. Smaller and more tailored priors could be used for each knot amplitude, for example based on a Welch estimate, however, this would mean that the same data is used twice: once for informing the prior, and once for obtaining the posterior. Posterior repartitioning (Chen et al. 2019; Prathaban et al. 2026) offers a solution to using the same data twice, but this is not explored in this work.

Inferred PSDs and CSDs

In Figure 2 we compare the inferred PSDs with the periodogram for the TDI $X$, $Y$ and $Z$ observables. When focusing on the median inferred PSD, we see that for all three channels, the blocked likelihood more closely matches the underlying periodogram: while the broadband structure of the individual PSDs from both methods are similar, the median estimate from the blocked likelihood more faithfully captures the nulls in the PSD due to the TDI transfer functions. This can be seen more clearly in Figure 3.

We believe the differences between the median inferred PSDs from the two likelihood approaches are due to poor sampler convergence when using the full likelihood. We found that the sampler struggled to find viable solutions, which is reflected in the large spread of individual PSD draws in Figure 2 and Figure 3. We see that for some draws, the nulls are captured well, but due to large uncertainties, the median estimate averages over the nulls to a smooth function. The sampler diagnostics support this interpretation. For the blocked likelihood, consecutive samples show the autocorrelation typical of an MCMC chain, where each step is a small move from the previous position. For the full likelihood, consecutive samples are effectively uncorrelated, and every proposed swap between temperatures is accepted. This indicates that the full likelihood provides little information to guide the walkers towards the high-probability region of the posterior.

The full likelihood is sampling over a 162 dimensional parameter space – 26 for each PSD (14 for the spline amplitudes and 12 for the three inverse Lorentzians) and 28 for each CSD (14 for both the real and imaginary components – while the blocked likelihood is sampling over at most 28 dimensions in a single analysis. The significant increase in the dimensionality, while keeping the sampler settings consistent, means the sampler explores a much smaller fraction of the posterior volume for the same number of samples when using the full likelihood. Scaling the autocorrelation time measured for the full likelihood with the number of dimensions, we estimate that Eryn would require a burn-in of $> 2 \times10^{4}$ steps, followed by $> 5 \times10^{4}$ steps of sampling, to obtain robust results; an order of magnitude more than used here. This highlights that alternative sampling approaches (see e.g. Williams et al. 2026) may be required if the full likelihood is used due to the significant increase in computational cost.

Looking specifically at Figure 3, we see that the combination of our PSD model and the blocked likelihood generally fits the nulls in the transfer function well. Despite wide priors, we find that the median inferred central frequency of each null is found within at most 0.3%, with the 60 mHz model performing optimally: 0.3%, 0.05% and 0.3% for the 30 mHz, 60 mHz and 90 mHz nulls respectively. Interestingly, we find that the tails of the inverse Lorentzian for the 30 mHz result slightly overestimate the noise when using the blocked likelihood. We found that this was primarily influenced by the inferred 60 mHz null: if narrow priors are chosen such that the 60 mHz null is misrepresented, the 30 mHz and 90 mHz perfectly describe the data. This implies that it may not be a problem when downsampling the data to a cadence of $\Delta t = 5\,\mathrm{s}$, where the 60 mHz null is removed due to downsampling (Burke et al. 2026a). Despite large uncertainties, with individual PSD draws varying significantly across each null, we found that the full likelihood result still infers median central frequencies of each null to within at most 1% of the truth. In general, we also see good agreement with a reference Welch estimate, both across the broadband and for each null, providing further confidence that our non-parametric fit is performing as expected.

In Figure 4, we assess how well each method whitens the data in each TDI channel. In both cases, we whitened by the median inferred PSD. If our inferred PSD is a good estimate of LISA’s noise, we expect the distribution of the whitened amplitudes should resemble a unit normal distribution (assuming the data is Gaussian and stationary). Agreement with a unit normal is a necessary, though not sufficient, indication of a good PSD estimate. As can be seen, we obtain the expected normal distribution for the blocked analysis. We see a small divergence from the expected normal distribution in the tails for the Z channel, but this will likely be resolved with stricter sampler settings. Since the median PSD inferred by the full likelihood is failing to capture the nulls due to poor convergence, we obtain a distribution with a smaller than expected unit variance. This highlights that it is not not whitening the data as well as the blocked approach.

Finally, in Figure 5 we compare the inferred CSDs with the cross periodograms for the TDI $XY$, $XZ$ and $YZ$ observables. We see a similar story as in Figure 2: the blocked likelihood more closely matches the cross periodograms. We see that the full likelihood obtains larger uncertainties, which is likely because the blocked likelihood is not accounting for errors on the PSDs. Nevertheless, if only a median estimate is used, consistent with ground-based GW detectors (Abac et al. 2026a), the blocked likelihood more closely resembles the periodogram.

Discussion

Modelling the noise distribution in the LISA observatory is crucial in order to search for GW signals, infer unbiased source properties, and maximise the scientific return of the LISA mission. Although simple methods exist (Welch 1967) that are commonly used throughout astronomy, more complex models are needed for LISA since the data will be non-stationary, with overlapping signals, and no noise-only segments of off-source data are available. In this work, we introduced a non-parametric model to fit the full and complex noise covariance matrix, with a combination of splines and inverse Lorentzians. Our model is completely agnostic and requires no prior knowledge of the LISA detector or the TDI transfer functions.

By analysing 7 days of simulated noise-only LISA data with realistic orbits, we demonstrated that our model fits the periodograms of the different TDI observables to good precision and whitens the data as expected. When comparing to Welch estimates we see good agreement, which is expected since our dataset contains no overlapping signals and any effects from non-stationarity are minimal. We emphasise that for realistic LISA scenarios, our model offers advantages beyond Welch since it can estimate the noise distribution on-source, and can be included within a global fit framework to jointly infer the astrophysical source properties and the noise model. However, we note that a fully model agnostic approach may be degenerate with LISA’s stochastic background (Adams & Cornish 2010).

We also considered two approaches to obtain the full noise covariance matrix: one where the full likelihood is split and solved in parallel with reduced dimensionality, and one where the full likelihood is solved in a single analysis. We found that although the split likelihood is an approximation that breaks down in highly correlated data, it is the only approach that gave converged results. We expect the full likelihood to give consistent results when significantly stricter settings are used. However, we note that in the majority of previous work (Katz 2022; Buscicchio et al. 2021; Pratten et al. 2023; Weaving et al. 2024; Hoy & Nuttall 2024; Cabourn Davies et al. 2025; Katz et al. 2025; Sharma et al. 2025), only the PSDs are used when modelling the noise, and the blocked approach provides robust and converged results for this case.

We envision that our non-parametric noise model can be used as part of the global fit for estimating the noise on-source, or used in conjunction with modelled approaches to identify, quantify, and constrain the frequency dependence for unknown noise artifacts once LISA flies.

Acknowledgments

We thank Sophie Bini, Neil Cornish, Peter Lott, Andrew Lundgren, Joe Metson, Julianna Ostrovska and Saulo Soares-de-Albuquerque-Filho for PSD related discussions in the context of ground-based GW observatories, which were helpful for this work. We also thank the UK LISA ground segment for comments on this work. C.H thanks the University of Portsmouth for support through the Dennis Sciama Fellowship. J.B. and O.B. acknowledge financial support from the Grant UKRI972 awarded via the UK Space Agency, and L.K.N thanks the UKRI Future Leaders Fellowship for support through the grant UKRI2745. We are also grateful for computational resources provided by the SCIAMA high performance computing cluster which is supported by the Institute of Cosmology and Gravitation (ICG) and the University of Portsmouth. For the purpose of open access, the author(s) has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

Data Availability

The analysis code behind this paper, including the non-parametric models for the PSDs and CSDs, is available through the lisa-noise package at DOI 10.5281/zenodo.23501674.

This paper has been typeset from a TEX/\LaTeX file prepared by the author.