The impact on astrometry by solar-wind effect in pulsar timing

Pulsars are highly magnetised, fast-rotating neutron stars, and at the same time exquisite tools allowing for a variety of astonishing astrophysical studies. In many of these experiments, pulsar astrometry is an essential input. In particular, accurate knowledge of pulsar distances is often a key determinant of the precision of the experiment. For instance, in gravity experiments with binary pulsars, external impact such as the Doppler effect caused by the relative pulsar–solar system motion (Shklovskii effect, Shklovskii, 1970) and gravitational potential in our Galaxy, needs to be corrected to accurately recover the intrinsic variation of the binary orbit (e.g., Liu et al., 2020; Kramer et al., 2021; Guo et al., 2021). The calculation of these effects requires the knowledge of a pulsar’s distance. In a pulsar-timing-array (PTA) experiment to detect gravitational waves (GWs), the precision of pulsar distances is required to be a fraction of the GW wavelength, typically a few parsecs or even smaller, so as to be able to use the ‘pulsar-term‘ signal within the PTA response to GWs (e.g., Lee et al., 2011; Babak et al., 2016). This capability will significantly enhance the sensitivity of a PTA, especially to continuous GWs from individual sources. Additionally, pulsar astrometry is crucial for understanding the structure of our Galaxy, e.g., for constructing the model of Galactic free electron (Cordes & Lazio, 2002; Yao et al., 2017) and magnetic field distribution (e.g., Han et al., 2018). Furthermore, astrometric measurements from pulsar timing can be used to compare the reference frame defined by the JPL planetary ephemeris with that employed by the International Celestial Reference Frame (e.g., Wang et al., 2017), which is fundamental to the comparison and translation of stellar distance measurements based on these two systems (Madison et al., 2013).

Typically, there are two primary ways to precisely measure pulsar astrometry. The first method involves the canonical very-long-baseline-interferometer (VLBI) technique which delivers precise localisation of the pulsar in the sky (e.g., Brisken et al., 2002). Accordingly, pulsar astrometry is achieved by tracking the apparent annual positional variation of the pulsar over a period of one or a number of years. This approach enables precise astrometric measurements for different classes of pulsars including canonical pulsars, millisecond pulsars (MSPs) and magnetars (Deller et al., 2016, 2019; Ding et al., 2020). For a pulsar at a typical distance of 1 kpc, the positions and the annual parallax (which is the inverse of the pulsar distance) can be measured with a precision of 0.01–0.1 milli-arcsecond (mas; Ding et al., 2023).

The other method involves modelling the annual parallax of the pulsar via precision pulsar timing (e.g., EPTA Collaboration et al., 2023a). This timing technique employs the pulsar as a precision clock in its inertial frame. It first measures the times-of-arrival (TOAs) of the pulsar signals at Earth, and then converts the measurements to the pulsar’s proper time. The time transformation formula, typically referred to as the timing model, consists of a sequence of time-delay terms that are functions of the pulsar’s rotational, astrometric, and orbital parameters. For timing experiments in the radio band, this model also accounts for external propagation effects, such as the frequency-dependent dispersion delay of the pulsar radio signal caused by the ionized interstellar medium. In pulsar timing, the position and annual parallax of the pulsar are modelled through geometric time delays. These delays are caused by the annual orbital motion of the Earth within the solar system and the curvature of the pulsar signal’s wavefront due to the finite distance of the pulsar, respectively (e.g., Edwards et al., 2006). For a canonical pulsar, its position can be measured with a precision of 10–100 milli-arcsecond (e.g., Wang et al., 2017). Meanwhile, the parallax time delay is only secondary to that caused by the Earth’s positional variation (Roemer delay). For a typical parallax of 1 mas, the parallax time delay amounts to several hundred nanoseconds. As a result, precise parallax measurements are feasible only for pulsars with high precision timing, predominantly MSPs. Typically, the position and parallax of an MSP can be measured with a precision of 10–100 micro-arcsecond ($\upmu$as) (e.g., EPTA Collaboration et al., 2023a).

However, pulsar astrometry derived from timing analysis may be affected by other effects that exhibit similar and partially correlated features in the timing data. Typically, timing parameters are determined simultaneously by fitting the entire timing (and sometimes the noise) model to the dataset. This methodology inherently captures the correlations among various effects present in the timing data. Hence, the measurements of astrometric parameters may be significantly biased if correlated effects are not properly accounted for. For instance, Madison et al. (2013) demonstrated how red noise in the timing data can compromise the accuracy of derived astrometric parameters. Their study showed that applying a whitening filter to the data noticeably reduces such biases, leading to more reliable astrometric measurements.

Here, we present a detailed investigation into the impact on astrometric measurements by the solar-wind effect, with a focus on the timing parallax. It is known that the solar wind produces free electrons, altering the local free electron density in the solar system and thus causing an additional frequency-dependent dispersion delay in the radio emissions from pulsars (e.g., Iwamoto et al., 1998). Substantial analyses have been carried out using pulsar timing as a sensitive tool to probe the solar-wind properties (e.g., Tiburzi et al., 2019; Shaifullah et al., 2020; Tiburzi et al., 2021). Typically, in pulsar timing the solar-wind effect is described by a spherical model parameterized by the number density of free electrons at a 1-AU distance from the Sun ($n_{\rm sw}$, e.g., Edwards et al., 2006). In many previous studies, the pulsar’s position and parallax were fitted in the timing analysis to obtain the pulsar’s distance, when $n_{\rm sw}$ was used as a fixed constant input (e.g., Reardon et al., 2021; EPTA Collaboration et al., 2023a). The fixed values that have been used for the solar wind amplitude vary significantly. For instance, the default value used in the tempo software package is 10 cm^-3, while in the tempo2 software package it is 4 cm^-3. Using the NANOGrav 11-Year dataset, Madison et al. (2019) suggested an average value of 7.9 cm^-3 among all pulsars. However, the solar-wind density at 1-AU distance is seen to be both time and ecliptic latitude dependent (Sokół et al., 2013; Provornikova et al., 2014; Tiburzi et al., 2021; Susarla et al., 2024), so applying the same $n_{\rm sw}$ value to different pulsars and datasets is probably not optimal. Moreover, as we point out later, if a fixed constant $n_{\rm sw}$ is used which significantly differs from the actual time-averaged value, it can introduce a significant bias in other timing parameters that are strongly correlated, such as the astrometric parameters. The correlation between the timing parallax and solar-wind parameters was already briefly discussed in Splaver et al. (2005).

The paper is organized as follows. In Section 2 we lay out mathematical framework for understanding the correlation between the astrometric and solar-wind effects. In Section 3 we present a detailed investigation into the impact on astrometric measurements by the solar-wind effect, using mock-data simulations across various parameter spaces. The results from reprocessing the EPTA ‘DR2full’ dataset guided by our findings is shown in Section 4. We conclude in Section 5 with a brief discussion.

The core of pulsar timing analysis is the time-transformation formula (or the timing model) that is a linear combination of a sequence of time-delay terms. Each term includes a collection of timing parameters to be determined during the timing analysis. The top-level timing model can be expressed as:

where $T_{\rm p}$ is the pulsar proper time, $t_{\rm topo}$ is the topocentric time. Here $\Delta_{\rm R\odot}$, $\Delta_{\rm p}$, $\Delta_{\rm sol}$ represent the time delay introduced by Earth motion in the solar system (the solar-system Roemer delay), annual parallax of the pulsar and free electron from the solar wind. For simplicity, we only explicitly write the terms directly related to the study in this paper. More details of the pulsar timing model can be found in e.g., Edwards et al. (2006).

The solar-system Roemer delay, $\Delta_{\rm R\odot}$, can be written as

where $\vec{r}(E)$ is the Sun–Earth vector, $\theta(E)$ is the angle between Sun–Pulsar and Sun–Earth vectors, $c$ is the speed of light, and $E$ is the eccentric anomaly of the Earth at a given epoch. Following Keplerian’s law, $\vec{r}(E)$ and $\theta(E)$ can be further expressed by:

where $a$, $e$, $A_{\rm T}$ are the semi-major axis, eccentricity and true anomaly of the Earth orbit, $\lambda$ and $\beta$ are the pulsar longitude and latitude in (geocentric) ecliptic coordinate.

The parallax delay, $\Delta_{\rm p}$, represents the curvature of wavefront of the pulsar signal caused by the finite distance to the pulsar. It is given by (e.g., Edwards et al., 2006):

Here, the parallax $\varpi$ is in unit of arcseconds and $E$ is the eccentric anomaly of the Earth’s orbit.

To first order, the delay introduced by the solar wind, $\Delta_{\rm sol}$, can be approximated by assuming a spherically symmetric free-electron distribution around the Sun. This delay is written as (e.g., Issautier et al., 1998; Edwards et al., 2006):

where ${\rm D}_{\rm 0}=2.410\times 10^{-16}\rm cm^{-3}~{}pc$ and $f$ is the observing radio frequency. The only model parameter here, $n_{\rm sw}$, is the mean electron density at 1-AU heliocentric radius. Though in reality, the solar wind exhibits complex time variability (Sokół et al., 2013; Provornikova et al., 2014; Shaifullah et al., 2020; Tiburzi et al., 2021; Susarla et al., 2024, e.g.,), this model is still a good approximation for examining the correlation of solar-wind effect with other timing parameters, and will be used in the subsequent part of the paper.

Clearly, all three delay terms above exhibit annual or bi-annual periodic variation modulated by $\theta(E)$, so some of their model parameters are likely to show some significant correlation. Figure 1 gives the waveform of $\Delta_{\rm p}$ and $\Delta_{\rm sol}$ over the course of a year, for various ecliptic latitudes of the pulsar. It can be seen that the two waveforms follow anti-correlated trends, which is most apparent when the amplitude of the solar-wind effect is around its maximum and the parallax delay reaches a minimum. The correlation coefficient between the waveforms of these two delay terms can be directly calculated via:

where

and

The overline in Eq. (7) represents time-averaging over the course of a year. For a given function of the eccentric anomaly of the Earth’s orbit $E$, its annual average is calculated by

As shown in Eq. (1), the timing model is a linear combination of a sequence of delay terms including $\Delta_{\rm p}$, $\Delta_{\rm sol}$. Note that the two model parameters, $\varpi$ and $n_{\rm sw}$, are simply the amplitude of the two terms, $\rho_{\Delta_{\rm p},\Delta_{\rm sol}}$ should directly reflect the correlation of the two parameters in a least-square fit to the timing model. Meanwhile, since $\Delta_{\rm p}$ and $\Delta_{\rm sol}$ follow an anti-correlated trend, an increase in $\Delta_{\rm p}$ (thus in $\varpi$) requires a decrease in $\Delta_{\rm sol}$ (thus increase in $n_{\rm sw}$) to keep $\chi^{2}$ around its minimum in a least-square fit. This means that $\varpi$ and $n_{\rm sw}$ are positively correlated, i.e., $\rho_{\varpi,n_{\rm sw}}=-\rho_{\Delta_{\rm p},\Delta_{\rm sol}}$. Thus, it is clear from Eq. (7) that $\rho_{\varpi,n_{\rm sw}}$ does not depend on the value of $n_{\rm sw}$ and $\varpi$ or the observing frequency, but on the waveform of the two delay terms over the course of the year. Both delay terms are functions of the ecliptic latitude $\beta$ of the pulsar, with their modulation increasing as $\beta$ becomes smaller. Hence, $\rho_{\varpi,n_{\rm sw}}$ is likely to exhibit a significant dependency on $\beta$.

Note that the correlation estimate above only takes into account $\Delta_{\rm p}$ and $\Delta_{\rm sol}$. In practice, during the timing analysis, a global fit is usually carried out for all model parameters. Since other parameters, such as the position of the pulsar, could also show a significant correlation with $\varpi$ and $n_{\rm sw}$, the $\rho_{\varpi,n_{\rm sw}}$ calculated from the timing analysis can be significantly different. Nonetheless, the correlation analysis can also be conducted through the covariance matrix in the least-square timing fit (Blandford & Teukolsky, 1976), and is detailed in Appendix A. Adding the two position parameters, ecliptic longitude $\lambda$ and latitude $\beta$ of the pulsar, leads to the expression shown in Eq. (25). The derivation including only $\varpi$ and $n_{\rm sw}$ gives Eq. (29), exactly the same as in Eq. (7).

Figure 2 shows the correlation coefficients $\rho_{\varpi,n_{\rm sw}}$, $\rho_{\lambda,n_{\rm sw}}$, $\rho_{\beta,n_{\rm sw}}$ calculated from Eq. (25)–(27) and Eq. (29), for different ecliptic latitudes of the pulsar. It can be seen that when the position parameters are also included in the fit, $\rho_{\varpi,n_{\rm sw}}$ is indeed significantly altered; it becomes larger for all $\beta$ values shown in the figure. The correlation between $\beta$ and $n_{\rm sw}$ exhibits a similar trend, and is already significant when $\beta$ is more than a few degree, while for $\lambda$ it is significant only when $\beta$ approaches 70 deg and above.

There are several additional factors that could influence the correlations between astrometric and solar-wind parameters when fitting to real data. Firstly, the estimates above assume single-frequency observations, whereas in reality, multi-frequency coverage could mitigate these correlations since the solar-wind effect is chromatic while the others are not. Secondly, the calculations above assume complete sampling of the time-delay features over the course of a year. In practice, however, for low ecliptic latitude pulsars where the solar wind time delay changes rapidly around solar conjunction, a typical cadence (weekly to monthly) of timing observations is likely insufficient to fully resolve the solar-wind feature, which may require daily cadence. Additionally, observing the pulsar when it is very close to the Sun in the sky is difficult. For most single-dish radio telescopes, the reflection and concentration of optical and infrared radiation from the Sun can drastically heat up the receiver, and the convergence of radio emission from the Sun can strongly increase the system temperature of the observation. The resulting lack of coverage around solar conjunction could further impact the correlations between astrometric and solar-wind parameters in the fit to real data. All these aspects will be discussed in the following sections.

In this section, we conduct a series of simulation studies to investigate the correlation between astrometric and solar-wind parameters, and their measurability under various circumstances. We utilized the fake plugin in the tempo2 software package to generate mock pulsar timing data, i.e., the TOAs. These simulations encompass a broad range of ecliptic latitudes to reflect the actual positional distribution of pulsars in the sky (Hobbs et al., 2004). In addition to astrometric and solar-wind parameters, the pulsar ephemeris used in the creation of our mock dataset also encompasses a set of other standard timing parameters which are the same in all simulations. This includes the pulsar’s spin frequency and its derivative, the dispersion measure (the integrated column density of free electrons between the pulsar and Earth) along with its first and second derivatives–whenever multi-frequency data were simulated–and the pulsar’s proper motion in relation to the solar system. All these parameters were fitted simultaneously to the simulated data during the timing analysis. The corresponding covariance matrix of the fit was calculated by the tempo2 software package which follows the same framework as used in Appendix A. Note that more details of the simulated data will be presented in each individual subsections below, as they can be different for various simulation studies.

Following the findings from simulation studies presented above, in this section we reprocessed data from the second data release of the European Pulsar Timing Array (EPTA), to incorporate modelling of $n_{\rm sw}$ in the timing analysis. We used the version ‘DR2full‘ dataset (see Appendix A of EPTA Collaboration et al., 2023a) that contains data from both the legacy and new-generation pulsar backends in the EPTA, mostly at gigahertz frequencies. The analysis focused on a selection of pulsars where $\varpi$ was measured and the anticipated scale of the solar-wind delay is comparable to the achieved timing precision, allowing for reasonable constraints on $n_{\rm sw}$. The deterministic timing model described in EPTA Collaboration et al. (2023a) were employed for the timing analysis. To account for the red noise presented in the data, we used the customized noise model determined by EPTA Collaboration et al. (2023b) for each pulsar. In short, three types of red-noise model were considered to be included: the chromatic red noise caused by DM and scattering variation, respectively, and the achromatic red-noise component. In all models, the red signal was assumed to be a stationary, stochastic process with a power-law spectrum in the form of

where $S(\mathscr{f})$, $A$, $\gamma$, $\mathscr{f}_{\rm r}$ are the power spectral density as a function of Fourier frequency $\mathscr{f}$, spectral amplitude, spectral index and the reference Fourier frequency (set to 1 yr^-1), respectively. Here, the spectral amplitude relates to the observing radio frequency as $A\propto\nu^{-\alpha}$, where the chromatic index $\alpha$ is $2,4,0$ for the case of DM variation, scattering variation and achromatic noise, respectively. The constant $C_{\rm 0}$ equals to $12\pi^{2}$ for scattering variation and achromatic red noise, and is $k^{2}_{\rm DM}$ for DM variation where $k_{\rm DM}=2.41\times 10^{-4}$ cm^-3 pc MHz² s^-1 is the DM constant. Again, the temponest software package was utilized to perform the Bayesian timing analysis and generate the posterior distributions of the model parameters (Lentati et al., 2014). Both $\varpi$ and $n_{\rm sw}$ were sampled simultaneously with the other timing and noise model parameters. The results of parameter measurements are summarized in Table 1, with the joint posterior distributions of $\varpi$ and $n_{\rm sw}$ presented in Figure 15.

In the analysis of EPTA Collaboration et al. (2023a), $n_{\rm sw}$ was fixed to 7.9 cm^-3 for all pulsars (Madison et al., 2019), except for PSR J1022+1001 where a fitted value of $11.1\pm 0.3$ cm^-3 was obtained. For PSRs J0030+0451, J0751+1807, J1730$-$2304, our timing analysis yields values very close to the previously fixed $n_{\rm sw}$. The measured positions and parallaxes are highly consistent with those from EPTA Collaboration et al. (2023a), though of slightly lower precision. This outcome is entirely expected due to the inclusion of $n_{\rm sw}$ in the timing fit and the resulting correlation between $\varpi$, $n_{\rm sw}$.

Meanwhile, for PSRs J1600$-$3053, J1713+0747, J1909$-$3744 and J1744$-$1134, the fitted $\varpi$ and $\beta$ values from our analysis are smaller than those obtained in EPTA Collaboration et al. (2023a), and their fitted $n_{\rm sw}$ values are also correspondingly less than 7.9 cm^-3 as anticipated from our simulations. In fact, the inclusion of $n_{\rm sw}$ in the timing fit delivers parallax measurements more consistent with those reported elsewhere, such as in Agazie et al. (2023) where the ‘DMX‘ model was employed to capture chromatic effects in the data including both DM and solar wind, and Reardon et al. (2021) where a fixed $n_{\rm sw}$ of 4 cm^-3 was used in the timing analysis. For instance, in the case of PSR J1713+0747, the fitted value of $n_{\rm sw}$ is only marginally constrained at $4.0\pm 1.3$ cm^-3. At the same time, the fitted $\varpi$ value now decreases from $0.88\pm 0.01$ to $0.83\pm 0.01$ mas, which is highly consistent with $0.833\pm 0.015$ mas in Agazie et al. (2023) and closer to $0.763\pm 0.021$ mas in (Reardon et al., 2021). The decrease of 0.05 mas generally matches the prediction in Figure 6 for a -3.9 cm^-3 offset in $n_{\rm sw}$. In all cases of Table 1, the fitted values of ecliptic longitude from our analysis are the same as those obtained in EPTA Collaboration et al. (2023a).

The relation between our obtained $n_{\rm sw}$ values to the ecliptic latitude of the pulsar matches the expectations from Susarla et al. (2024). For pulsars at very low ecliptic latitudes ($\beta\lesssim 5$ deg: PSRs J0030+0451, J0751+1807, J1730$-$2304) the solar-wind density is higher because the slow (and denser) solar wind is the dominant component along the line-of-sight (LOS). For those at higher ecliptic longitudes (here PSRs J1600$-$3053, J1713+0747, J1744$-$1134 and J1909$-$3744), the value is lower since the LOS could pass through either primarily the fast solar wind (around solar minima) or a mixture of the two components (around solar maxima). For PSRs J0030+0451, J1730$-$2304 and J1744$-$1134, our measured values from EPTA ‘DR2full‘ dataset show good consistency with those reported by Susarla et al. (2024) using low-frequency ($\sim 200$ MHz) observations with LOFAR. Nonetheless, the interpretation of these measurements should be taken with caution as pointed out later.

The noise-model parameters measured from our analysis are highly consistent with those obtained in EPTA Collaboration et al. (2023b). Notably, the achromatic red-noise measurements remain unchanged after the inclusion of $n_{\rm sw}$ in the timing fit. Therefore, we do not expect significant impact on the gravitational-wave search results presented in EPTA Collaboration et al. (2023c) by fixing the $n_{\rm sw}$ value.

In this paper, we investigated the impact on measurement of astrometry in pulsar timing by the solar wind effect, which could pose a potential issue for precisely measuring pulsar astrometry and distance. Using both theoretical calculation and mock-data simulations, we demonstrated a significant correlation between the pulsar position and parallax parameters, and the solar-wind density parameter, $n_{\rm sw}$. This correlation strongly depends on the ecliptic latitude of the pulsar. We showed that using a fixed $n_{\rm sw}$ with an arbitrary value in the timing analysis could introduce significant bias in the estimated pulsar position and parallax, leading to inaccuracies in the pulsar distance measurement. The significance of the bias strongly depends on the ecliptic latitude of the pulsar and the timing precision of the dataset. For pulsars with favourable timing precision and ecliptic latitude, the astrometric and solar-wind parameters can be measured jointly with other timing parameters already using single-frequency data, which will avoid the bias in the estimate of astrometric parameters. The correlation between astrometric and solar-wind parameters can be mitigated by using data with multi-frequency coverage, which also greatly improves the measurement precision of these parameters. This is in particular significant for pulsars at a medium or high ecliptic latitude ($\beta\gtrsim 30$ deg). Furthermore, we reprocessed EPTA ‘DR2full‘ dataset to include the solar-wind effect in the timing fit for a selection of pulsars, which resulted in significant measurements of both parallax and solar-wind density. For pulsars where $n_{\rm sw}$ differed from the fixed value of 7.9 cm^-3 in the EPTA DR2 analysis, the parallax measurement became more consistent with values reported elsewhere. For sources overlapping with the analysis in Susarla et al. (2024) using low-frequency observations, our $n_{\rm sw}$ measurements align well with theirs.

As pointed out by various previous work, the spherical model in Eq. 6 is not complete in describing the solar-wind effect. Most importantly, it neglects the time variation in the solar-wind density which can be caused by both the long-term (11-year) solar cycles and short-term stochastic process in the fast solar wind component (Shaifullah et al., 2020; Tiburzi et al., 2021). For instance, recently Susarla et al. (2024) observed complex time variability in the solar-wind intensity for several PTA pulsars using data collected with LOFAR. Therefore, our measurements using the model in Eq. (6) should be seen as representing a time average of the solar-wind density over many years. Caution needs to be taken when using these measurements to derive other physical properties of the solar wind. Nonetheless, it was shown that the EPTA DR2 data, mostly at frequencies higher than 1 GHz, is likely not to be sensitive to such time variation (Niţu et al., 2024). The model should also serve the primary purpose of our simulation studies to understand the correlation between astrometric and solar-wind parameters. It is anticipated that being able to model the solar-wind time variability could mitigate the degeneracy between astrometric and solar-wind parameters.

Our simulations have shown that the measurability of $\varpi$ and $n_{\rm sw}$ strongly depends on several variables, including the pulsar’s ecliptic latitude, observing frequency, timing precision of the data, etc.. At 1.4 GHz, the typical frequency of most precision timing observations, for a pulsar at a very low ecliptic latitude of $\beta\lesssim 15$ deg with a timing precision of the order 1 $\upmu$s, it should always be recommended to include solar-wind modelling in the timing analysis to avoid significant bias in the fitted astrometric parameters. This also applies to pulsars at higher ecliptic latitude of 15 deg$\lesssim\beta\lesssim 30$ deg, if the timing precision is of the order 100 ns. For pulsars with $\beta\gtrsim 30$ deg, it is possible at this frequency to encounter a situation where $\varpi$ is well measured with the data but $n_{\rm sw}$ may not. In this scenario, while fitting for $n_{\rm sw}$ still remains a conservative option, using a given $n_{\rm sw}$ may also be feasible, which may introduce a bias in the astrometric parameters that is tolerable within the measurement uncertainty. If this approach is employed in the timing analysis, based on the results from Susarla et al. (2024), a fixed value (e.g., 4 cm^-3) lower than the gross average of many PTA pulsars (7.9 cm^-3) obtained in Madison et al. (2019), would likely provide a more realistic representation of the actual $n_{\rm sw}$ for pulsars falling in this ecliptic latitude range.

Our work shows complexity in obtaining robust astrometry measurements in pulsar timing analysis, due to the significant correlation between the astrometric and solar-wind parameters. In addition, Madison et al. (2013) has also demonstrated that timing measurement of astrometry can be contaminated by red noise in the timing data, if it is not properly accounted for. Therefore, to obtain unbiased astrometric measurements with pulsar timing, it would be recommended to conduct both the timing and noise analysis at the same time. Meanwhile, measurements from independent experiments, such as the VLBI observations (e.g., Ding et al., 2023) or the scintillation studies (e.g., Reardon et al., 2019), could be used to verify the timing results and even help to break the degeneracy among the model parameters.

In line with previous studies (e.g., Lee et al., 2012; Lam et al., 2018; Donner et al., 2020), our results highlight the importance of achieving high timing precision at higher radio frequencies in pulsar timing analysis to correct for chromatic effects, such as those caused by the solar wind. However, this precision is not easily achieved and requires dedicated efforts. Since pulsars are typically steep-spectrum radio sources (e.g., Jankowski et al., 2018), their flux density is significantly higher at lower frequencies, making it easier to achieve sufficient timing precision for modeling chromatic effects without placing a strong demand on system sensitivity. In contrast, at higher radio frequencies, the sharp drop in flux density necessitates the use of the largest radio telescopes or extended integration times to achieve comparative timing precision. Despite these challenges, using the most sensitive telescopes for high-frequency timing observations aligns with recommendations in other studies. Notably, timing precision at canonical timing frequencies (e.g., 1.4 GHz) for many MSPs observed with the largest radio telescopes is often limited by intrinsic white noise caused by the instability of pulsar emission on short timescales, known as “jitter noise” (e.g., Liu et al., 2011; Parthasarathy et al., 2021; Wang et al., 2024). In such cases, conducting timing observations at higher frequencies is more beneficial, as it allows for a precision that less sensitive telescopes cannot achieve. Therefore, using more sensitive telescopes at higher radio frequencies is advantageous for multiple reasons.

Finally, our simulation studies have focused on assessing the impact of the solar-wind effect on astrometric measurements in precision pulsar timing, particularly under the assumption of canonical timing observations around gigahertz radio frequencies with weekly cadence. We emphasise the significant value of ongoing low-frequency ($<800$ MHz) pulsar timing programs, such as those conducted with the Canadian Hydrogen Intensity Mapping Experiment (CHIME) Telescope (CHIME/Pulsar Collaboration et al., 2021) and LOFAR (Tiburzi et al., 2019), for studying chromatic effects in timing analysis such as the solar wind. Meanwhile, due to limitations in achievable timing precision and baseline, these low-frequency data alone may not suffice for measuring pulsar astrometry or achieving competitive precision. Therefore, combining low-frequency and canonical pulsar timing data would provide the optimal solution for simultaneously and robustly measuring both pulsar astrometry and the solar-wind properties.

We thank Joris Verbiest and Paulo Freire for reading the manuscript and providing constructive comments. KL and MK acknowledge funding by the MPG-CAS LEGACY programme. JA acknowledges support from the European Commission (ARGOS-CDS; Grant Agreement number: 101094354). A.C. acknowledges financial support provided under the European Union’s Horizon Europe ERC Starting Grant "A Gamma-ray Infrastructure to Advance Gravitational Wave Astrophysics" (GIGA, Grant Agreement: 101116134). The Effelsberg 100-m telescope is operated by the Max-Planck-Institut für Radioastronomie. The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from “Programme National de Cosmologie and Galaxies” (PNCG), and “Programme National Hautes Energies” (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radioastronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK’s Science and Technology Facilities Council (STFC).

The EPTA ‘DR2full‘ dataset are available at: https://zenodo.org/records/8300645.