# Modeling non stationary noise in pulsar timing array data analysis

Mikel Falxa, J Antoniadis, D.J Champion, Ismaël Cognard, G Desvignes, Lucas Guillemot, H Hu, G Janssen, J Jawor, R Karuppusamy, et al.

► To cite this version:

Mikel Falxa, J Antoniadis, D.J Champion, Ismaël Cognard, G Desvignes, et al. Modeling non stationary noise in pulsar timing array data analysis. Physical Review D, 2024, 109 (12), pp.123010. 10.1103/PhysRevD.109.123010. [https://doi.org/10.1103/PhysRevD.109.123010](https://doi.org/10.1103/PhysRevD.109.123010). [arXiv:2009.04586](https://arxiv.org/abs/2009.04586)

HAL Id: hal-04584563

https://hal.science/hal-04584563v1

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Modeling non stationary noise in pulsar timing array data analysis

Mikel Falxa,<sup>1, 2, *</sup> J. Antoniadis,<sup>3, 4</sup> D. J. Champion,<sup>4</sup> I. Cognard,<sup>1, 5</sup> G. Desvignes,<sup>4</sup> L. Guillemot,<sup>1, 5</sup>

H. Hu,<sup>4</sup> G. Janssen,<sup>6, 7</sup> J. Jawor,<sup>4</sup> R. Karuppusamy,<sup>4</sup> M. J. Keith,<sup>8</sup> M. Kramer,<sup>4</sup> K. Lackeos,<sup>4</sup> K. Liu,<sup>9, 4</sup>

J. W. McKee,<sup>10, 11</sup> D. Perrodin,<sup>12</sup> S. A. Sanidas,<sup>13</sup> G. M. Shaifullah,<sup>14, 15, 12</sup> and G. Theureau<sup>1, 5, 16</sup>

<sup>1</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace,
Université d'Orléans/CNRS, 45071 Orléans Cedex 02, France

<sup>2</sup>Dipartimento di Fisica "G. Occhialini", Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

<sup>3</sup>FORTH Institute of Astrophysics, N. Plastira 100, 70013, Heraklion, Greece

<sup>4</sup>Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>5</sup>Observatoire Radioastronomique de Nançay, Observatoire de Paris,
Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

<sup>6</sup>ASTRON, Netherlands Institute for Radio Astronomy,
Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands

<sup>7</sup>Department of Astrophysics/IMAPP, Radboud University Nijmegen,
P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

<sup>8</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Manchester M13 9PL, U

<sup>9</sup>Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China

<sup>10</sup>E. A. Milne Centre for Astrophysics, University of Hull,
Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

<sup>11</sup>Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM),
University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

<sup>12</sup>INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius (CA), Italy

<sup>13</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Manchester, M13 9PL

<sup>14</sup>Dipartimento di Fisica "G. Occhialini", Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

<sup>15</sup>INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

<sup>16</sup>Laboratoire Univers et Théories LUTh, Observatoire de Paris,
Université PSL, CNRS, Université de Paris, 92190 Meudon, France

(Dated: May 7, 2024)

Pulsar Timing Array (PTA) collaborations recently reported evidence for the presence of a gravitational wave background (GWB) in their datasets. The main candidate that is expected to produce such a GWB is the population of supermassive black hole binaries (SMBHB). Some analyses showed that the recovered signal may exhibit time-dependent properties, i.e. non-stationarity. In this paper, we propose an approximated non-stationary Gaussian process (GP) model obtained from the perturbation of stationary processes. The presented method is applied to the second data release of the European pulsar timing array to search for non-stationary features in the GWB. We analyzed the data in different time slices and showed that the inferred properties of the GWB evolve with time. We find no evidence for such non-stationary behavior and the Bayes factor in favor of the latter is $B_S^{NS} = 1.5$. We argue that the evolution of the GWB properties most likely comes from the improvement of the observation cadence with time and better characterization of the noise of individual pulsars. Such non-stationary GWB could also be produced by the leakage of non-stationary features in the noise of individual pulsars or by the presence of an eccentric single source.

# I. INTRODUCTION

Supermassive black holes (SMBHs) are incredibly massive, over a million times heavier than the sun, and are typically located at the centers of galaxies. When these black holes come together, they form binary systems that emit strong gravitational waves as they approach each other [1]. The combined effect of many such binaries produces a continuous and stationary noise signal in the nanohertz frequency band [2], known as the gravitational wave background (GWB). Pulsar timing array

(PTA) collaborations search for the induced effect of this GWB in the timing measurements of millisecond pulsars. They specifically target the characteristic spatial correlations it would exhibit between pairs of pulsars, known as the Hellings-Downs (HD) correlation pattern [3]. Recent studies have shown strong evidence for the presence of such a signal [4–6]. However, it does not exactly match what was expected for a population of supermassive black hole binaries (SMBHB) in circular orbits. In particular, it was observed that the inferred GWB properties differed when datasets of different lengths were analyzed. Several hypotheses have been proposed to explain this discrepancy, one of which is the possible non-stationarity of the GWB.

The European Pulsar Timing Array (EPTA) collab-

<footnote>* mikel.falxa@unimib.it</footnote>
oration has released its second dataset, combining the timing data of 25 pulsars over a span of up to 24.7 years of observation <fnref n="1" />. Several radio telescopes were used to collect the data: the Effelsberg Radio Telescope in Germany, the Lovell Telescope in the UK, the Nançay Radio Telescope in France, the Westerbork Synthesis Radio Telescope in the Netherlands, the Sardinia Telescope in Italy, and the Large European Array for Pulsars. In their analysis [6], the EPTA searches for the presence of a GWB in the dataset by evaluating the Bayes factor, comparing common HD correlated noise against common uncorrelated noise denoted as $B_{CURN}^{HD}$. The significance and properties of the GWB vary notably across the dataset. When analyzing the entire dataset, the Bayes factor is only $B_{CURN}^{HD} = 4$, but it is shown that using only the latest 10.33 years of data, the Bayes factor rises to $B_{CURN}^{HD} = 60$. The justification for trimming the dataset is the evolution of data quality. Legacy data have a worse cadence of observation and fewer multi-frequency band observations [7]. These two qualities are crucial for performing a good characterization of all sources of noise present in the timing data, essential for an accurate observation of the GWB [8]. Radio observatories are equipped with receiver-backend systems where the receiver collects the raw radio data that goes through a real-time processing stage in the backend systems, providing analyzable data [7]. The latest 10.33 years of data were obtained from new-generation backends with greatly improved observing abilities. However, the evolution of properties of the observed signal was not only reported by the EPTA collaboration. The Parkes Pulsar Timing Array (PPTA) performed a time-slice analysis where their dataset was analyzed in slices of 3 or 9 years and showed that the amplitude of the signal varied with time [5]. The reason behind this time-evolving property remains unexplained as the GWB is expected to be stationary [1, 2]. One advanced hypothesis is non-stationarity of noise present in the pulsars themselves that would contaminate our estimate of the signal. The two dominant sources of time-correlated noise in pulsar timing are the red noise (RN), due to stochastic variation of the pulsar's spin rate, and the dispersion measure (DM) noise, due to the fluctuation of the electron density in the interstellar medium during the path that the pulses take to Earth [8] (this noise has an amplitude proportional to $\nu^{-2}$ with $\nu$ being the incoming photon's frequency). The noise contributions are typically assumed to be stationary. Other potential sources of non-stationarity of astrophysical origin are eccentric SMBHs. They produce gravitational waves (GW) with significantly varying frequency content within one orbital period of the binary [9, 10]. The presence of one or more eccentric sources could influence a stationary GWB and make it non-stationary.

In this paper, we propose a way to model non-stationary noise. Previous works explored the idea of modeling non-stationary signals in GW analysis using wavelets [11, 12]. Here, we develop a model which is a perturbation of the Gaussian process (GP) methods that are commonly used in PTA data analysis [13]. The main motivation behind this model is to evaluate the significance of non-stationary features in noise through the evaluation of the Bayes factor by allowing time dependence of the noise spectrum properties. Strong variations in time result in strong correlations between different frequencies of the spectrum that are not accounted for by ordinary Fourier decomposition [12]. In general, the data are cut into smaller intervals where the noise is considered approximately stationary, and a Fourier transform is applied in each window to monitor the evolution of the frequency content. The assumption of local stationarity presupposes that the properties of the spectrum evolve slowly in the analyzed interval [14].

This paper is organized as follows. In the first section, we will present a non-stationary GP model. The ability of this model to correctly infer parameters is then tested on mock datasets in which non-stationary noise is injected. Finally, we will search for the possible presence of a non-stationary GWB in the EPTA dataset and argue that nearby eccentric SMBH binaries could introduce non-stationary features in the GWB.

## II. MODEL

### A. Stationary gaussian process

A stationary process is a stochastic process whose properties, typically its mean and variance, remain constant over time. In practice, we often encounter what we refer to as weak-sense stationary processes, where the mean and correlation functions remain invariant under shifts in time, i.e., they are time-translation symmetric [15]. Given the one-sided power spectral density (PSD) $S_0(f)$ of the process, the correlation function $C(t, t')$ can be expressed using the Wiener-Khinchine integral:

$$C(t, t') = \int_{0}^{\infty} df S_{0}(f) \cos(2\pi f(t - t')) \quad (2.1)$$

which is a function of the time difference $t - t'$.

Alternatively, we can define a zero-mean stationary Gaussian Process (GP) expressed as a finite Fourier sum with normally distributed coefficients [13, 16]

$$\begin{align}
n_0(t) &= \sum_{i=1}^{N} X_i \sin(2\pi f_i t) + Y_i \cos(2\pi f_i t) \tag{2.2} \\
&= \sum_{i=1}^{N} \vec{\omega}_i \cdot \vec{\phi}_i(t) \nonumber
\end{align}$$

<footnote n="1">¹The EPTA data is publicly available at https://epta.pages.in2p3.fr/epta-dr2/</footnote>
where $\vec{\phi}_i(t) = [\sin(2\pi f_i t), \cos(2\pi f_i t)]$ and $\vec{\omega}_i = [X_i, Y_i]$ with $X_i, Y_i \sim \mathcal{N}(0, \sigma_i^2)$.

The basis $\vec{\phi}_i(t)$ of cosine and sine functions is evaluated at a discrete set of frequencies $f_i$. For stationary signals, the random weights $\vec{\omega}_i$ are uncorrelated because the underlying Fourier basis of the process is orthogonal. Thus, using $\langle \vec{\omega}_i, \vec{\omega}_j \rangle = \langle \omega^2 \rangle_{ij} = \sigma_i^2 I_2 \delta_{ij}$ (where $\langle a, b \rangle$ denotes the ensemble averaged inner product between $a$ and $b$, and $I_2$ is the identity matrix of rank 2), we have:

$$ 
\begin{aligned}
C(t, t') &= \langle n(t), n(t') \rangle \\
&= \sum_{ij} \vec{\phi}_i^\top(t) \langle \omega^2 \rangle_{ij} \vec{\phi}_j(t') \\
&\equiv \sum_i \sigma_i^2 \cos(2\pi f_i (t - t'))
\end{aligned} 
\quad (2.3) $$

which is an approximation of the integral in Equation 2.1 for $\sigma_i^2 = S_0(f_i) \Delta f_i$ with $\Delta f_i = f_{i+1} - f_i$ [17].

## B. Non-stationary gaussian process

By definition, a non-stationary process exhibits spectral properties that vary with time [18, 19]. Let us consider a time-dependent evolutionary power spectral density defined as:

$$ S(f, t) = S_0(f) \times g^2(f, t), \quad (2.4) $$

where $S_0(f)$ represents a stationary PSD, and $g(f, t)$ is an arbitrary function of time that introduces perturbations to the stationary PSD.

Using this approach, we extend the definition of the zero-mean GP presented in Equation 2.2 by introducing time-varying coefficients $X_i(t)$ and $Y_i(t)$, where $X_i(t), Y_i(t) \sim \mathcal{N}(0, \sigma_i^2(t))$ and $\sigma_i^2(t) = S_0(f_i) \times g^2(f_i, t) \Delta f_i$. By applying the identity $\mathcal{N}(0, a^2\sigma^2) = a\mathcal{N}(0, \sigma^2)$, the function $g(f, t)$ factors out of the Gaussian weights, yielding:

$$ 
\begin{aligned}
n(t) &= \sum_{i=1}^{N} X_i(t_0)g(f_i, t) \sin(2\pi f_i t) + Y_i(t_0)g(f_i, t) \cos(2\pi f_i t) \\
&= \sum_{i=1}^{N} \vec{\omega}_i \cdot \vec{\Phi}_i(t)
\end{aligned} 
\quad (2.5) $$

with $\Phi_i(t) = g(f_i, t) \times \phi_i(t)$.

This expression is identical to the stationary case except the basis $\Phi_i(t)$ can now be modulated over time through the function $g(f_i, t)$. This modulation of the basis gives a description of the non-stationarity of the signal (see Figure 1). In this context, the relationship $\langle \vec{\omega}_i, \vec{\omega}_j \rangle = \langle \omega^2 \rangle_{ij}$ does not necessarily maintain a diagonal form because the weights are not expected to be

uncorrelated. However, we proceed by making an assumption of diagonality. In other words, we assume that the basis $\Phi_i(t)$ is the natural (eigen)basis of the process. In that case, the covariance matrix becomes:

$$ 
\begin{aligned}
C(t, t') &= \sum_{ij} \vec{\Phi}_i^\top(t) \langle \omega^2 \rangle_{ij} \vec{\Phi}_j(t') \\
&\approx \sum_i^N \sigma_i^2 g(f_i, t) g(f_i, t') \cos(2\pi f_i (t - t'))
\end{aligned} 
\quad (2.6) $$

This kernel can be seen as a specific instance of non-stationary generalized spectral kernels (see [20, 21] and appendix). It is simplified by considering stationary frequency components $f_i$ accompanied by time-varying weights (a consequence of the assumption that $\langle \vec{\omega}_i, \vec{\omega}_j \rangle = \langle \omega^2 \rangle_{ij} = \sigma_i^2 I_2 \delta_{ij}$ is diagonal). This assumption is reasonable if the evolution of frequencies $f_i$ of the basis $\Phi_i(t)$ is negligible during the considered time-span of observation $T$ (i.e. stationary phase). Such drifts in $f_i$ (typically chirp-like signals) would produce correlations between consecutive frequencies that are not accounted for in the expression of this kernel. Nevertheless, modeling the time dependent variance of each weight with this expression still captures information about non-stationary behaviors through the function $g(f, t)$.

The previous expression is valid for slowly evolving spectra that can be considered locally stationary. The notion of local stationarity is defined in [14] as when the signal is approximately stationary over some sliding interval $\tau$. In our case, because we want to model slowly evolving spectra, the stationarity timescale $\tau$ should be comparable to the total time of observation $T$. This condition translates to the function $g(f, t)$ as having most of its Fourier spectrum concentrated at frequencies around $1/T$. Essentially, we approximate non-stationary kernels by concentrating the 2D spectral measure $S(f, f')$ along the diagonal where $f = f'$ (see appendix V A).

In the upcoming section, we propose conducting simulations to assess the impact of this approximation on parameter recovery. The primary goal is to estimate whether neglecting correlations between Gaussian weights in Equation 2.6 significantly influences the inference of parameters within a Bayesian framework.

## C. Colored noise

In PTA, pulsar noise is typically described with a power-law spectrum $S_0(f) \propto Af^{-\gamma}$, where $A$ represents the amplitude and $\gamma$ denotes the spectral index [8, 22, 23]. Let us consider a time-dependent power spectral density with time-varying amplitude and spectral index:

$$ S(f, t) = \frac{A^2(t)}{12\pi^2} \left(\frac{f}{f_{yr}}\right)^{-\gamma(t)} f_{yr}^{-3}, \quad (2.7) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Example of a simulated non-stationary coloured noise for $log_{10} A(t_0) = -14$ and $\gamma(t_0) = 2$ observed in the timing residuals of a PTA pulsar with 15 years of observation. The noise was simulated using a first order Taylor expanded power-law (see Equation 4.1) with $a = 0.8$ and $b = -1.2$ meaning that we have a noise with increasing amplitude and decreasing spectral index (i.e. the noise whitens with time). The top panel shows the evolutionary PSD $S(f, t)$ and the bottom panel shows the corresponding time domain representation.</figcaption></figure>

with $f_{yr}$ the frequency corresponding to one year. Any arbitrary functions of time for $A(t)$ and $\gamma(t)$ can be used as long as we ensure $S(f, t) > 0$. We choose to express the amplitude and spectral index as:

$$ 
\begin{aligned}
\log_{10} A(t) &= \log_{10} A(t_0) + \sum_k a_k P_k(t - t_0), \\
\gamma(t) &= \gamma(t_0) + \sum_k b_k P_k(t - t_0),
\end{aligned}
\qquad (2.8) $$

where $t_0$ denotes the reference time, $a_k$ and $b_k$ are constant coefficients, and $P_k$ are polynomial functions (such as Chebyshev, spline, etc.). We can factor out the stationary PSD $S_0(f) \propto A^2(t_0)f^{-\gamma(t_0)}$ from a time-dependent function $g(f, t)$:

$$ 
\begin{aligned}
S(f, t) &= S(f, t_0) \times 10^{2 \sum_k a_k P_k(t - t_0)} \left(\frac{f}{f_{yr}}\right)^{-\sum_k b_k P_k(t - t_0)} \\
&= S_0(f) \times g^2(f, t),
\end{aligned}
\qquad (2.9) $$

providing an expression for the function $g(f, t)$ as presented in the previous section. According to Equation 2.6, we then have:

$$ 
\begin{cases}
C(t, t') = \sum_i \sigma_i^2 g(f_i, t) g(f_i, t') \cos(2\pi f_i (t - t')) \\
\sigma_i^2 \propto A^2(t_0) f_i^{-\gamma(t_0)} \Delta f_i \\
g(f, t) = 10^{\sum_k a_k P_k(t - t_0)} \left(\frac{f}{f_{yr}}\right)^{-\sum_k \frac{b_k}{2} P_k(t - t_0)}
\end{cases}
\qquad (2.10) $$

In the next sections of this paper, we will only utilize the first-order term of this expression. We will assume slowly varying properties of the signal.

### III. INJECTION RECOVERY TEST

#### A. Likelihood

To perform inference of the parameters, we need to define a likelihood. The GP model previously described assumes that the noise in the data is Gaussian. Hence, for a time series $s(\vec{t})$ of length $N$ measured at times $\vec{t}$, we will consider a Gaussian likelihood of the form:

$$ \mathcal{L}(s|\vec{\theta}) = \frac{1}{|2\pi C(\vec{\theta})|^{\frac{1}{2}}} \exp\left\{-\frac{1}{2} \sum_i \sum_j s(t_i) C^{-1}(\vec{\theta}, t_i, t_j) s(t_j)\right\} \qquad (3.1) $$

where $C^{-1}$ is the inverse of the covariance matrix $C$, $\vec{\theta}$ the vector of the model parameters and $|\cdot|$ denotes the determinant.

In the Bayesian framework, we define prior probability distributions $\pi(\vec{\theta})$ for the parameters $\vec{\theta}$ of the model. Combining the prior distributions with the likelihood, we can update our knowledge about $\vec{\theta}$, defining the posterior distribution $p(\vec{\theta}|s)$:

$$ p(\vec{\theta}|s) = \mathcal{L}(s|\vec{\theta})\pi(\vec{\theta})/\mathcal{Z}, \qquad (3.2) $$

with $\mathcal{Z}$ the so-called evidence of the model acting as a normalizing constant.

The key quantity to evaluate the significance of model A against model B in Bayesian analysis is the Bayes factor $\mathcal{B}_B^A$. It is defined as the ratio of the evidences

$$ \mathcal{B}_B^A = \frac{\mathcal{Z}_A}{\mathcal{Z}_B}. \qquad (3.3) $$

In this paper, the posterior distribution is sampled using Markov Chain Monte Carlo (MCMC) sampler PTMCMC [24] or nested sampling library dynasty [25].

#### B. Time-varying filter

In order to test our model, we want to simulate a colored noise with evolving properties in time using a
method that is independent of the one developed in the previous section. Various methods are available to generate such noise, for example, one can utilize auto-regressive models [12]. Here we chose to generate non-stationary data with a time-varying filter that allows one to control the amplitude and spectral index of the noise at each instant [14, 18].

We first generate an evenly sampled white noise series $y \sim N(0, 1)$ of size $N_s + 1$ with a cadence $\Delta t$ and a total duration $N_s \Delta t$. Then, the time series $y$ is filtered using a time-dependent filter $\hat{h}(f, t)$. The filtering is performed in the frequency domain. For a colored noise, we define a filter in Fourier domain $\hat{h}(f, t)$ at time $t$ as:

$$ \hat{h}(f, t) = \begin{cases} 0 & \text{for } f < 1/N_s \Delta t \\ g(f, t) H_0 \left( \frac{f}{f_{yr}} \right)^{-\alpha_0/2} & \text{for } f \ge 1/N_s \Delta t \end{cases} \quad (3.4) $$

where $H_0$ is the gain of the filter, and $\alpha_0$ is the slope. The gain and slope can be controlled over time using the function $g(f, t)$.

For a time series $y$ and its Fourier transform $\hat{y}(f)$, the output $y'(t)$ of the filter is:

$$ y'(t) = \mathcal{F}^{-1}(\hat{h}(f, t)\hat{y}(f)), \quad (3.5) $$

with $\mathcal{F}^{-1}$ the inverse Fourier transform.

Finally, we fit a quadratic polynomial of the form $a + bt + ct^2$ to the generated time series. The obtained polynomial is then subtracted from the time series. This final step mimics the timing model fit performed in PTA data [7, 8], which accommodates the quadratic spin-down (deceleration) of the pulsar rotation while also eliminating part of the moving average component attributed to the non-stationarity.

### C. Setting the model and priors

To test the model, we aim to estimate the Bayes factor for the non-stationary case versus the stationary case given the generated data. As mentioned earlier, we want to restrict ourselves to first-order variations in time of the power-law spectrum amplitude and spectral index. We evaluate the covariance matrix in Equation 2.6 on a discrete set of frequencies $f_i = i/\tau$ with $i = [1, ..., N_f]$ an integer, $N_f$ the size of the basis and $\tau$ the considered stationarity time-scale. In this context, the function $g(f_i, t)$ can be expressed as:

$$ \log_{10} g(f_i, t) = a \frac{(t - t_0)}{T} - \frac{1}{2} b \frac{(t - t_0)}{T} \log_{10}(f_i / f_{yr}). \quad (3.6) $$

where $T$ is the total time of observation, $f_{yr} = 1/(1yr)$ is a reference frequency and $t_0$ is the initial time.

In this work, we chose to fix $t_0$ at the start time of the entire dataset to directly measure the evolution of amplitude and spectral index from starting values $\log_{10} A(t_0)$ and $\gamma(t_0)$.

To perform inference of the parameters, we need to set prior probability distributions. For a slowly evolving signal observed within a time interval $T$, we assume that the stationarity time-scale $\tau \simeq T$ so the chosen function $g(f_i, t)$ must have most of its Fourier spectrum concentrated around the lowest frequency $1/T$. This ensures that the 2D spectral measure is concentrated on the diagonal $f = f'$, minimizing leakage between different frequencies (see section V A). If not, $g(f_i, t)$ is rapidly evolving, and the assumption of local stationarity is no longer valid.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Norm of the Fourier transform of $g(f_i, t)$ expressed in Equation 3.6 for $a = 1$, $b = 1$ and three ratios $f_i/f_{yr}$. The x-axis corresponds to the harmonics of $1/T$. We show $|\hat{g}(n)|/|\hat{g}(0)|$ for better visibility and comparison.</figcaption></figure>

In Figure 2, we display the Fourier spectrum of $g(f_i, t)$ on the interval $[t_0, t_0 + T]$, denoted as $\hat{g}(n)$ where $n$ corresponds to the harmonics of $1/T$. We can observe that for parameter values $a = b = 1$, the condition is roughly fulfilled. However, as we decrease the ratio $f_i/f_{yr}$ in the expression of $g(f_i, t)$, the condition is less respected. This implies that the Fourier transform $\hat{g}(n)$ of the function $g(f_i, t)$ has a wider spread at lower $f_i$, leading to more significant changes occurring within timescales shorter than $T$. Thus, we choose $a = b = 1$ to represent typical parameter values where the assumption of a slowly evolving spectrum begins to break down for $T \sim 10yr$. We decide to set a normal distributed prior $N(0, 1)$ for $a$ and $b$.

<table><caption>TABLE I. Prior probability distributions used for the parameters of a non stationary coloured noise.</caption><thead><tr><th>Parameter</th><th>Prior probability</th></tr></thead><tbody><tr><td>log<sub>10</sub> A(t<sub>0</sub>)</td><td>Uniform(-18, -10)</td></tr><tr><td>γ(t<sub>0</sub>)</td><td>Uniform(1, 6)</td></tr><tr><td>a</td><td>N(0, 1)</td></tr><tr><td>b</td><td>N(0, 1)</td></tr></tbody></table>
(see Table I) to confine their values around the interval [-1, 1]. The inference of parameters may not be precise enough to accurately estimate the evolutionary PSD of the process within the Bayesian scheme (resulting in wide posterior uncertainty). In this sense, this model remains an approximation. Nevertheless, we should be able to capture the non-stationarity of the process through the evaluation of the Bayes factor if a and b deviate too much from 0. This last point is crucial to test whether we are actually dealing with non-stationary noise.

### D. Simulation results

Following the procedure presented in section III B, we generate 1000 simulations with 500 samples and 20 years of data containing non-stationary red noise with a low white noise level of $10^{-9}$ s. For each simulation, the parameters a and b are drawn from their respective prior probability distribution, while the initial gain $H_0$ and slope $\alpha_0$ are fixed at $10^{-7}$ and 3, respectively. Inference is performed using the expressions of the likelihood and posterior probability presented in section III A. The prior distributions we use are given in Table I. We employ nested sampling to explore the parameter space for both stationary and non-stationary models to estimate the Bayes factors [25]. We denote the Bayes factor for non-stationarity against stationarity of noise as $\mathcal{B}_S^{NS}$.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Normalized histogram of the differences $\Delta\theta$ between injected and recovered maximum likelihood parameter value.</figcaption></figure>

Figure 3 illustrates the precision of non-stationary parameter recovery. We present a histogram showing the differences in each simulation between the injected and recovered maximum likelihood values of parameters a and b. The uncertainties are narrower than the prior distribution, indicating that the model is informative. However, the histogram for b appears broader than parameter a, which makes an accurate estimate of the evolution of the spectral index challenging.

Figure 4 demonstrates the ability to accurately predict non-stationarity through the evaluation of the Bayes fac-

<figure><img src="image_5.png" /><figcaption>FIG. 4. Normalized histogram of the Bayes factor $\mathcal{B}_S^{NS}$ non-stationary vs stationary noise obtained from each simulation. We define the level of non-stationarity $\sqrt{a^2 + b^2}$ combining parameters a and b. We plot the distribution obtained from their recovered maximum likelihood values for the corresponding simulation and Bayes factor in the x-axis.</figcaption></figure>

tor $\mathcal{B}_S^{NS}$. We present a histogram of the Bayes factors evaluated for each simulation. The distribution peaks at low Bayes factors, which corresponds to the fact that we draw a and b from $N(0, 1)$, hence they are often near zero. If a and b are small, the stationary model should be equally preferred. Additionally, we observe that the Bayes factor rapidly increases if a and b significantly deviate from zero. However, in this simulation, due to the large uncertainties on the recovery of b, the information it provides is limited, and the recovered maximum likelihood value does not effectively convey this lack of information (see section V B). In Figure 4, the non-stationarity level $\sqrt{a^2 + b^2}$ and the increase in the Bayes factor are predominantly influenced by the value of parameter a.

## IV. RESULTS

The results section will focus on the search for a GWB in the second data release of the EPTA collaboration [7]. The following analyses use the same methods, priors, and noise models presented in [6, 8]. We use the 24.7 years of dataset referred to as DR2full in [6]. The search is for a GWB with a power-law spectrum modeled as stochastic noise spatially correlated between pairs of pulsars, following the so-called Hellings-Downs correlation pattern [3]. The individual pulsar noise models are constructed from a combination of RN and DM noise, also modeled with power-law spectra [8]. These pulsar noise models are built using ENTERPRISE [26].

First, we study the GWB using the stationary formulation of the covariance matrix, as described in [6], but across different time slices of the dataset. Then, we investigate the GWB using the non-stationary formulation of the covariance that was developed in the previous sec-
tions

A. Time-slice analysis of GWB

Time-slice analyses of the GWB were previously con-
ducted by the Parkes Pulsar Timing Array collaboration
for their 3rd data release [5], where the dataset was seg-
mented into several slices and analyzed as a stationary
GWB in each slice. In this section, we perform a simi-
lar analysis on three different slices of our dataset. For
T, the total time span of observation, and t₀, the initial
time, we analyze the slices [t₀, T/2], [T/4, 3T/4], [T/2,
T], as illustrated in Figure 5.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Timeline representation (to scale) of the three ana-
lyzed slices represented by the colors blue, orange and green
with respect to the full time of observation T.</figcaption></figure>

In all three slices, we search for a stationary GWB with
a power-law spectrum on a 30-frequency linearly spaced
Fourier basis with the lowest frequency 2/T, allowing the
spectral index γ and amplitude log₁₀ A to vary (see [6]
for details). We denote the recovered GWB posteriors in
slices 1, 2, and 3 as GWB₁, GWB₂, and GWB₃ respectively.
The results are presented in Figure 6.

The GWB₁ posterior appears to be much broader than
the others. This is expected because the cadence of early
observations was lower than those of today (some pulsars
may have a cadence of 1 point per month in the early data
versus 1 every 3 days in the latest). The level of white
noise is directly influenced by the cadence of observa-
tion<fnref n="2" />. The interplay between the powerlaw spectrum of
red noise and the flat white noise tail may strongly influ-
ence the estimate of the red noise parameters. Moreover,
the quality of the backends processing raw data from the
radio observatories has improved over time. Their ability
to simultaneously observe incoming photons at multiple
frequency bands (larger bandwidth) allows for a better
characterization of the DM noise due to pulses traveling
through the interstellar medium [8]. The latter induces
a delay in the time of arrival of the pulses with an am-
plitude that is inversely proportional to the square of

the incoming photon's frequency. Its characterization re-
quires good multiband observing capabilities; otherwise,
it could be mistaken for RN and contaminate GWB pa-
rameter estimation. Improved cadence and noise charac-
terization should provide better constraints on the GWB
posteriors, which could explain why GWB₁ appears less
informative than GWB₂ and GWB₃. Moreover, if some of
the individual pulsar noises are non-stationary, they may
leak into the GWB signal. This last point will be ex-
plored in the next subsection. It is important to note
that this first slice, GWB₁, corresponds to the first data
release from EPTA, where no GWB was less significant.
The recovered posteriors are similar to what was observed
in [28].

The GWB₂ and GWB₃ posteriors nearly overlap and closely resemble the recovered posteriors using the latest 10 years of data in [6]. We observe a slight shift of the spectral index γ towards lower values. Figure 7 displays the 1D posteriors of γ for the three slices, indicating an estimated decrease in γ of approximately 0.66 from the median values. Regarding the amplitude, the median values show little to no evolution. However, evaluating the evolution of γ and log₁₀ A, especially for GWB₁, is challenging due to the large uncertainties.

One caveat of this analysis is that we lose information
about the low-frequency range after performing the slices.
Indeed, the lowest sampled frequency bin is at frequency
1/T, and increasing T raises its value. In PTA, the lowest
frequencies contain important information when search-
ing for a GWB. Additionally, this analysis prevents us
from estimating a Bayes Factor in favor of non-stationary
behavior against stationarity. The model presented in
this paper should address both issues.

B. Non-stationary powerlaw spectrum GWB

We will now analyze the data using a GWB model with
HD correlations and an evolutionary power-law spec-
trum. In section IIB, we developed an expression for
such a model.

The previous time-slice analysis only reveals small vari-
ations in the amplitude and spectral index of the power-
law. For this reason, we chose to limit ourselves to the
simplest possible model where the evolution of amplitude
and spectral index is approximated to first order in time.
We have:

$$
\begin{equation}
\begin{aligned}
\log_{10} A(t) &\sim \log_{10} A(t_0) + a(t - t_0)/T, \\
\gamma(t) &\sim \gamma(t_0) + b(t - t_0)/T,
\end{aligned}
\tag{4.1}
\end{equation}
$$

where $t_0$ is the minimum time of the dataset, $a$ the pa-
rameter controlling the evolution of log-amplitude and $b$
the parameter controlling the evolution in $\gamma$.

We performed the analysis using a normal prior $\mathcal{N}(0, 1)$ on $a$ and $b$ and show the recovered posteriors in Figure 8. We find a GWB with $\log_{10} A(t_0) = -14.4_{-0.2}^{+0.2}$,

<footnote n="2">²For an evenly sampled dataset with measurement uncertainties σ and red noise, the PSD can be expressed as S(f) = Af⁻¹ + 2σ²Δt where Δt is the cadence of observation that directly influences the tail of the spectrum [27].</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 6. 2D posterior distributions of GWB parameters analyzed in three time slices. The contours on the 2d histograms are the 68% and 95% credible regions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. 1D marginalized posterior distribution of GWB spectral index $\gamma$ analyzed in three time slices.</figcaption></figure>

$\gamma(t_0) = 3.7_{-0.5}^{+0.5}$, $a = 0.0_{-0.3}^{+0.3}$, and $b = -0.7_{-0.6}^{+0.7}$, fitting a spectrum with a constant amplitude and decreasing spectral index. This behavior aligns with the estimated evolution of the parameters using the medians of the 1D posteriors in the slice analysis. Still, the uncertainties on the estimate of parameter $b$ are very large and $b = 0$ lies within the $1 - \sigma$ credible interval. For this GWB model, we estimate a Bayes factor of non-stationarity against stationarity $\mathcal{B}_S^{NS} = 1.5$, indicating no evidence for a non-stationary GWB. The apparent evolution of the spectral properties of this process is likely due to improvements in receiver-backend systems and the enhanced cadence

of observation over the years. As demonstrated in the previous section, the inferred posteriors of the GWB parameters in the oldest slice of our dataset are very uninformative. This limitation may impede our ability to characterize the time-dependent features of the spectrum effectively.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Corner plot of non-stationary GWB parameters inferred from the EPTA dataset. The contours on the 2d histograms are the 68% and 95% credible regions. The dashed lines on the 1d histograms show the 14/50/86% percentiles.</figcaption></figure>

### C. Non-stationary individual noise

Earlier we mentioned the possible non-stationarity of individual pulsar noise. In the introduction, we presented the two dominating sources of time-correlated noise in pulsar timing: the RN, due to stochastic variation of the pulsar's spin rate, and the DM noise, due to the variation of the interstellar medium electron density throughout the travel of the pulses to Earth [8]. Another source of frequency dependent noise is identified, chromatic noise, with an amplitude proportional to $\nu^{-4}$ with $\nu$ the incoming photon frequency. This type of noise is assumed to be stationary in PTA data analysis. Incorrect modeling of pulsar noise could deteriorate the inference of the GWB parameters.

In this section, we compare the stationary formulation of the custom noise models presented in [8] to the non-stationary formula of section 4.2 in order to assess the stationarity of pulsar noises. Every noise spectrum is modeled as a power-law spectrum. We only test the first order model with function $g(f, t)$ given by
Equation 3.6 to allow each noise component to vary with time. We perform MCMC analysis using product-space sampling [29] to estimate the Bayes factor $\mathcal{B}_S^{NS}$. Results are gathered in Table II.

According to the distribution of Bayes factors, it appears that the pulsar noise do not exhibit first order variations in time, with the exception of J1713+0747 for which $\log_{10} \mathcal{B}_S^{NS} = 4.3$, hence favoring the non-stationary model.

Given the noise parameter values in Table III, it seems that the recovered non-stationarity originates from the DM noise with a relatively well constrained $a_{DM} - 0.7_{-0.2}^{+0.2}$ indicating that the DM amplitude decreases with time and $b_{DM} = 0.4_{-0.6}^{+0.6}$ indicating an increasing spectral index (though with large posterior uncertainties). For the RN, the $a_{RN}$ and $b_{RN}$ parameters are consistent with zero within the $1-\sigma$ credible interval suggesting stationary behavior.

However, J1713+0747 is known to have challenging noise properties [30–32] and non-stationary DM events were already reported for this pulsar in [33]. In [8], two exponential dip events were reported in the timing residuals of J1713+0747. An exponential dip is a deterministic DM-like effect that could be induced by extreme scattering events in the interstellar medium or a sudden change in the pulse profile of the pulsar, introducing transitory frequency dependent timing delays that can last up to a few months [34, 35]. In Table III, we see that the inclusion of the exponential dips in our model removes the apparent non-stationary behavior of J1713+0747, now showing $a_{DM} = -0.2_{-0.2}^{+0.3}$ and $b_{DM} = -0.4_{-0.7}^{+0.7}$. The Bayes factor for non-stationary noise versus stationary noise after the inclusion of exponential dip events is now $\log_{10} \mathcal{B}_S^{NS} = -0.6$.

Although not the best description, the non-stationary model was able to hint towards the presence of time-dependent features in the timing residuals of pulsar J1713+0747 through the evaluation of the Bayes factor. We tried searching again for a GWB in the EPTA dataset, now considering the non-stationarity (exponential dips) of J1713+0747, and found no difference in the estimate of the Bayes factor in favor of HD or in the inferred GWB parameters with respect to [6]. The stationary model of individual pulsar noise seems to be enough an approximation when searching for a GWB.

## D. Eccentric single source as non-stationary GWB

If a binary system of SMBH is eccentric, it will radiate polychromatic GWs, contrary to the monochromatic circular binaries [36]. On average, the GW signal can be decomposed as a sum of frequencies, which are the harmonics of the orbital frequency of the binary system. The waveform of an eccentric GW shows that the instantaneous frequency of the signal changes significantly within one orbital period [9, 10].

In [38], it was shown that single binaries could be mis-

<table><caption>TABLE II. Bayes factors for non-stationary against stationary noise models. The pulsar with $\log_{10} \mathcal{B}_S^{NS} > 2$ is shown in bold font. We observe that, in general, pulsars do not exhibit evidence of non-stationary noise, except for J1713+0747, for which we show the two values of $\log_{10} \mathcal{B}_S^{NS}$ obtained with (left) and without (right) including the exponential dip events in the model. The case of J1713+0747 is discussed in the text.</caption><thead><tr><th>Pulsar</th><th>log<sub>10</sub> B<sub>S</sub><sup>NS</sup></th></tr></thead><tbody><tr><td>J0030+0451</td><td>0.0</td></tr><tr><td>J0613-0200</td><td>-0.6</td></tr><tr><td>J0751+1807</td><td>-0.3</td></tr><tr><td>J0900-3144</td><td>-0.3</td></tr><tr><td>J1012+5307</td><td>-1.2</td></tr><tr><td>J1022+1001</td><td>-0.2</td></tr><tr><td>J1024-0719</td><td>-0.3</td></tr><tr><td>J1455-3330</td><td>0.0</td></tr><tr><td>J1600-3053</td><td>0.0</td></tr><tr><td>J1640+2224</td><td>-0.2</td></tr><tr><td><b>J1713+0747</b></td><td><b>4.3</b>/−0.6</td></tr><tr><td>J1730-2304</td><td>1.1</td></tr><tr><td>J1738+0333</td><td>0.1</td></tr><tr><td>J1744-1134</td><td>-0.5</td></tr><tr><td>J1751-2857</td><td>-0.2</td></tr><tr><td>J1801-1417</td><td>0.3</td></tr><tr><td>J1804-2717</td><td>0.2</td></tr><tr><td>J1843-1113</td><td>-0.4</td></tr><tr><td>J1857+0943</td><td>-0.4</td></tr><tr><td>J1909-3744</td><td>0.7</td></tr><tr><td>J1910+1256</td><td>0.0</td></tr><tr><td>J1911+1347</td><td>-0.4</td></tr><tr><td>J1918-0642</td><td>-0.3</td></tr><tr><td>J2124-3358</td><td>-0.3</td></tr><tr><td>J2322+2057</td><td>0.0</td></tr></tbody></table>

<table><caption>TABLE III. Table of median noise parameters (RN and DM) inferred from pulsar J1713+0747 timing data for the non-stationary model (NS), non-stationary model with exponential dips (NSed) and the stationary model with exponential dips (Sed). The displayed quantiles are 0.14 and 0.86.</caption><thead><tr><th>Parameter</th><th>Median (NS)</th><th>Median (NSed)</th><th>Median (Sed)</th></tr></thead><tbody><tr><td>γ<sub>DM</sub></td><td>1.0<sup>+0.4</sup><sub>-0.4</sub></td><td>1.9<sup>+0.6</sup><sub>-0.6</sub></td><td>1.6<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>log<sub>10</sub> A<sub>DM</sub></td><td>-13.0<sup>+0.1</sup><sub>-0.1</sub></td><td>-13.6<sup>+0.2</sup><sub>-0.2</sub></td><td>-13.7<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>a<sub>DM</sub></td><td>-0.7<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.2<sup>+0.3</sup><sub>-0.2</sub></td><td>-</td></tr><tr><td>b<sub>DM</sub></td><td>0.4<sup>+0.6</sup><sub>-0.6</sub></td><td>-0.4<sup>+0.7</sup><sub>-0.7</sub></td><td>-</td></tr><tr><td>γ<sub>RN</sub></td><td>3.2<sup>+0.8</sup><sub>-0.8</sub></td><td>3.2<sup>+0.8</sup><sub>-0.8</sub></td><td>3.5<sup>+0.6</sup><sub>-0.6</sub></td></tr><tr><td>log<sub>10</sub> A<sub>RN</sub></td><td>-13.9<sup>+0.4</sup><sub>-0.5</sub></td><td>-14.0<sup>+0.4</sup><sub>-0.5</sub></td><td>-14.3<sup>+0.3</sup><sub>-0.3</sub></td></tr><tr><td>a<sub>RN</sub></td><td>-0.3<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.2<sup>+0.5</sup><sub>-0.5</sub></td><td>-</td></tr><tr><td>b<sub>RN</sub></td><td>-0.2<sup>+0.9</sup><sub>-0.9</sub></td><td>-0.3<sup>+0.9</sup><sub>-0.9</sub></td><td>-</td></tr></tbody></table>

taken for a GWB with a power-law spectrum. For eccentric binaries with orbital periods comparable to the total time of observation of a PTA, the instantaneous change in frequency of the GW signal can be seen as a GWB with evolving spectral index $\gamma(t)$ (see Figure 9). Indeed, the $\gamma$ controls the balance between low and high frequencies
<figure><img src="image_10.png" /><figcaption>FIG. 9. One period of an eccentric signal (Earth term only) with eccentricity 0.6 and orbital period of 15 years We see that the instantaneous frequency significantly changes within one period. The waveform was generated using **ENTERPRISE** extensions [37].</figcaption></figure>

of the power-law spectrum and the presence of an eccentric signal would introduce time-dependent features. GW sources of higher eccentricity and higher amplitude would introduce more non-stationarity because of their stronger frequency evolution or higher signal to noise ratio (SNR). In this section, we want to investigate what we would see when a stationary GWB is influenced by an eccentric single source with low SNR.

We create a simulation with 10 pulsars (isotropically distributed in the sky), 15 years of observation and only white noise at a level of $10^{-7}$s. The simulation contains a stationary GWB with $\log_{10} A = -14.2$ and $\gamma = 3$. We generate two cases: (i) only a stationary GWB is present in the data (ii) the same stationary GWB is accompanied by an eccentric source located at the Virgo galaxy cluster with amplitude $h = 10^{-13.5}$ and orbital period of 15 years (here corresponding to SNR $\simeq 1.8$). We produce 8 datasets with different eccentricities of the source $e$ ranging from 0.1 to 0.8. We analyzed these datasets with the non stationary model presented in IIB, considering HD correlations. In Figure 10, we show that the presence of an eccentric source with $e = 0.5$ introduces non stationarity in the inferred GWB parameters. For case (i), parameters $a$ and $b$ are near zero meaning that we are dealing with a stationary signal. In case (ii), where the eccentric signal is introduced, they deviate significantly from zero. Note that in our simulation, the parameter $b$ is negative (i.e. $\gamma$ is decreasing) showing that some specific realizations of eccentric binaries can somewhat reproduce the behavior observed in the real data.

Figure 11

The evolution of the recovered Bayes factor $\mathcal{B}_S^{NS}$ and parameters $a$ and $b$, as a function of the eccentricity $e$ of the injected GW signal, is illustrated in Figure 11. It is evident that the model effectively captures the non-stationary features induced by the presence of the eccentric

<figure><img src="image_11.png" /><figcaption>FIG. 10. Corner plot of non-stationary GWB parameters. In orange, we show the results for a simulated stationary GWB with $\log_{10} A = -14.2$ and $\gamma = 3$. In blue, we show the results for the same realisation of GWB "contaminated" by an eccentric single source with 15 years orbital period and eccentricity $e_0 = 0.5$. The dashed line are the medians. The displayed values of parameters correspond to the medians of the orange histograms. The contours on the 2d histograms are the 68% and 95% credible regions.</figcaption></figure>

source when the latter is loud enough. The Bayes factor begins to favor the non-stationary model for eccentricities higher than 0.4. For this analysis, we have utilized the expression of $g(f, t)$ from Equation 3.6. In future work, it might be beneficial to design a $g(f, t)$ that better suits this problem.

V. DISCUSSION AND CONCLUSION

The derived kernel serves the purpose of modeling non-stationarity in noise. It resembles the stationary kernels commonly used in PTA analyses but incorporates the capability to capture the time-dependent variance of the Gaussian weights. In its current form, this kernel remains an approximation and is valid only for small deviations from a stationary model. A significant challenge lies in understanding how correlations between frequencies can be accurately modeled, leading to an exact and non-approximated version of the non-stationary kernel. Further exploration in this area is necessary to improve our understanding and refine the methodology.

This approach necessitates a prior assumption about the expression of the evolutionary PSD. One significant advantage is its consistent expression of the covariance
<figure><img src="image_12.png" /><figcaption>FIG. 11. (top) Bayes factor for non-stationary against stationary GWB models as a function of the eccentricity of the injected source, (bottom) recovered posterior median of non-stationary parameters *a* and *b* as a function of eccentricity. The shaded region around *a* and *b* represents the 1 − σ credible interval.</figcaption></figure>

matrix for the entire dataset, unlike the time-slice analy-
sis, which necessitates cutting the data into several seg-
ments. This can be particularly valuable for applications
like PTA, where noise is most pronounced in the lowest
frequency bins.

In section IIID, we showed that our model can oc-
casionally produce broad posterior uncertainties for the
non-stationary parameters, and the stationary solution
may fall within the 1 − σ credible interval. Consequently,
accurately characterizing the evolutionary PSD can be
challenging. Thus, it becomes crucial to assess the Bayes
factor against the stationary model to determine whether
the preference lies with a genuinely non-stationary signal.

The slice analysis reveals a clear evolution in the in-
ferred GWB properties. The first epoch appears much
less informative than the following two. This can be at-
tributed to an evolution of the data quality over time or
the leakage of non-stationary features of individual pulsar
noises into the GWB estimate. We conducted individual
analyses for each pulsar to investigate non-stationary fea-
tures and found that only J1713+0747 seemed to exhibit
time-dependent DM noise amplitude variations. Non-
stationary DM events were already reported for this pul-
sar in previous studies. Exponential dip events offer the
best description of these events.

The non-stationary model developed indicates that
there is no evidence for a non-stationary GWB in the
second data release of the EPTA data. Based on the re-
sults of the slice analysis, we infer that the apparent evo-

lution of the inferred parameter values may stem from
a poor characterization of individual pulsar noise con-
taminating the observed GWB. However, there is still
a possibility that the presence of GW signals induced
by eccentric SMBHBs in the data could produce simi-
lar non-stationary behaviors. The last section explored
this hypothesis and showed that PTAs could detect non-
stationary features of the GWB if an eccentric source is
present in the data.

With the current dataset, it remains inconclusive
whether the observed evolution solely results from a
change in backend quality or represents a genuine non-
stationary feature of the GWB. Furthermore, if the noise
from individual pulsars exhibits non-stationarity, it could
potentially contaminate the GWB estimate, making it
appear as non-stationary. Further investigation into this
matter is warranted in future studies. Additionally, the
upcoming release of the International Pulsar Timing Ar-
ray (IPTA) third data release, with enhanced observation
cadence and improved frequency coverage, is expected to
provide better noise characterization and definitive in-
sights into this issue.

ACKNOWLEDGEMENTS

I thank Stas Babak, Philippa Cole, Aurélien
Chalumeau and Hippolyte Quelquejay for the feedback
on the manuscript. I also thank Philippe Bacon and
Quentin Baghi for the interesting discussions.

M.F. acknowledges financial support from MUR PRIN,
Grant No. 2022-MYL2X, project 'GRAPE', funded by
the European Union - Next Generation EU. The Nançay
radio Observatory is operated by the Paris Observa-
tory, associated with the French Centre National de
la Recherche Scientifique (CNRS), and partially sup-
ported by the Region Centre in France. IC, LG, GT
acknowledge financial support from "Programme Na-
tional de Cosmologie and Galaxies" (PNCG), and "Pro-
gramme National Hautes Energies" (PNHE) funded by
CNRS/INSU-IN2P3-INP, CEA and CNES, France. IC,
LG, GT acknowledge financial support from Agence Na-
tionale de la Recherche (ANR-18-CE31-0015), France.
The Sardinia Radio Telescope (SRT) is funded by the De-
partment of University and Research (MIUR), the Ital-
ian Space Agency (ASI), and the Autonomous Region of
Sardinia (RAS) and is operated as a National Facility by
the National Institute for Astrophysics (INAF). DP ac-
knowledges support from the INAF Large Grant "Grav-
itational Wave Detection using Pulsar Timing Arrays"
(P.I. D. Perrodin, 2023). JA acknowledges support from
the European Commission (ARGOS-CDS; Grant Agree-
ment number: 101094354). GMS sacknowledge finan-
cial support provided under the European Union's H2020
ERC Consolidator Grant "Binary Massive Black Hole As-
trophysics" (B Massive, Grant Agreement: 818691). The
Westerbork Synthesis Radio Telescope is operated by the
Netherlands Institute for Radio Astronomy (ASTRON)
with support from the Netherlands Foundation for Scientific Research (NWO).

# APPENDIX

## A. 2D spectral measure

The generalized Wiener-Khintchine relations are defined in [18] as

$$ 
\begin{aligned}
C(t, t') &= \int df df' S(f, f') e^{2\pi i (tf - t'f')} \\
S(f, f') &= \int dt dt' C(t, t') e^{-2\pi i (tf - t'f')}
\end{aligned}
\quad (5.1) $$

From the definition of the non-stationary covariance matrix in Equation 2.6, we can estimate the 2D spectral density:

$$ S(f, f') = \int dt dt' \sum_i \sigma_i^2 g_i(t) g_i(t') \left\{ e^{-2\pi i t(f+f_i)} e^{-2\pi i t'(-f-f_i)} + e^{-2\pi i t(f_i-f)} e^{-2\pi i t'(-f'-f_i)} \right\} \quad (5.2) $$

with $\hat{g}_i(f)$ the Fourier transform of the function $g_i(t)$, we have

$$ S(f, f') = \sum_i \frac{\sigma_i^2}{2} \left[ \hat{g}_i(f+f_i)\hat{g}_i(-f'-f_i) + \hat{g}_i(-f+f_i)\hat{g}_i(f'-f_i) \right] \quad (5.3) $$

Note that if $\hat{g}_i(f) = \delta_i(f)$ (a Dirac delta) and applying Equation 5.1, we recover the stationary covariance matrix of Equation 2.3, giving a 2D spectral measure $S(f, f')$ that is only defined on the diagonal $f = f'$ and evaluated at discrete frequencies $f_i$. The spread of the function $\hat{g}_i(f)$ determines how much we deviate from the stationary covariance matrix.

## B. Simulations

In section III B, we presented a method to generate a non-stationary colored noise. Here, we show an example of one realisation of non-stationary noise using this method. In Figure 13, we show the corner plots with the injected versus recovered values of parameter. We can see that the model performs well, except for parameter $b$, where the posterior uncertainties are quite large.

## C. PP-plots

In [39], a validation method is presented. It necessitates both data-generation and model-fitting software.

<figure><img src="image_13.png" /><figcaption>FIG. 12. Time-domain plot of a non-stationary noise generated with the method presented in III B.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 13. Corner plot of non-stationary parameters inferred from the simulation in Figure 12. The solid black lines show the injected values of parameters. The contours on the 2d histograms are the 68% and 95% credible regions. The dashed lines on the 1d histograms show the 14/50/86% percentiles.</figcaption></figure>

In this context, our model-fitting software is PTMCMC [24], which has been previously utilized. During the fitting process, we sample the posterior distribution employing our model, specifically the non-stationary kernel, using the fitting software. In each simulation, we possess knowledge of the injected parameter values $\tilde{\theta}_0$, allowing us to calculate the corresponding quantiles $q(\tilde{\theta}_0)$ with respect to the posterior distribution derived from the fitting software. The methodology in [39] demonstrates that if
the fitting software operates accurately and the gener-
ated data adheres to the same posterior distribution for
parameters as the constructed model, then the quantiles
$q(\tilde{\theta}_0)$ must be distributed as Uniform(0, 1).

In Figure 14 we show the PP plot for the quantiles
$q(\tilde{\theta}_0)$ obtained from simulations. They seem to follow the
relationship $y=x$, but we notice discrepancies, especially
for parameter $b$ which is the one with larger uncertainty.

<figure><img src="image_15.png" /><figcaption>FIG. 14. PP plot comparing the fraction of simulations within confidence interval (CI) and CI. The grey shaded region represent the $1 - \sigma$ error from $y=x$.</figcaption></figure>

[1] A. Sesana, A. Vecchio, and C. N. Colacino, Monthly Notices of the Royal Astronomical Society **390**, 192–209 (2008).

[2] E. S. Phinney, “A practical theorem on gravitational wave backgrounds,” (2001), arXiv:astro-ph/0108028 [astro-ph].

[3] B. Allen and J. D. Romano, Physical Review D **108** (2023), 10.1103/physrevd.108.043026.

[4] G. Agazie et al., The Astrophysical Journal Letters **951**, L8 (2023).

[5] D. J. Reardon et al., The Astrophysical Journal Letters **951**, L6 (2023).

[6] J. Antoniadis et al. (EPTA, InPTA), Astron. Astrophys. **678**, A50 (2023), arXiv:2306.16214 [astro-ph.HE].

[7] J. Antoniadis et al. (EPTA, InPTA), Astron. Astrophys. **678**, A48 (2023).

[8] J. Antoniadis et al. (EPTA, InPTA), Astron. Astrophys. **678**, A49 (2023), arXiv:2306.16225 [astro-ph.HE].

[9] S. R. Taylor, E. A. Huerta, J. R. Gair, and S. T. McWilliams, The Astrophysical Journal **817**, 70 (2016).

[10] A. Susobhanan, A. Gopakumar, G. Hobbs, and S. R. Taylor, Physical Review D **101** (2020), 10.1103/phys-revd.101.043022.

[11] J. Ellis and N. Cornish, Physical Review D **93** (2016), 10.1103/physrevd.93.084048.

[12] N. J. Cornish, Phys. Rev. D **102**, 124038 (2020).

[13] R. van Haasteren and M. Vallisneri, Phys. Rev. D **90**, 104012 (2014).

[14] S. Mallat, G. Papanicolaou, and Z. Zhang, The Annals of Statistics **26**, 1 (1998).

[15] A. Papoulis and S. U. Pillai, Probability, Random Variables, and Stochastic Processes, 4th ed. (McGraw Hill, Boston, 2002).

[16] C. E. Rasmussen, “Gaussian processes in machine learning,” in Advanced Lectures on Machine Learning, edited by O. Bousquet, U. von Luxburg, and G. Rätsch (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 63–71.

[17] R. van Haasteren and M. Vallisneri, Monthly Notices of the Royal Astronomical Society **446**, 1170 (2014), https://academic.oup.com/mnras/article-pdf/446/2/1170/13767575/stu2157.pdf.

[18] M. B. Priestley, Journal of the Royal Statistical Society. Series B (Methodological) **27**, 204 (1965).

[19] H. Hong, Mechanical Systems and Signal Processing **190**, 110131 (2023).

[20] Y.-L. K. Samo and S. Roberts, “Generalized spectral kernels,” (2015), arXiv:1506.02236 [stat.ML].

[21] S. Remes, M. Heinonen, and S. Kaski, in Proceedings of the 31st International Conference on Neural Information Processing Systems (NIPS’17) (Curran Associates Inc., Red Hook, NY, USA, 2017) p. 4645–4654.

[22] G. Agazie and other, The Astrophysical Journal Letters **951**, L10 (2023).

[23] D. J. Reardon et al., The Astrophysical Journal Letters **951**, L7 (2023).

[24] J. Ellis and R. van Haasteren, “jellis18/ptmcmcsampler: Official release,” (2017).

[25] J. Skilling, Bayesian Analysis **1**, 833 (2006).

[26] J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, “Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite,” Zenodo (2020).

[27] C. J. Moore, S. R. Taylor, and J. R. Gair, Classical and Quantum Gravity **32**, 055004 (2015).

[28] Lentati, L. et al., Monthly Notices of the Royal Astronomical Society **453**, 2576 (2017).
(2015), https://academic.oup.com/mnras/article-pdf/453/3/2576/4900884/stv1538.pdf.

[29] S. Hee, W. J. Handley, M. P. Hobson, and A. N. Lasenby, Monthly Notices of the Royal Astronomical Society **455**, 2461 (2015), https://academic.oup.com/mnras/article-pdf/455/3/2461/9377568/stv2217.pdf.

[30] Agazie, Gabriella et al., The Astrophysical Journal Letters **951**, L10 (2023).

[31] M. T. Lam, Research Notes of the AAS **5**, 167 (2021).

[32] J. Singha et al., Monthly Notices of the Royal Astronomical Society: Letters **507**, L57 (2021), https://academic.oup.com/mnrasl/article-pdf/507/1/L57/54645627/slab098.pdf.

[33] L. Lentati et al., Monthly Notices of the

Royal Astronomical Society **458**, 2161 (2016), https://academic.oup.com/mnras/article-pdf/458/2/2161/18241404/stw395.pdf.

[34] R. M. Shannon et al., The Astrophysical Journal Letters **828**, L1 (2016).

[35] W. A. Coles et al., The Astrophysical Journal **808**, 113 (2015).

[36] P. C. Peters and J. Mathews, Phys. Rev. **131**, 435 (1963).

[37] S. R. Taylor, P. T. Baker, J. S. Hazboun, J. Simon, and S. J. Vigeland, "enterprise extensions," (2021), v2.4.3.

[38] J. Antoniadis et al. (EPTA, InPTA), (2023), arXiv:2306.16226 [astro-ph.HE].

[39] S. R. Cook, A. Gelman, and D. B. Rubin, Journal of Computational and Graphical Statistics **15**, 675 (2006).