A Gaussian-processes approach to fitting for time-variable spherical solar wind in pulsar timing data

Propagation effects are one of the main sources of noise in high-precision pulsar timing. For pulsars below an ecliptic latitude of $5^\circ$, the ionised plasma in the solar wind can introduce dispersive delays of order 100 microseconds around solar conjunction at an observing frequency of 300 MHz. A common approach to mitigate this assumes a spherical solar wind with a time-constant amplitude. However, this has been shown to be insufficient to describe the solar wind. We present a linear, Gaussian-process piecewise Bayesian approach to fit a spherical solar wind of time-variable amplitude, which has been implemented in the pulsar software run_enterprise. Through simulations, we find that the current EPTA+InPTA data combination is not sensitive to such variations; however, solar wind variations will become important in the near future with the addition of new InPTA data and data collected with the low-frequency LOFAR telescope. We also compare our results for different high-precision timing datasets (EPTA+InPTA, PPTA, and LOFAR) of three millisecond pulsars (J0030$+$0451, J1022$+$1001, J2145$-$0450), and find that the solar-wind amplitudes are generally consistent for any individual pulsar, but they can vary from pulsar to pulsar. Finally, we compare our results with those of an independent method on the same LOFAR data of the three millisecond pulsars. We find that differences between the results of the two methods can be mainly attributed to the modelling of dispersion variations in the interstellar medium, rather than the solar wind modelling.


INTRODUCTION
Pulsar timing consists in recording the times of arrival (ToAs) of highly stable pulses emitted by pulsars, which are then compared with predictions from long term models.These models take into account the pulsar's behaviour, as well as factors such as astrometric effects (e.g.spin frequency and its derivatives, position, etc.), binary companions, or the dispersive delays induced by the ionised medium through which the pulsar signal propagates (for more details see e.g.Edwards et al. 2006).The left-over signal after subtracting the model from the observed ToAs is often referred to as 'timing residuals', or just 'residuals', and is expected to be only white noise if the model is optimal.
Due to the high stability of pulsar rotation, and particularly of the population of recycled millisecond pulsars (MSPs; Backer et al. 1982), high-precision pulsar timing constitutes a great tool for a large variety of scientific investigations (see e.g.Manchester 2017).Huge efforts are currently being concentrated on gravitational wave searches using decades of observations of a large sample of MSPs with multiple telescopes, referred to as Pulsar Timing Arrays (PTAs; e.g.Tiburzi 2018).PTA experiments are expected to be primarily sensitive to the gravitational wave background (GWB) in the nanohertzfrequency regime, most likely originating from supermassive black hole binary mergers (e.g.Burke-Spolaor et al. 2019).Three major collaborations have historically been involved in the search for the GWB, namely the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the Parkes Pulsar Timing Array (PPTA; Manchester et al. 2013), and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav;Demorest et al. 2013); these are also the founding members of the International Pulsar Timing Array (IPTA; Verbiest et al. 2016) collaboration, having recently been joined by the Indian Pulsar Timing Array (InPTA; Joshi et al. 2022).Recently, the three aforementioned collaborations (with the InPTA working along the EPTA) all coherently reported marginal evidence for a GWB signal (EPTA Collaboration et al. 2023c;Reardon et al. 2023;Agazie et al. 2023).While neither of these currently meet the requirements for being defined as a clear detection, further investigations into the analysis methods, as well as the upcoming combined IPTA dataset, are expected to improve on these GWB measurements.
The signal induced by the GWB in pulsar data is expected to be extremely weak, even when correlated over tens of pulsars (e.g.Siemens et al. 2013;Janssen et al. 2015).Therefore, other effects in pulsar timing data must be carefully considered, as they can obscure, or even mimic, a GWB signal (Tiburzi et al. 2016).One of the strongest sources of 'noise' in this context are the dispersive delays introduced in the ToAs by the interaction between the radio waves and the ionised medium through which the pulsar signal propagates on its way to the observer (Lentati et al. 2016).These dispersive delays are modelled as having an inverse-squared dependency with the observing frequency,  obs .The dispersive delay on a ToA is expressed as where  D ≃ 2.41×10 −4 MHz −2 pc cm −3 s −1 is a dispersion constant (Manchester & Taylor 1972), and the dispersion measure DM is defined as the free electron number density,  e , integrated over the line of sight to the pulsar, ℓ; the DM is usually quoted in pc cm −3 .As illustrated by Eq. 1, the dispersion delay is stronger at lower observing frequencies.The noise due to the dispersion delay is dominated by the effects of the turbulent and inhomogeneous ionised interstellar medium (IISM) along the line-of-sight, which can induce fluctuations in the DM of order 10 −3 pc cm −3 over a timescale of years (e.g.Keith et al. 2013;Jones et al. 2017;Donner et al. 2020).Several possible mitigating strategies for the turbulent IISM contribution are used throughout the PTAs, such as modelling it as a chromatic rednoise Gaussian process in Fourier space (e.g.Lentati et al. 2014), or using a time-domain piecewise binned model (e.g. the DMX model; NANOGrav Collaboration et al. 2015).
For pulsars with a line-of-sight that passes close to the Sun, the delay induced by the propagation through the solar wind (SW hereafter) is also noticeable in the current quality of pulsar data, inducing DM fluctuations as high as 10 −3 -10 −4 pc cm −3 .Furthermore, Tiburzi et al. (2016) showed that, if not carefully considered, the influence of the SW may mimic a GWB signal in PTA-like data, as the SW can create spatial correlations among the pulsars.To account for the influence of the SW in pulsar timing data, PTA collaborations generally use a simple 'spherical SW' model, based on the assumption that the number density of ionised electrons varies under a spherically symmetric law away from the Sun, according to the inverse square law (Edwards et al. 2006), i.e.
where  is the position vector from the Sun to the point of interest affected by the SW, generally given in astronomical units (au);  SW e is the amplitude of the number density at 1 au, and we generally refer to it as the 'amplitude of the SW' in this work.Under the spherical SW assumption, this amplitude  SW e is space invariant (does not depend on the vector ); thus integrating the number density as per Eq. 2 over the line-of-sight gives the DM contribution as (e.g.You et al. 2007) DM SW  sph ≃ 4.85 × 10 −6  SW e cm −3 where  is the solar elongation angle of the pulsar, i.e. the pulsarobserver-Sun angle, which is minimum at the solar conjunction of the pulsar.Replacing the DM given as per Eq. 3 into Eq. 1 gives the time delay of a ToA, due to a spherical SW, which we factorise as where we have summarised some of the physical constants and geometrical dependence in the variable In the standard approach to mitigating the effect of the SW, the amplitude  SW e in the described spherically symmetric model is constant in time; in e.g. the recent EPTA dataset (EPTA Collaboration et al. 2023a), this amplitude is generally kept fixed at 7.9 cm −3 , as per Madison et al. (2019).Furthermore, data taken when the pulsar appears < 5 • away from the Sun in the sky are considered to be poorly described by this simple model and commonly removed (Verbiest et al. 2016).
No model has yet been developed that fully captures the observed impact of the SW, especially on low-frequency pulsar data, as shown by e.g.Tiburzi et al. (2019).You et al. (2007) proposed a model based on the coronal magnetograms derived by the Wilcox Solar Observatory, and the bi-modal nature of the SW, considering the contributions from both a slow (equatorial) and a fast (polar) solar stream (e.g.Coles 1996).You et al. (2007) argued that this two-phase model performed better than the spherical SW for their PPTA observations.However, Tiburzi et al. (2019) found that the spherically symmetric SW model with time-dependent amplitude performed systematically better than the two-phase model in removing the SW contribution in their longer-span and lower-frequency LOFAR data.The apparent discrepancy between these two analyses is thought to be due to either (i) the increased DM precision of the lower observing frequency of the data used in Tiburzi et al. (2019) as compared with that of You et al. (2007); or (ii) a difference in performance of the two-phase model with the heliospheric latitude of the pulsar, since the two papers investigated data from different pulsars.
A clear improvement to current general models is to allow the amplitude of the SW to vary each year (i.e.be time-variable) in the spherically symmetric model.Tiburzi et al. (2021), hereafter referred to as T21, have already shown this to be beneficial and more adequate for low-frequency data taken with the European interferometer LOFAR (Van Haarlem et al. 2013).T21 fit a time-variable spherical SW to each pulsar, and a clear temporal variation was observed in several of the analysed pulsars.Note that previously, Madison et al. (2019) had reported little evidence for long-term variations in the SW density using the NANOGrav 11-yr dataset.
In this work, we present a Gaussian-process piecewise approach, implemented as part of the pulsar analysis toolkit run_enterprise (Keith et al. 2022), which allows for an automatic time-variable spherical SW fitting in a Bayesian framework, simultaneously with all the other pulsar timing parameters and noise models.Recently, Hazboun et al. (2022) developed several comprehensive Bayesian algorithms that allow fitting for SW across several pulsars simultaneously.They used a uniform-prior piecewise model that globally fits for a spherical SW amplitude in each temporal bin (of e.g. 3 months) of the data, also exploring variations in the exponent of the 1/ 2 law of Eq. 2; further, they also explore a model based on globally fitting for continuous Fourier-basis, time-dependent variations in the SW.In this work, we present a different approach to SW Bayesian fitting.We use a simple piecewise algorithm based on a Gaussian process with a Normal-distribution prior to fit for a time-variable amplitude spherical SW in each pulsar individually.The choice of prior allows us to marginalise over the individual yearly amplitudes, reducing the problem to a single additional hyperparameter in the width of the Gaussian prior.This, together with the mathematical implementation that keeps the parameter estimation algebraically linear, makes any additional computational cost negligible.
This article is structured as follows: in Section 2, the main properties of the datasets are summarised; in Section 3 we describe the algorithm implemented in our pipeline and simulations; in Section 4 the results of our analysis are discussed; and in Section 5 we summarise our conclusions.

DATASETS
For this work, we selected three millisecond pulsars which pass in close proximity (≤ 5.31 • ) to the Sun during conjunction, and are included in the study by T21.These are PSRs J0030+0451, J1022+1001, and J2145−0450, and they are present both in LO-FAR and PTA observations.Here we consider data from the recent EPTA+InPTA data combination (EPTA Collaboration et al. 2023a;Tarafdar et al. 2022), as well as from the Second Data Release (DR2) of PPTA (Reardon et al. 2021).Table 1 shows the ecliptic latitudes of the three pulsars of interest, equivalent to the sky-angle between the pulsar and the Sun at conjunction.It also shows which data were available for each pulsar.The datasets, each with different frequency coverage and properties (discussed in more detail in Sections 2.1-2.3), are used separately to construct simulations, and thus assess their sensitivity to the SW.Further, the results of our pipeline are compared across all the datasets, as well as with the independent measurements in T21.

EPTA (DR2) + InPTA (DR1) combined dataset
The Second Data Release (DR2) of the EPTA collaboration is used in this work, augmented with the First Data Release (DR1) of the InPTA (EPTA Collaboration et al. 2023a;Tarafdar et al. 2022).Here, we refer to this combined dataset as EPTA+InPTA.
There are five European radio telescopes which provided data for the EPTA DR2 dataset, namely: the 100-m Effelsberg Telescope (in Germany), Jodrell Bank Observatory's 76-m Lovell Telescope (in the United Kingdom), Nançay Radio Observatory's large Radio Telescope (NRT; in France), the Astronomical Observatory of Cagliari's 64-m Sardinia Radio Telescope (SRT; in Italy), and the Westerbork Synthesis Radio Telescope (WSRT; in the Netherlands).Moreover, these telescopes were also used collectively, on a monthly cadence, as the Large European Array for Pulsars (LEAP), equivalent to a 194-m sixth interferometric telescope in the EPTA (Bassa et al. 2016).The EPTA DR2 dataset contains observations of 25 millisecond pulsars, up to 25 years in length.The large majority of these observations were taken at frequencies in the 'L-band' (1-2 GHz) and above, with bandwidths of up to 512 MHz, while there were a limited number of observations centred at lower frequencies of 350 MHz.For a detailed description of the properties of the EPTA DR2 telescopes and observations, see Chen et al. (2021), andEPTA Collaboration et al. (2023a).
The InPTA dataset includes observations taken with the upgraded Giant Meterwave Radio Telescope (uGMRT; in India) over a period of 3.5 years (for a detailed description see Tarafdar et al. 2022).The uGMRT took simultaneous observations at two frequency bands, referred to as 'B3' (300-500 MHz), and 'B5' (1260-1460 MHz), respectively.These simultaneous observations at different frequencies, as well as the wide bandwidths, make the InPTA observations extremely valuable for measuring the DM influence in pulsar data, including the effect of the SW.

LOFAR dataset
For PSRs J0030+0451 and J2145−0750, the same LOFAR datasets as in T21 are used; while for PSR J1022+1001, the T21 dataset was supplemented with the most recent ∼ 2 years of observations for this work.These data were taken with subsets of the International LOFAR telescope (Van Haarlem et al. 2013;Stappers et al. 2011), namely six of the German stations, the Swedish station, and the LOFAR Core in the Netherlands.More than 100 pulsars are regularly observed using this setup.The extremely low observing frequency of the LOFAR instruments, covering a range of roughly 110-190 MHz, makes this telescope incredibly well-suited for studying chromatic effects on pulsar data such as the SW.For more details on the LOFAR datasets, see e.g.Porayko et al. (2019), Donner et al. (2019), Tiburzi et al. (2019).

PPTA (DR2) dataset
We also make use of the open-access PPTA DR2 timing data (Reardon et al. 2021) to validate and compare with the EPTA+InPTA results.This PPTA dataset is described in detail in Kerr et al. (2020) and was taken with the Australian 64-m Murriyang Parkes Radio Telescope.It spans ∼ 14 years (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) of observations of 26 millisecond pulsars at three frequency bands, roughly centred at 700, 1400 and 3100 MHz.The observing cadence for each pulsar was approximately two-three weeks.The ToAs are available in both sub-banded and band-averaged form, of which we used the latter in this work, for simplicity.In particular, the PPTA DR2 data for PSRs J1022+1001 and J2145−0450 were used to compare SW results with those of the LOFAR and EPTA+InPTA ToAs.The data of PSR J2145−0450 was also supplemented with archival observations of the Parkes telescope since 1994, as published in Manchester et al. (2013); these were however much less sensitive to the SW.

Timing model
In this work, we use open-source pulsar software, namely the Bayesian fitting software enterprise (Ellis et al. 2019), which has been integrated into the pulsar analysis toolkit run_enterprise (Keith et al. 2022), together with the pulsar timing software tempo2 (Edwards et al. 2006).
The Bayesian method allows for simultaneous fitting of deterministic pulsar parameters (such as the spin frequency and its derivatives, position, known binary companions, etc.), which are generally marginalised over, unless of particular interest; and of parameters characterising stochastic timing noise.In this framework, white noise in timing data is described by fitting the parameters 'EFAC' and 'EQUAD' which can rescale the ToA error bars, and account for additional white noise, respectively (Edwards et al. 2006).Red timing noise is described using Fourier-domain Gaussian-process modelling as described in Lentati et al. (2014).The sinusoidal Fourier-basis components at each Fourier frequency  are multiplied by a set of amplitudes described by a Gaussian process; we refer to the covariance of these amplitudes as the power spectral density.In our fitting, we model the power spectral density prior as a power-law, as often used in the pulsar community (Lentati et al. 2014;Van Haasteren et al. 2011).However, for this study we choose to characterise this power-law by an amplitude at a reference frequency of 0.1 yr −1 , different from the usual choice of 1 yr −1 ; this choice ensures the amplitude is more robust against the yearly periodicity of the solar wind.We therefore characterise the power spectral density model as where A and  correspond to the Bayesian hyperparameters characterising the red timing noise, while the individual Gaussian-process amplitudes are marginalised over.This is used for both achromatic red timing noise, and chromatic DM noise, which has an inversesquare dependence on observing frequency, as per Eq. 1.For the former, the power amplitude is often written as A ≡  2 red /(12 2 ), and log 10 (  red ) is chosen as the fitting hyperparameter together with the exponent  ≡  red .For the DM noise, the Fourier components are also proportional to the square-inverse of the observing frequency  obs ; the prior power amplitude is expressed as A ≡  2 DM , and log 10 (  DM ) and  ≡  DM are the fitting hyperparameters.

Solar wind fitting using Gaussian processes
We present an implementation of time-variable amplitude spherical SW modelling which has been incorporated in the run_enterprise package and can therefore be included in the simultaneous Bayesian fitting of pulsar timing data, together with the other deterministic pulsar parameters, and the white and red noise parameters.
An approach based on Gaussian processes is used, similar in concept to that used to describe the red noise model as presented in Lentati et al. (2014).The SW signal is expressed, in the time domain, as the sum of independent components, which correspond to the unit-amplitude spherical SW for each solar conjunction, multiplied by a set of amplitudes equivalent to the quantity  SW e from Eq. 3. We choose a set of simple piecewise linear ('triangular') functions, with centre points at the solar conjunction times of the pulsar, to mathematically represent each yearly SW variation; this follows the idea presented in Keith et al. (2013) for modelling DM noise in general.In practice, the SW contribution to the vector of ToAs is written as where  is a column vector containing the spherical SW amplitudes equivalent to  SW e for each of the  c solar conjunctions within the data.V is a matrix of size ( c ×  t ), where  t is the number of ToAs, i.e. the size of the column vector  sw .We define an element of the matrix V, such that where  ∈ {1, . . .,  c } is the matrix row and  ∈ {1, . . .,  t } the column;  ,  is defined as per Eq. 5 for the solar elongation angle   corresponding to the ToA   .  is the time of the th solar conjunction in the data, and Λ is the triangular function, such that by definition Fig. 1 illustrates the steps in creating the SW components in this model.
We assume the distribution of SW amplitudes to be characterised by a Gaussian function, of mean value  0 , and a standard deviation  sw .In practice, we therefore express the vector of amplitudes  as a sum where  0 is the mean SW amplitude (i.e. the equivalent of the standard tempo2 parameter 'NE_SW'), and 1 is an  c -point column vector of ones. is therefore the column vector of SW amplitude variations away from the mean.We assume the elements of  to have mean zero, and their variation with each solar conjunction to be characterised by a Gaussian process, and therefore by a corresponding covariance matrix.In this analysis, we use a simple constant-variance, such that the covariance matrix of the SW amplitude variations  is a diagonal matrix with all diagonal elements equal to the variance  2 sw .The quantity  sw therefore sets the prior for the amplitudes , and is a hyperparameter referred to as 'SW_sigma' in the fitting code, and in  the output parameter file.Note that one can also choose a different, more complicated covariance informed by physical processes, e.g. one that follows the 11-yr solar cycle.While this may be explored in the future, we believe that the simple approach taken here is suited for capturing the year-by-year variability in amplitude in the current datasets.
Similarly to the approach taken in Fourier-basis red-noise fitting, we marginalise over the Gaussian-process amplitudes  and only keep  sw as a Bayesian fitting hyperparameter.As this model keeps the parameter estimation algebraically linear, including the additional time-variable spherical SW model when fitting a pulsar dataset is computationally inexpensive, as it only adds one extra hyperparameter to the entire timing model.Furthermore, tempo2 can be used after the Bayesian fit performed with the run_enterprise software package to explicitly find the amplitudes of the spherical SW at each solar conjunction, if of interest.This is done using the implemented constraining of the least-squares fitting, which takes into account the fitted  sw ; for an explanation of how this works, see Appendix A in Keith et al. (2013).

General setup
To test the capabilities of the SW model presented in this work, we created sets of simulated ToAs of both uniform and PTA-like cadence.For all simulations, we started from a set of 'idealised' ToAs (characterised by zero residuals) for a chosen pulsar, and added realistic levels of noise, informed by typical values in the observations available.This was done using tempo2 plugins, which generate the types of signal discussed in Section 3.1, namely: • addGaussian to add white noise; • addRedNoise to add achromatic red noise, characterised by a power-law prior of user-specified amplitude and slope; • addDmVar to add DM noise, also characterised by a power-law prior of user-specified amplitude and slope; • addArbitraryDM to add the DM influence from a simulated SW (see Section 3.2.2 for more details on simulating realistic SW DM series).
Using the above types of signals, multiple ToA sets were produced, serving various testing purposes which we describe in Section 4. Multiple Gaussian-process realisations were created for the same characteristic hyperparameters, such that we were also able to check the robustness of our pipeline with repeat measurements.

Simulating the solar wind
We simulate the SW influence on a pulsar as a DM time-series that is then added to the total simulated ToAs using the tempo2 plugin addArbitraryDM.We consider the variation with respect to the unit-amplitude ( SW e = 1) spherical SW, i.e. the ratio of the DM series to   , and sample it as a Gaussian process characterised by an exponential-squared kernel (see e.g.Rasmussen & Williams 2006), of the form where  represents the 'distance' between two observing times,   is the kernel amplitude, and  is the metric, i.e. the scale of the correlation within the signal described by this kernel.The exponentialsquared kernel is a somewhat arbitrary choice, but this simple stationary kernel is widely used and well-suited for describing smooth functions characterised by a single overall metric.In practice, we choose  = 400 d; the exact value of this is not strictly relevant, but a value slightly above 1 yr ensures that the simulated signal at each solar conjunction is roughly independent of that at the other solar conjunctions.
The amplitude   broadly characterises the overall variation away from the mean of the signal.To find a realistic value of   , we use the SW DM series of PSR J0034−0534 as presented in Fig. 3 of T21.We refer the reader to Tiburzi et al. (2019Tiburzi et al. ( , 2021) ) for details on how this was obtained.In short, using 5 yr of LOFAR observations, the total DM contribution at each average epoch was estimated.The contribution of the IISM to this DM series was then modelled by a cubic spline with each piece corresponding to a solar conjunction; the SW contribution was simultaneously modelled assuming a spherical SW as in Eq. 3, while the amplitude was allowed to vary year-by-year.In this work, the DM series after subtracting the cubic spline of the IISM model was used, and we refer to this as DM SW Lo .We define a time series  0 ≡ DM SW Lo /  , such that  0 only encompasses the estimated variation of the SW with respect to a spherical SW of time-invariant, unit-amplitude.Note that at a time   of the th solar conjunction in the data, the value of  0 is equivalent to the SW amplitude   .The time series  0 is then assumed to be described by a smooth function, characterised by an exponential-squared kernel as in Eq. 11, of metric  = 400 d.To find the representative amplitude   , we employ the python library george to fit this exponential-squared kernel to the  0 series, and find a maximum-likelihood value of   ≈ 3.7 cm −3 .
With the fully parameterised exponential-squared kernel, we are able to draw as many new samples of time-series  (of the same type as  0 ) as needed, for any pulsar with an available ephemeris.The unit-amplitude spherical SW (i.e.  ) of any known pulsar is easily obtainable from the tempo2 plugin 'general2' using the pulsar's ephemeris, such that samples of SW DM series for any known pulsar are simulated as DM SW sim ≡   .Fig. 2 shows an example of two simulated SW DM series obtained in this way.Note that because our method allows the spherical SW model to vary in time, the resulting simulated DM series is not strictly symmetric with respect to the solar conjunction peaks.The sampled series  is not a discrete set of amplitudes for each solar conjunction, but rather a continuous function that varies with each observation time in a way informed by real data.

RESULTS & DISCUSSION
In this section, we present and discuss the results of our tests of the pipeline and of the sensitivity of available data to SW variations.In all cases, we use the software package run_enterprise to simultaneously fit for the standard pulsar deterministic parameters, as well as achromatic and DM red noise with power-law priors.The deterministic parameters are marginalised over, and recovered later if necessary using least-squares fitting in tempo2; the hyperparameters for the red-noise power-law priors (log 10 (  red ),  red , log 10 (  DM ),  DM ), and the square-root of the variance in SW yearly amplitudes ( sw ) are sampled using a Markov Chain Monte Carlo technique through the python package emcee (Foreman-Mackey et al. 2013).In this context, we refer to 'time-invariant' SW fitting as only fitting for an overall mean amplitude  0 (i.e.'NE_SW' in tempo2) of a spherical SW model.This is in contrast to the model described in Section 3.1.2,in which we fit for a 'time-variable' spherical SW.In practice, this time-variable spherical SW is described by yearly amplitude variations  away from an overall mean amplitude  0 , and we fit for the  elements, as well as for  0 .

Testing the solar wind fitting code
To initially probe the capabilities of our pipeline, we simulate a 20-yr long set of uniform, high-cadence (one observation every 10 days) ToAs of PSR J0034−0534.This testing dataset includes simultaneous observations of three common frequency bands, centred at 300 MHz, 1440 MHz and 2400 MHz, with an rms of 1 µs on each ToA.While these data are 'optimistic' compared with typical pulsar datasets, they provide a useful first test to gain insight into the performance of the SW fitting pipeline.To create the simulated dataset, we start from idealised ToAs (of zero residuals), and 'inject' white noise variations at the 1-µs level, Gaussian-process achromatic red noise and DM noise with known power-law priors informed by typical real levels of noise seen in the observations, as well as a realistic time-variable SW, using the method described in Section 3.2.2.Multiple realisations of the noise and SW Gaussian-process samples are used for robustness and repeat-measurement checks on the results of the pipeline.
In this section, we aim to explore the proficiency of our pipeline.Firstly, the injected signals are compared to the recovered ones using either the time-invariant or the time-variable SW method.The comparison of the results of the two fitting methods is discussed in Section 4.1.1.Secondly, we use the measurements of 2000 SW amplitudes in order to validate the size of the uncertainties produced by our pipeline; this is described in Section 4.1.2.

Comparison between time-invariant and time-variable fitting
One of the benefits of using simulations is that the quantity of red noise in the pulsar ToAs is known, and we can therefore compare the properties of the recovered signals to those of the injected ones.Here, we run both the time-invariant and the time-variable fitting pipelines on the same set of simulated ToAs, which include both achromatic and DM red noise.
In general, we expect the recovered DM power-law to be affected by whether the SW influence is well modelled.Specifically, if only a time-invariant mean amplitude is used to model the SW, any yearby-year variation is likely to be absorbed by the DM noise model instead.This leads to excess power in the measured DM noise at high Fourier frequencies (at 1 yr −1 and further harmonics), which in turn causes the power-law to appear flatter than the real power-law process.Therefore, unmodelled SW influences may bias the interpretation of the DM spectrum.Furthermore, this could also bias the estimation of other timing parameters, and especially their uncertainties, as the DM power-law would not accurately represent the actual noise present in the data.
To investigate the performance of our pipeline compared to the time-invariant fitting in this context, we inspect the low-frequency (300 MHz) residuals after subtracting the maximum-likelihood SW contribution and fitted achromatic noise.Fig. 3 shows these residuals, overlaid with the injected IISM DM variations.When fitting only for a time-invariant SW, we can see additional 'spikes' in the residuals around the solar conjunction times, a result of absorbing unmodelled SW into the DM variations, whilst the time-variable SW fitting is fully consistent with the injected IISM contribution.
The average power spectral densities of these recovered DM influences, as well as the injected power-law prior, are shown in Fig. 4. As expected, in the time-invariant SW fit, the absorbed SW influence creates additional power at high frequencies, and therefore flattens the power-law shape compared to the injected prior; while the timevariable fit recovers the injected power-law well.The shape of the average power spectral density in the time-invariant fit shows the effect of the SW 1/yr-frequency and its harmonics, modulated by subtracting a mean-amplitude SW, which creates the observed dips of smaller width.
The injected and recovered red noise hyperparameters for 15 realisations of the simulation are shown in Fig. 5, as defined in the power-law in Eq. 6.The achromatic red noise hyperparameters are in agreement with the injected values and equally well recovered in both cases; this may be expected, as the SW influence is an intrinsically chromatic effect.The time-variable SW pipeline also recovers the DM power-law well, though this is not the case for the time-invariant SW model.Notably the time-invariant model leads to a spectral exponent that is significantly flatter, with the recovered slope being of order 10-sigma smaller than the 'true' value.This appears to be a systematic rather than statistical effect, i.e. not dependent on the specific Gaussian-process sample used in the simulation, but rather due to the fitting model.The recovered SW model from each fitting method can also be compared to the injected SW influence.Fig. 6 shows this for the DM series (top panel), as well as the spherical SW amplitude  SW e (bottom panel).The modelled solar-conjunction amplitudes estimated by our pipeline follow the injected SW well, being within 2-sigma of all injected values.The time-invariant fit (shown only in the bottom panel of Fig. 6) only models the SW influence with a constant-in-time amplitude, and the uncertainty in the mean amplitude recovered is too small to account for the yearly variation.
We also consider the timing residuals after removing all the modelled signals, including the achromatic red noise and DM noise; these are shown in Fig. 7. Ideally these fully subtracted residuals should be as close to pure uncorrelated noise as possible.The size of this left-over noise, which can be characterised by its root-mean-square (rms), is also relevant in general when searching for signals in the data -such as, for example, the search for pulsar binary companions, or a GWB signature, for which the lower the noise the better the chance of detection.Of the two fitting modes used in this test, the yearly variable SW appears to be a better model, as strongly supported by the Bayesian evidences, which yield a natural log-Bayes factor of 160 in favour of the time-variable pipeline; the Bayesian evidences are computed using the nested-sampler dynesty (Speagle 2020).This effect can also be seen in the fully subtracted residuals (Fig. 7), where the residuals obtained using the time-variable SW fitting qualitatively appear more 'white', lacking the unmodelled SW 'spikes'.However, note that the reduced- 2 values of the two whitened time series are indistinguishable (∼ 1), likely due to the large number of data points not at solar conjunctions.The level of white noise, quantified here by the rms, is better by ∼ 1 µs for the residuals obtained from the time-variable fitting than those from the time-invariant fitting.

Recovered uncertainties
To check the robustness of our results and their uncertainties, 100 sets of ToAs are simulated, all with the same general properties as the dataset used in Section 4.1.1.Each of these simulated datasets is created from different Gaussian-process realisations of the same power-law prior achromatic red noise and DM noise, as well as different samples of SW variations.The results of running these 100 independent datasets through our pipeline are consistent with the findings in Section 4.1.1 with respect to the injected values.
The 2000 solar conjunctions available in total in the 100 datasets are used to study the statistical properties of the estimated uncertainties on the SW amplitudes.Each measurement of the SW amplitude is normalised with respect to the injected value and the associated uncertainty, such that the equivalent 'standardised variable' is computed.The population of these standardised variables is expected to be normally distributed, with a mean of zero and a standard deviation of one, if the distribution of the measured variable is Gaussian.In the case of a set of uncorrelated measurements, {  }, assumed to be Gaussian distributed, the standardised variable is simply defined as (  − )/ for each measurement, where  and  are the mean and standard deviation of the Gaussian distribution.In this analysis, however, the measurements of SW amplitudes are correlated through the mean  0 , such that we use linear algebra to compute the standardised variable values, as follows.
If  is the column vector of the  = 2000 measured SW amplitudes, and  inj is the equivalent column vector of injected SW amplitudes, then we can write where  is a column vector representing the random variation of  around  inj , expected to be normally distributed with mean zero and variance according to the parameter covariance matrix C of the elements in .If C is known, the vector  =  −  inj can be 'whitened' to disentangle the correlations between different amplitude measurements, and therefore to obtain the column vector of standardised variables .This approach is similar to e.g. the whitening presented in Coles et al. (2011) for general pulsar data, and we summarise the steps below.Using the Cholesky decomposition, the parameter covariance matrix is written as C = LL ⊤ , where L is a lower triangular matrix with a real and positive diagonal.The whitening process then yields that the column vector of standardised variables can be computed as Recall that in our pipeline, we fit for the mean SW amplitude  0 and the deviations from it at each solar conjunction, i.e.  =  −  0 1, where in this case 1 is an -point column vector of ones.The relevant column vector of measured SW parameters is therefore given by the ( +1)-point column vector  = ( 0 , ) ⊤ , with the associated parameter covariance matrix, C  .The parameter covariance matrix   is a direct result of the timing analysis, and can be obtained directly from tempo2 post-fitting.The whitening covariance matrix C can then be computed as where M is the transformation matrix from the (+1)-point column vector  to the -point column vector , i.e.  = M , and therefore has the shape where the first column contains exclusively '1's and the rest of the matrix has the shape of an ( × ) identity matrix.Therefore, the standardised variable  can be computed as per Eq. 13 using the Cholesky decomposition of the parameter covariance matrix C, which can be estimated from the output covariance matrix of the analysis, C  , using Eq.14 and the shape of the transformation matrix, M, as given in Eq. 15.
Fig. 8 shows a histogram of the (normalised) distribution of the measurements of  for the 2000 solar-conjunctions amplitudes fitted with our pipeline.The distribution of this rescaled uncertainty is well described by a Normal distribution, with an Anderson-Darling test equivalent -value of 0.55 (e.g.D'Agostino 1986).The Normal distribution has mean of −0.05(2) and a standard deviation of 1.09(2), not meaningfully distinguishable from a unit-variance zero-mean Normal distribution; we therefore conclude that our uncertainties are well estimated within our assumptions.

Solar wind fitting in PTA-type simulations
In this section, we aim to establish whether PTA-like data are sensitive to yearly variations in SW away from the mean amplitude.For this purpose, EPTA+InPTA data of PSR J1022+1001 are used as a basis for simulated ToAs which are then analysed through our pipeline.We use the real cadence and uncertainties of this dataset to simulate current observations (Section 4.2.1), as well as median properties of relevant telescopes to approximate future observations of the EPTA+InPTA (Section 4.2.2).

Current EPTA + InPTA observations
To simulate a set of ToAs with properties characteristic to the EPTA+InPTA dataset, we start from the real data of PSR J1022+1001, which shows the highest SW influence of the 25 pulsars included in the EPTA DR2 (as this pulsar has the closest approach to the Sun, having an ecliptic latitude of −0.06 • ).Idealised ToAs (of zero residuals) are created such that the observing cadence, frequencies, and error bars of the real data are preserved into the simulated ToAs.As before, realistic levels of noise and SW variations, based on observations, are added to these idealised ToAs.
The comparison of the measured red-noise hyperparameters with the injected values is presented in Fig. 9 for both the time-variable and the time-invariant SW fitting.It can be seen that there is no significant difference between the results of the two fitting modes in this case, and that the best-fit values are mostly consistent with the 'true' values within their uncertainties.However, the precision of these measurements is, as one would expect, worse than in our previous simulations of uniform cadence and low-frequency, discussed in Section 4.1.1.Fig. 10 shows the fitted SW amplitudes through our pipeline, compared to the injected Gaussian process sample.While our results are generally consistent with the injected values, the uncertainties are large enough that the results for this dataset are also consistent with a constant-in-time SW.This outcome is also confirmed through an examination of the residuals after subtracting both the time-variable fit- ting model and the time-invariant model, which reveal no significant difference between the properties of the two: the reduced- 2 ≃ 0.97 in both cases; while the rms values are 1.37 µs and 1.33 µs for the time-invariant and time-variable fitting, respectively.We therefore conclude that the current EPTA+InPTA dataset is not sensitive enough to measure yearly changes in the SW amplitude.However, we also note that there is no noticeable disadvantage of including this additional model in the fitting, as it adds no significant time to the computational run, nor does it increase the overall left-over white noise level.Moreover, when there is not sufficient data around a solar conjunction to measure the variation in SW amplitude, the total amplitude value defaults to the mean (as can be seen especially in the early data in Fig. 10).The lack of sensitivity to measuring SW changes of even the more recent EPTA+InPTA data is indicative of the small number of observations at low frequencies (≲ 1000 MHz) -where the SW influence is larger and therefore easier to quantify -, as well as of the relatively small fractional bandwidths of the observations.However, this is expected to change in the future, as more uGMRT (InPTA) data, which includes high-quality simultaneous observations at low frequencies (∼ 400 MHz), is combined into the EPTA dataset.Note that the current data combination, as used in the simulations in this section, includes some InPTA data, particularly around the solar conjunction around MJD 58360; these InPTA observations are indeed seen to have a plausibly beneficial effect on SW amplitude fitting for that year.Furthermore, the upcoming inclusion of the high-cadence and very low-frequency LOFAR data into the EPTA dataset will unlock next-level sensitivity to the SW.

Simulated 'future' EPTA+InPTA observations
With the results of the realistic (EPTA+InPTA)-like simulations in mind (as presented in Section 4.2.1),we aim to investigate the sensitivity of this dataset after 10 years of observations at the current observing setup.For this purpose, we create sets of simulations based on the median of recent properties of the LEAP observing system (which is the most sensitive L-band 'observatory' in the EPTA dataset), and on the median of recent properties of uGMRT observations.The properties of these simulated observations are summarised in Table 2, showing the number of frequency sub-bands for each observing system, the total frequency range, the (uniform) observing cadence, and the (uniform) rms characterising the system white noise.Note that the observations with the two uGMRT frequency bands (B3 and B5) are simultaneous, which is particularly useful for studying chromatic effects.Uniform-cadence observations are used in this case for simplicity; while the sensitivity of any pulsar data to the SW indeed should depend on the density of observations close to its solar conjunction, a detailed study on how this cadence affects measurements and on determining the best observing strategy is beyond the scope of the current work, and is left as a future investigation.
Here, we consider the improvement of the noise budget of the observations when using the time-variable amplitude SW fitting instead of the time-invariant SW fitting.Fig. 11 shows the residuals in both cases, after subtracting all the fitted signals, including the  11, but the residuals are plotted against solar angle, i.e. the sky distance between the pulsar line-of-sight and the Sun.Again, the plots on the left correspond to the time-invariant SW amplitude fit, while the plots on the right show the results of the time-variable fit.The top plots show the residuals for every individual data point, while the bottom plots show the rms of these residuals in bins of 5 • solar angle.Of particular interest are the points at low solar angle, i.e.where the pulsar is close to solar conjunction and the SW influence is the largest.red noise processes.The time-variable fit generally performs better in whitening the pulsar ToAs.Qualitatively, there are obvious signatures of the SW that are not removed in the time-invariant case, particularly near the solar conjunction close to MJD 62740, visible both in low-frequency uGMRT data and, to a lesser extent, in Lband LEAP observations.Quantitatively, while the reduced- 2 ≃ 1 in both cases as before, the overall left-over noise level is reduced by using our fitting compared to the time-invariant mean: for the full simulated dataset, the rms is reduced by nearly 1.3 µs; while just for the L-band LEAP data, the rms is improved by 160 ns.The plots of Fig. 12 show the same LEAP residuals, but as a function of the solar angle of the pulsar, i.e. its angular distance to the Sun as projected on the sky, for both the time-invariant and time-variable fitting.These plots show that the main difference between the two fitting methods is, as expected, due mainly to those observations closest to the Sun, although even some observations at tens of degrees of solar angle are affected.
Overall, we conclude that in the near future, as more InPTA observations will be combined in the EPTA datasets, the sensitivity to the SW in pulsar data will increase.Moreover, by fitting this data for a time-variable spherical SW, the noise budget is likely to improve, even for L-band observations.This is particularly interesting in the context of e.g.detecting the GWB by PTAs, where a decrease in noise of just a few hundred nanoseconds could be valuable in reaching the target detection significance in light of the recent GWB results (EPTA Collaboration et al. 2023a,b,c).

Real data
We also test our pipeline on real data, of PSRs J0030+0451, J1022+1001, and J2145−0750, respectively.Firstly, our pipeline's SW amplitude estimates from independent datasets taken by the observing systems EPTA+InPTA, PPTA, and LOFAR are compared in Section 4.3.1.Secondly, the recovered DM time-series from the fits on LOFAR data are compared to those given in T21, which used the same (while slightly shorter) LOFAR datasets, but a different method; we are therefore able to directly compare SW results between our pipeline and another, independent method; this is discussed in Section 4.3.2.

Comparison between results from different datasets
Not all telescope datasets were available for each of the three pulsars chosen for this analysis; the data used is summarised in Tables 1  and 3.All datasets are run through the same pipeline, fitting for deterministic parameters, the SW model, white noise, achromatic red noise, and DM noise simultaneously.We compare the yearly SW amplitudes estimated from each separate dataset, for each pulsar.This is shown in Figs. 13,14,and 15, while the mean values, as well as the median and minimum of the SW amplitude uncertainties for each dataset are also presented in Table 3.The values of SW amplitudes from the independent LOFAR-data analysis from Fig. 6 of T21 are also included here for comparison, but are mostly discussed in the following Section 4.3.2.
Firstly, the results for PSR J0030+0451 are shown in Fig. 13.There is agreement within uncertainties between the results of the EPTA and LOFAR datasets.As is perhaps expected, the much-lower frequency LOFAR data performs significantly better in estimating the SW amplitudes.
Secondly, Fig. 14 shows the results for PSR J1022+1001.In general, the SW amplitude measurements are consistent between the different observing systems.While the EPTA dataset is a few years longer than the others, the early data are sparse enough that only rough estimates of the amplitudes are possible.Broadly, the PPTA data allow for more precise measurements of the SW than the EPTA data, likely owing to the additional lower-frequency and wider-bandwidth observations in the PPTA dataset.Further, the LOFAR data are seen to be much more sensitive to these measurements than the PPTA data.The collection of InPTA data has started relatively recently, and as such only two amplitude measurements are supplemented by this in the EPTA+InPTA data combination; these are the last data point and that near MJD 58360.The increased sensitivity of the amplitude measurement which included both EPTA and InPTA data for the solar conjunction near MJD 58360 hints at the advantage of the additional InPTA data.
Thirdly, the SW measurements for PSR J2145−0750 are shown in Fig. 15.The early data (before 2004) from the PPTA represent archival data with a low cadence, such that no variation away from the mean-amplitude SW can be measured.The precision of the LO-FAR measurements is significantly better than that of the PPTA measurements in the case of PSR J2145−0750 as well.We also note that the recent, roughly 4 years of InPTA data show an improved precision compared to the PPTA dataset for this pulsar, likely due to their high-cadence simultaneous observations.The observed dip in SW amplitude occurring around the year 2009 (around MJD 55000) for PSR J2145−0750 is coincident with the independently observed low in solar activity (e.g.Hathaway 2015), such that our measurements appear to track the broad behaviour of the solar cycle.The same dip can also be seen, in fact, for one solar conjunction in PSR J0030+0451; however, this is not obviously present in the results of PSR J1022+1001.This suggests that, on one hand, the SW signature in pulsars further away from the ecliptic plane is likely to follow the broad 11-yr solar cycle behaviour.On the other hand, the shape of the variation away from the mean seems harder to predict for pulsars of very low ecliptic latitude (such as PSR J1022+1001).This may be caused, for example, by the slow SW contribution dominating the fast SW around the solar equator, even during the solar cycle maxima, thus showing less of the variation between the two solar cycle stages.The presence of persistent streamers near the solar equator could also induce the kind of DM variations seen in pulsars at very low heliospheric latitudes.Furthermore, the mean SW amplitude of observations of PSR J1022+1001, of ∼ (10 ± 1) cm −3 is significantly larger than that of PSR J2145−0750, of ∼ (6 ± 1) cm −3 .The very similar range of observing times for these two pulsars suggests that this divergence may be due, as before, to the different ecliptic latitudes, and the corresponding heliospheric latitudes and magnetosphere areas probed.Therefore we conclude that, in general, a global fit describing the SW properties of multiple pulsars is likely to be advantageous only if the diverse regions of the SW probed by pulsars at varying heliospheric latitudes are carefully considered.Tiburzi et al. (2021) Finally, the results of our pipeline are compared with the independent analysis in T21, using the same observational LOFAR data available for PSRs J0030+0451, J1022+1001, and J2145−0750.For our results, the data are processed through our pipeline, as with all the previous analyses presented in this work, and the yearly SW amplitudes are estimated.The full DM series are also reconstructed from the fitted parameters, namely as the sum between the SW model, the chromatic red-noise power-law Gaussian process, and the deterministic DM model; the latter is expressed as a polynomial

LOFAR results: Comparison with
where  is a ToA expressed with respect to a chosen epoch ('DME-POCH' in tempo2), and the three polynomial coefficients (DM 0 , DM 1 , and DM 2 ) are equivalent to the 'DM', 'DM1', and 'DM2' fitting parameters in tempo2.
Conversely, the analysis in T21 was based on first obtaining a DM value per observation; for more details on how this was performed, we refer the reader to Tiburzi et al. (2019Tiburzi et al. ( , 2021)).To estimate the SW from the DM series, a reference value (equivalent to DM 0 ) was subsequently subtracted from the DM series.The remaining contributions of the IISM and SW were disentangled and simultaneously modelled using a Bayesian framework, for each specific segment of data corresponding to a solar conjunction.The IISM contribution was modelled as a cubic polynomial for each solar-conjunction segment, and continuity between different segments was insured.In short, the total DM series in the T21 analysis was modelled as a sum between the reference value DM 0 , the IISM cubic-polynomial, and the spherical SW of yearly variable amplitude.Note that the DM series used here for comparison with our full reconstructed DM is the one initially obtained from the observations by T21, prior to the SW analysis.Figs. 16,17,and 18 show the results of our analyses, as well as those of T21 for comparison, for the three PSRs J0030+0451, J1022+1001, and J2145−0750.Note that more recent LOFAR data were available for PSR J1022+1001 since the T21 work, which are included in our analysis; the additional data which appears here within one solar conjunction either side were discarded in T21 as there were not enough ToAs to provide a robust estimate by the criteria chosen in that work.
From these plots, we conclude that our estimates of the SW amplitudes are for the most part consistent with those of T21, while using an independent Bayesian approach.We note that where there are discrepancies in the SW amplitudes, they appear to be caused by the difference in the estimated IISM contribution.For example, at the solar conjunction near MJD 58000 of PSR J1022+1001 (Fig. 17), our IISM estimate shows an overall flat behaviour, while the T21 analysis presented a higher, cubic-varying DM estimate.This emphasises an advantage of our pipeline with respect to the T21 method: the IISM contribution to DM is estimated for the entire dataset, rather than in segments, which in general ensures a smoother, and likely more realistic behaviour.We also note that, for the specific measurement near MJD 58000 for PSR J1022+1001, the observed difference in IISM estimates of roughly 10 −4 pc cm −3 would only correspond to a direct difference in the SW amplitude ( SW e ) of roughly 10 −2 cm −3 ; the actual measured difference is, however, of order ∼ 1 cm −3 .Nevertheless, we believe that, since this pulsar has a very small ecliptic latitude, the tails of the SW influence may impact the shape and peak of the SW model to a larger extent than would be trivially expected.
Lastly, it is important to acknowledge that the spherical SW model, even while allowing its amplitude to change each year, is not sufficient to fully describe the observed SW influence on these data, as also observed previously in e.g.Tiburzi et al. (2016).This is seen in our results, as the IISM contribution to the DM clearly absorbed some asymmetric features close to the solar conjunction, which are highly likely due to the solar influence.Therefore, any analysis assuming a spherically symmetric SW, which includes that presented in this work, would not allow for an entirely robust studying of the SW influence in pulsar data.The pipeline developed and presented here, however, can create a basis that can be straight-forwardly adapted in future work to improve on this, and include various additional models, such as e.g.different piecewise components for before and after the solar conjunctions, allowing the tails of the SW variation (i.e. the pulsar ingress and egress from the Sun) to change independently of each other.

CONCLUSIONS
We have implemented a computationally inexpensive, linear Gaussian-process piecewise Bayesian approach to fit pulsar ToAs for a spherical SW of yearly time-variable amplitude, simultaneously to all the other pulsar timing and noise parameters; this is available through the pulsar analysis toolkit run_enterprise.We have explored the functionality of this fitting pipeline using simulations, and found that it performs as expected, and better than the currently widely used time-invariant spherical SW, particularly in data of low frequency and large fractional bandwidths.
Using simulations, we found that the current EPTA+InPTA dataset is not yet sensitive enough to measure variations of the spherical SW amplitude.However, assuming that the current observing strategy of uGMRT will continue, future EPTA+InPTA data will have increased sensitivity to the SW, such that using our time-variable fitting could improve the rms (white-noise levels) of the residuals by at least ∼ 200 ns at L-band, which may help with the noise budget for extremely sensitive experiments such as the search for the GWB.Indeed, recovering the SW influence in EPTA+InPTA pulsar data will depend on the cadence of observations near the solar conjunctions, but a study of these effects is beyond the scope of this work, and is left as a future investigation.
We also applied our pipeline to real data of three pulsars, i.e.J0030+0451, J1022+1001, and J2145−0450, that are known to show the influence of the SW in their ToAs, and were also part of the previous study by T21.The SW amplitudes found from individual fitting of pulsar data from the EPTA+InPTA and/or the PPTA were compared with the results from LOFAR data through our pipeline, as well as the independent results from T21.This showed that SW amplitudes found from different datasets were mostly consistent with each other and that, as expected, where lower-frequency and wider-bandwidth data were present, the uncertainties were reduced.Furthermore, the variation in the fitted SW amplitudes for the pulsar of the highest ecliptic latitude in this study (J2145−0450; elat = 5.31 • ) roughly followed the 11-yr solar cycle; while for the pulsar of very low ecliptic latitude (J1022+1001; elat = −0.06• ), the amplitude variation generally did not seem to correlate with the long-term solar cycle.

Fit ToAs Tiburzi+2021
Figure 16.The recovered DM series and yearly SW amplitudes from our analysis of the LOFAR ToAs (green) are compared with the same quantities as presented by T21 for the same LOFAR data (purple), for PSR J0030+0451.While the top plot shows the broad picture of the DM series, the middle plot is a zoomed in version of the same data for a clear comparison.For both our results (in green), and the T21 results (in purple), the data points represent the full estimated DM series -which includes the polynomial terms, the Gaussian-process power-law ('PL'), and the SW.The continuous lines represent only the estimated IISM contribution; this was fitted as an average DM 0 plus multiple consecutive cubic terms in T21, and as a quadratic plus a Gaussian process power-law in our analysis.The bottom plot simply illustrates the SW amplitudes, as also shown in Fig. 13, but only for LOFAR data.
This hints that a global SW fit may be more beneficial if e.g. it is performed in slices of ecliptic latitude.A larger study would clarify if this is indeed a wide systematic effect.
The DM series, including the SW effect, were also compared for the same LOFAR data as recovered from our pipeline and as used in T21.We found that where there was a difference between our SW amplitude and that in T21, this was likely due to the estimated IISM background; we believe our simultaneous fitting of deterministic and noise pulsar components produces a more plausible shape for this background than the consecutive cubic fits as used in T21.However, we conclude that even with the addition of the time-variable amplitude, the spherical SW model is not enough to fully account for the SW influence in LOFAR data.The left-over SW absorbed in the background DM series suggests that perhaps using an asymmetric model, with different components pre-and post-solar conjunction may be beneficial.The piecewise framework presented in this work could be modified for this, or for an alternate way to further study SW models.

Figure 1 .
Figure 1.Graphical illustration of the mathematical representation in Eqs.7-9.The quantities in the top three plots (  , Λ, and ) are multiplied to give the DM (bottom plot).Note that we do not show the dependence on observing frequency for simplicity, and DM =  D  2 obs V.Each colour represents a different component, corresponding to each solar conjunction,   .The data points represent the ToAs   .

Figure 2 .
Figure 2.An example of two realisations of simulated SW using an exponential-squared kernel with amplitude informed by real data.The data points are at uniformly sampled ToA values.The top plot represents the simulated series , which is then multiplied by the pulsar-specific geometrical factor   to create a DM series, shown in the bottom plot.

Figure 3 .Figure 4 .
Figure 3.A single realisation of simulated residuals for the 300 MHz band, after subtracting the fitted SW variations and achromatic red noise when using the time-invariant (left) and the time-variable (right) SW amplitude model.The solid black line shows the delays at 300 MHz for the injected IISM DM variations.Solar conjunction times are shown by the vertical gray lines.

Figure 5 .Figure 6 .
Figure 5. Measurements obtained from 15 realisations of the uniform simulation as described in Section 4.1.1 are shown.Specifically, these plots show the power-law log-amplitude (log ) and slope () for the DM (top plots) and achromatic red noise (bottom plots).The error bars represent one standard deviation.The vertical dashed lines show the 'true' values used in the simulations.In each plot, 15 realisations of the simulation are shown for each of the two fitting modes: the top red circle points represent the measurements obtained from the time-invariant SW fit; while the bottom green squares correspond to the time-variable SW fitting as described in this work.

Figure 7 .
Figure 7. Post-fit residuals after subtracting the entire timing model, including achromatic and DM red noise, are shown for both the time-invariant (left) and time-variable (right) SW amplitude model.Only the data at the 300 MHz observing frequency is shown for clarity, as it shows the largest difference between the two models.The root-mean-square (rms) values for both the 300 MHz and the full data ('all') are also presented.

Figure 8 .
Figure 8.The histogram of the standardised variable () obtained from the 2000 recovered solar conjunction amplitude values and uncertainties.The best-fit Gaussian distribution is also shown, with a mean of −0.05(2) and a standard deviation of 1.09(2).

Figure 9 .Figure 10 .
Figure 9. Measurements obtained from 15 realisations of the (EPTA+InPTA)like simulation as described in Section 4.2.1 are shown.See the caption of Fig. 5 for details.

Figure 11 .Figure 12 .
Figure11.All plots show post-fit residuals after removing the fitted achromatic and DM Gaussian-processes red noise.The plots on the left correspond to the results after using the time-invariant SW amplitude fit, while those on the right correspond to the time-variable fit.The top plots show the full simulated data, with the uGMRT-type data in orange and the LEAP-like data in blue.The bottom plots show the same information, but zoomed in and just for the LEAP data.The weighted root-mean-square (Wrms) is given in each case for all the data shown in the plot.

Table 3 .
Table summarising the uncertainties in the measured SW amplitudes ( SW ), as well as the mean values ( 0 ) from the real datasets of the three pulsars discussed, as shown in Figs. 13, 14, and 15.The median and minimum values of these measured uncertainties (u.) are quoted in columns 3 and 4, respectively.The values in square brackets represent the median estimated only using 'newer' data, where relevant, i.e. after 2004 for PPTA, and after 2005 for EPTA.Column 5 summarises the mean SW amplitude values, with the values in brackets representing one standard deviation.

Figure 13 .Figure 14 .
Figure13.The yearly SW amplitudes for PSR J0030+0451 are shown, resulting from the EPTA (blue dots) and LOFAR (green squares) data, as well as the corresponding values presented in T21 (purple diamonds).The cadences of the observations used are also shown, where 'E' stands for the EPTA data, and 'L' for the LOFAR data.Note that the two x-axes are the same, but presented in both MJD and Year for convenience.

Figure 15 .
Figure15.The yearly SW amplitudes for PSR J2145−0750 are shown, resulting from the InPTA (blue dots), PPTA (orange crosses) and LOFAR (green squares) data, as well as the corresponding values presented in T21 (purple diamonds).The cadences of the observations used are also shown, where 'In' stands for the InPTA data, 'P' for the PPTA data, and 'L' for the LOFAR data.Note that the two x-axes are the same, but presented in both MJD and Year for convenience.

Table 1 .
Ecliptic latitudes (elat) and summary of which datasets are available for the three pulsars used in this work; a tick (✓) signifies the data are available, whereas a cross (✗) means they are not.By EPTA availability we refer specifically to the Data Release 2 dataset.The EPTA and InPTA data are combined for PSR J1022+1001.

Table 2 .
The properties of the uniform-cadence simulations based on median real-data values.